#
source:
trunk/third/gmp/tests/rand/statlib.c
@
18191

Revision 18191, 20.3 KB checked in by ghudson, 22 years ago (diff) |
---|

Line | |
---|---|

1 | /* statlib.c -- Statistical functions for testing the randomness of |

2 | number sequences. */ |

3 | |

4 | /* |

5 | Copyright 1999, 2000 Free Software Foundation, Inc. |

6 | |

7 | This file is part of the GNU MP Library. |

8 | |

9 | The GNU MP Library is free software; you can redistribute it and/or modify |

10 | it under the terms of the GNU Lesser General Public License as published by |

11 | the Free Software Foundation; either version 2.1 of the License, or (at your |

12 | option) any later version. |

13 | |

14 | The GNU MP Library is distributed in the hope that it will be useful, but |

15 | WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY |

16 | or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public |

17 | License for more details. |

18 | |

19 | You should have received a copy of the GNU Lesser General Public License |

20 | along with the GNU MP Library; see the file COPYING.LIB. If not, write to |

21 | the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, |

22 | MA 02111-1307, USA. |

23 | */ |

24 | |

25 | /* The theories for these functions are taken from D. Knuth's "The Art |

26 | of Computer Programming: Volume 2, Seminumerical Algorithms", Third |

27 | Edition, Addison Wesley, 1998. */ |

28 | |

29 | /* Implementation notes. |

30 | |

31 | The Kolmogorov-Smirnov test. |

32 | |

33 | Eq. (13) in Knuth, p. 50, says that if X1, X2, ..., Xn are independent |

34 | observations arranged into ascending order |

35 | |

36 | Kp = sqr(n) * max(j/n - F(Xj)) for all 1<=j<=n |

37 | Km = sqr(n) * max(F(Xj) - (j-1)/n)) for all 1<=j<=n |

38 | |

39 | where F(x) = Pr(X <= x) = probability that (X <= x), which for a |

40 | uniformly distributed random real number between zero and one is |

41 | exactly the number itself (x). |

42 | |

43 | |

44 | The answer to exercise 23 gives the following implementation, which |

45 | doesn't need the observations to be sorted in ascending order: |

46 | |

47 | for (k = 0; k < m; k++) |

48 | a[k] = 1.0 |

49 | b[k] = 0.0 |

50 | c[k] = 0 |

51 | |

52 | for (each observation Xj) |

53 | Y = F(Xj) |

54 | k = floor (m * Y) |

55 | a[k] = min (a[k], Y) |

56 | b[k] = max (b[k], Y) |

57 | c[k] += 1 |

58 | |

59 | j = 0 |

60 | rp = rm = 0 |

61 | for (k = 0; k < m; k++) |

62 | if (c[k] > 0) |

63 | rm = max (rm, a[k] - j/n) |

64 | j += c[k] |

65 | rp = max (rp, j/n - b[k]) |

66 | |

67 | Kp = sqr (n) * rp |

68 | Km = sqr (n) * rm |

69 | |

70 | */ |

71 | |

72 | #include <stdio.h> |

73 | #include <stdlib.h> |

74 | #include <math.h> |

75 | |

76 | #include "gmp.h" |

77 | #include "gmpstat.h" |

78 | |

79 | /* ks (Kp, Km, X, P, n) -- Perform a Kolmogorov-Smirnov test on the N |

80 | real numbers between zero and one in vector X. P is the |

81 | distribution function, called for each entry in X, which should |

82 | calculate the probability of X being greater than or equal to any |

83 | number in the sequence. (For a uniformly distributed sequence of |

84 | real numbers between zero and one, this is simply equal to X.) The |

85 | result is put in Kp and Km. */ |

86 | |

87 | void |

88 | ks (mpf_t Kp, |

89 | mpf_t Km, |

90 | mpf_t X[], |

91 | void (P) (mpf_t, mpf_t), |

92 | unsigned long int n) |

93 | { |

94 | mpf_t Kt; /* temp */ |

95 | mpf_t f_x; |

96 | mpf_t f_j; /* j */ |

97 | mpf_t f_jnq; /* j/n or (j-1)/n */ |

98 | unsigned long int j; |

99 | |

100 | /* Sort the vector in ascending order. */ |

101 | qsort (X, n, sizeof (__mpf_struct), mpf_cmp); |

102 | |

103 | /* K-S test. */ |

104 | /* Kp = sqr(n) * max(j/n - F(Xj)) for all 1<=j<=n |

105 | Km = sqr(n) * max(F(Xj) - (j-1)/n)) for all 1<=j<=n |

106 | */ |

107 | |

108 | mpf_init (Kt); mpf_init (f_x); mpf_init (f_j); mpf_init (f_jnq); |

109 | mpf_set_ui (Kp, 0); mpf_set_ui (Km, 0); |

110 | for (j = 1; j <= n; j++) |

111 | { |

112 | P (f_x, X[j-1]); |

113 | mpf_set_ui (f_j, j); |

114 | |

115 | mpf_div_ui (f_jnq, f_j, n); |

116 | mpf_sub (Kt, f_jnq, f_x); |

117 | if (mpf_cmp (Kt, Kp) > 0) |

118 | mpf_set (Kp, Kt); |

119 | if (g_debug > DEBUG_2) |

120 | { |

121 | printf ("j=%lu ", j); |

122 | printf ("P()="); mpf_out_str (stdout, 10, 2, f_x); printf ("\t"); |

123 | |

124 | printf ("jnq="); mpf_out_str (stdout, 10, 2, f_jnq); printf (" "); |

125 | printf ("diff="); mpf_out_str (stdout, 10, 2, Kt); printf (" "); |

126 | printf ("Kp="); mpf_out_str (stdout, 10, 2, Kp); printf ("\t"); |

127 | } |

128 | mpf_sub_ui (f_j, f_j, 1); |

129 | mpf_div_ui (f_jnq, f_j, n); |

130 | mpf_sub (Kt, f_x, f_jnq); |

131 | if (mpf_cmp (Kt, Km) > 0) |

132 | mpf_set (Km, Kt); |

133 | |

134 | if (g_debug > DEBUG_2) |

135 | { |

136 | printf ("jnq="); mpf_out_str (stdout, 10, 2, f_jnq); printf (" "); |

137 | printf ("diff="); mpf_out_str (stdout, 10, 2, Kt); printf (" "); |

138 | printf ("Km="); mpf_out_str (stdout, 10, 2, Km); printf (" "); |

139 | printf ("\n"); |

140 | } |

141 | } |

142 | mpf_sqrt_ui (Kt, n); |

143 | mpf_mul (Kp, Kp, Kt); |

144 | mpf_mul (Km, Km, Kt); |

145 | |

146 | mpf_clear (Kt); mpf_clear (f_x); mpf_clear (f_j); mpf_clear (f_jnq); |

147 | } |

148 | |

149 | /* ks_table(val, n) -- calculate probability for Kp/Km less than or |

150 | equal to VAL with N observations. See [Knuth section 3.3.1] */ |

151 | |

152 | void |

153 | ks_table (mpf_t p, mpf_t val, const unsigned int n) |

154 | { |

155 | /* We use Eq. (27), Knuth p.58, skipping O(1/n) for simplicity. |

156 | This shortcut will result in too high probabilities, especially |

157 | when n is small. |

158 | |

159 | Pr(Kp(n) <= s) = 1 - pow(e, -2*s^2) * (1 - 2/3*s/sqrt(n) + O(1/n)) */ |

160 | |

161 | /* We have 's' in variable VAL and store the result in P. */ |

162 | |

163 | mpf_t t1, t2; |

164 | |

165 | mpf_init (t1); mpf_init (t2); |

166 | |

167 | /* t1 = 1 - 2/3 * s/sqrt(n) */ |

168 | mpf_sqrt_ui (t1, n); |

169 | mpf_div (t1, val, t1); |

170 | mpf_mul_ui (t1, t1, 2); |

171 | mpf_div_ui (t1, t1, 3); |

172 | mpf_ui_sub (t1, 1, t1); |

173 | |

174 | /* t2 = pow(e, -2*s^2) */ |

175 | #ifndef OLDGMP |

176 | mpf_pow_ui (t2, val, 2); /* t2 = s^2 */ |

177 | mpf_set_d (t2, exp (-(2.0 * mpf_get_d (t2)))); |

178 | #else |

179 | /* hmmm, gmp doesn't have pow() for floats. use doubles. */ |

180 | mpf_set_d (t2, pow (M_E, -(2 * pow (mpf_get_d (val), 2)))); |

181 | #endif |

182 | |

183 | /* p = 1 - t1 * t2 */ |

184 | mpf_mul (t1, t1, t2); |

185 | mpf_ui_sub (p, 1, t1); |

186 | |

187 | mpf_clear (t1); mpf_clear (t2); |

188 | } |

189 | |

190 | static double x2_table_X[][7] = { |

191 | { -2.33, -1.64, -.674, 0.0, 0.674, 1.64, 2.33 }, /* x */ |

192 | { 5.4289, 2.6896, .454276, 0.0, .454276, 2.6896, 5.4289} /* x^2 */ |

193 | }; |

194 | |

195 | #define _2D3 ((double) .6666666666) |

196 | |

197 | /* x2_table (t, v, n) -- return chi-square table row for V in T[]. */ |

198 | void |

199 | x2_table (double t[], |

200 | unsigned int v) |

201 | { |

202 | int f; |

203 | |

204 | |

205 | /* FIXME: Do a table lookup for v <= 30 since the following formula |

206 | [Knuth, vol 2, 3.3.1] is only good for v > 30. */ |

207 | |

208 | /* value = v + sqrt(2*v) * X[p] + (2/3) * X[p]^2 - 2/3 + O(1/sqrt(t) */ |

209 | /* NOTE: The O() term is ignored for simplicity. */ |

210 | |

211 | for (f = 0; f < 7; f++) |

212 | t[f] = |

213 | v + |

214 | sqrt (2 * v) * x2_table_X[0][f] + |

215 | _2D3 * x2_table_X[1][f] - _2D3; |

216 | } |

217 | |

218 | |

219 | /* P(p, x) -- Distribution function. Calculate the probability of X |

220 | being greater than or equal to any number in the sequence. For a |

221 | random real number between zero and one given by a uniformly |

222 | distributed random number generator, this is simply equal to X. */ |

223 | |

224 | static void |

225 | P (mpf_t p, mpf_t x) |

226 | { |

227 | mpf_set (p, x); |

228 | } |

229 | |

230 | /* mpf_freqt() -- Frequency test using KS on N real numbers between zero |

231 | and one. See [Knuth vol 2, p.61]. */ |

232 | void |

233 | mpf_freqt (mpf_t Kp, |

234 | mpf_t Km, |

235 | mpf_t X[], |

236 | const unsigned long int n) |

237 | { |

238 | ks (Kp, Km, X, P, n); |

239 | } |

240 | |

241 | |

242 | /* The Chi-square test. Eq. (8) in Knuth vol. 2 says that if Y[] |

243 | holds the observations and p[] is the probability for.. (to be |

244 | continued!) |

245 | |

246 | V = 1/n * sum((s=1 to k) Y[s]^2 / p[s]) - n */ |

247 | |

248 | void |

249 | x2 (mpf_t V, /* result */ |

250 | unsigned long int X[], /* data */ |

251 | unsigned int k, /* #of categories */ |

252 | void (P) (mpf_t, unsigned long int, void *), /* probability func */ |

253 | void *x, /* extra user data passed to P() */ |

254 | unsigned long int n) /* #of samples */ |

255 | { |

256 | unsigned int f; |

257 | mpf_t f_t, f_t2; /* temp floats */ |

258 | |

259 | mpf_init (f_t); mpf_init (f_t2); |

260 | |

261 | |

262 | mpf_set_ui (V, 0); |

263 | for (f = 0; f < k; f++) |

264 | { |

265 | if (g_debug > DEBUG_2) |

266 | fprintf (stderr, "%u: P()=", f); |

267 | mpf_set_ui (f_t, X[f]); |

268 | mpf_mul (f_t, f_t, f_t); /* f_t = X[f]^2 */ |

269 | P (f_t2, f, x); /* f_t2 = Pr(f) */ |

270 | if (g_debug > DEBUG_2) |

271 | mpf_out_str (stderr, 10, 2, f_t2); |

272 | mpf_div (f_t, f_t, f_t2); |

273 | mpf_add (V, V, f_t); |

274 | if (g_debug > DEBUG_2) |

275 | { |

276 | fprintf (stderr, "\tV="); |

277 | mpf_out_str (stderr, 10, 2, V); |

278 | fprintf (stderr, "\t"); |

279 | } |

280 | } |

281 | if (g_debug > DEBUG_2) |

282 | fprintf (stderr, "\n"); |

283 | mpf_div_ui (V, V, n); |

284 | mpf_sub_ui (V, V, n); |

285 | |

286 | mpf_clear (f_t); mpf_clear (f_t2); |

287 | } |

288 | |

289 | /* Pzf(p, s, x) -- Probability for category S in mpz_freqt(). It's |

290 | 1/d for all S. X is a pointer to an unsigned int holding 'd'. */ |

291 | static void |

292 | Pzf (mpf_t p, unsigned long int s, void *x) |

293 | { |

294 | mpf_set_ui (p, 1); |

295 | mpf_div_ui (p, p, *((unsigned int *) x)); |

296 | } |

297 | |

298 | /* mpz_freqt(V, X, imax, n) -- Frequency test on integers. [Knuth, |

299 | vol 2, 3.3.2]. Keep IMAX low on this one, since we loop from 0 to |

300 | IMAX. 128 or 256 could be nice. |

301 | |

302 | X[] must not contain numbers outside the range 0 <= X <= IMAX. |

303 | |

304 | Return value is number of observations actally used, after |

305 | discarding entries out of range. |

306 | |

307 | Since X[] contains integers between zero and IMAX, inclusive, we |

308 | have IMAX+1 categories. |

309 | |

310 | Note that N should be at least 5*IMAX. Result is put in V and can |

311 | be compared to output from x2_table (v=IMAX). */ |

312 | |

313 | unsigned long int |

314 | mpz_freqt (mpf_t V, |

315 | mpz_t X[], |

316 | unsigned int imax, |

317 | const unsigned long int n) |

318 | { |

319 | unsigned long int *v; /* result */ |

320 | unsigned int f; |

321 | unsigned int d; /* number of categories = imax+1 */ |

322 | unsigned int uitemp; |

323 | unsigned long int usedn; |

324 | |

325 | |

326 | d = imax + 1; |

327 | |

328 | v = (unsigned long int *) calloc (imax + 1, sizeof (unsigned long int)); |

329 | if (NULL == v) |

330 | { |

331 | fprintf (stderr, "mpz_freqt(): out of memory\n"); |

332 | exit (1); |

333 | } |

334 | |

335 | /* count */ |

336 | usedn = n; /* actual number of observations */ |

337 | for (f = 0; f < n; f++) |

338 | { |

339 | uitemp = mpz_get_ui(X[f]); |

340 | if (uitemp > imax) /* sanity check */ |

341 | { |

342 | if (g_debug) |

343 | fprintf (stderr, "mpz_freqt(): warning: input insanity: %u, "\ |

344 | "ignored.\n", uitemp); |

345 | usedn--; |

346 | continue; |

347 | } |

348 | v[uitemp]++; |

349 | } |

350 | |

351 | if (g_debug > DEBUG_2) |

352 | { |

353 | fprintf (stderr, "counts:\n"); |

354 | for (f = 0; f <= imax; f++) |

355 | fprintf (stderr, "%u:\t%lu\n", f, v[f]); |

356 | } |

357 | |

358 | /* chi-square with k=imax+1 and P(x)=1/(imax+1) for all x.*/ |

359 | x2 (V, v, d, Pzf, (void *) &d, usedn); |

360 | |

361 | free (v); |

362 | return (usedn); |

363 | } |

364 | |

365 | /* debug dummy to drag in dump funcs */ |

366 | void |

367 | foo_debug () |

368 | { |

369 | if (0) |

370 | { |

371 | mpf_dump (0); |

372 | #ifndef OLDGMP |

373 | mpz_dump (0); |

374 | #endif |

375 | } |

376 | } |

377 | |

378 | /* merit (rop, t, v, m) -- calculate merit for spectral test result in |

379 | dimension T, see Knuth p. 105. BUGS: Only valid for 2 <= T <= |

380 | 6. */ |

381 | void |

382 | merit (mpf_t rop, unsigned int t, mpf_t v, mpz_t m) |

383 | { |

384 | int f; |

385 | mpf_t f_m, f_const, f_pi; |

386 | |

387 | mpf_init (f_m); |

388 | mpf_set_z (f_m, m); |

389 | mpf_init_set_d (f_const, M_PI); |

390 | mpf_init_set_d (f_pi, M_PI); |

391 | |

392 | switch (t) |

393 | { |

394 | case 2: /* PI */ |

395 | break; |

396 | case 3: /* PI * 4/3 */ |

397 | mpf_mul_ui (f_const, f_const, 4); |

398 | mpf_div_ui (f_const, f_const, 3); |

399 | break; |

400 | case 4: /* PI^2 * 1/2 */ |

401 | mpf_mul (f_const, f_const, f_pi); |

402 | mpf_div_ui (f_const, f_const, 2); |

403 | break; |

404 | case 5: /* PI^2 * 8/15 */ |

405 | mpf_mul (f_const, f_const, f_pi); |

406 | mpf_mul_ui (f_const, f_const, 8); |

407 | mpf_div_ui (f_const, f_const, 15); |

408 | break; |

409 | case 6: /* PI^3 * 1/6 */ |

410 | mpf_mul (f_const, f_const, f_pi); |

411 | mpf_mul (f_const, f_const, f_pi); |

412 | mpf_div_ui (f_const, f_const, 6); |

413 | break; |

414 | default: |

415 | fprintf (stderr, |

416 | "spect (merit): can't calculate merit for dimensions > 6\n"); |

417 | mpf_set_ui (f_const, 0); |

418 | break; |

419 | } |

420 | |

421 | /* rop = v^t */ |

422 | mpf_set (rop, v); |

423 | for (f = 1; f < t; f++) |

424 | mpf_mul (rop, rop, v); |

425 | mpf_mul (rop, rop, f_const); |

426 | mpf_div (rop, rop, f_m); |

427 | |

428 | mpf_clear (f_m); |

429 | mpf_clear (f_const); |

430 | mpf_clear (f_pi); |

431 | } |

432 | |

433 | double |

434 | merit_u (unsigned int t, mpf_t v, mpz_t m) |

435 | { |

436 | mpf_t rop; |

437 | double res; |

438 | |

439 | mpf_init (rop); |

440 | merit (rop, t, v, m); |

441 | res = mpf_get_d (rop); |

442 | mpf_clear (rop); |

443 | return res; |

444 | } |

445 | |

446 | /* f_floor (rop, op) -- Set rop = floor (op). */ |

447 | void |

448 | f_floor (mpf_t rop, mpf_t op) |

449 | { |

450 | mpz_t z; |

451 | |

452 | mpz_init (z); |

453 | |

454 | /* No mpf_floor(). Convert to mpz and back. */ |

455 | mpz_set_f (z, op); |

456 | mpf_set_z (rop, z); |

457 | |

458 | mpz_clear (z); |

459 | } |

460 | |

461 | |

462 | /* vz_dot (rop, v1, v2, nelem) -- compute dot product of z-vectors V1, |

463 | V2. N is number of elements in vectors V1 and V2. */ |

464 | |

465 | void |

466 | vz_dot (mpz_t rop, mpz_t V1[], mpz_t V2[], unsigned int n) |

467 | { |

468 | mpz_t t; |

469 | |

470 | mpz_init (t); |

471 | mpz_set_ui (rop, 0); |

472 | while (n--) |

473 | { |

474 | mpz_mul (t, V1[n], V2[n]); |

475 | mpz_add (rop, rop, t); |

476 | } |

477 | |

478 | mpz_clear (t); |

479 | } |

480 | |

481 | void |

482 | spectral_test (mpf_t rop[], unsigned int T, mpz_t a, mpz_t m) |

483 | { |

484 | /* Knuth "Seminumerical Algorithms, Third Edition", section 3.3.4 |

485 | (pp. 101-103). */ |

486 | |

487 | /* v[t] = min { sqrt (x[1]^2 + ... + x[t]^2) | |

488 | x[1] + a*x[2] + ... + pow (a, t-1) * x[t] is congruent to 0 (mod m) } */ |

489 | |

490 | |

491 | /* Variables. */ |

492 | unsigned int ui_t; |

493 | unsigned int ui_i, ui_j, ui_k, ui_l; |

494 | mpf_t f_tmp1, f_tmp2; |

495 | mpz_t tmp1, tmp2, tmp3; |

496 | mpz_t U[GMP_SPECT_MAXT][GMP_SPECT_MAXT], |

497 | V[GMP_SPECT_MAXT][GMP_SPECT_MAXT], |

498 | X[GMP_SPECT_MAXT], |

499 | Y[GMP_SPECT_MAXT], |

500 | Z[GMP_SPECT_MAXT]; |

501 | mpz_t h, hp, r, s, p, pp, q, u, v; |

502 | |

503 | /* GMP inits. */ |

504 | mpf_init (f_tmp1); |

505 | mpf_init (f_tmp2); |

506 | for (ui_i = 0; ui_i < GMP_SPECT_MAXT; ui_i++) |

507 | { |

508 | for (ui_j = 0; ui_j < GMP_SPECT_MAXT; ui_j++) |

509 | { |

510 | mpz_init_set_ui (U[ui_i][ui_j], 0); |

511 | mpz_init_set_ui (V[ui_i][ui_j], 0); |

512 | } |

513 | mpz_init_set_ui (X[ui_i], 0); |

514 | mpz_init_set_ui (Y[ui_i], 0); |

515 | mpz_init (Z[ui_i]); |

516 | } |

517 | mpz_init (tmp1); |

518 | mpz_init (tmp2); |

519 | mpz_init (tmp3); |

520 | mpz_init (h); |

521 | mpz_init (hp); |

522 | mpz_init (r); |

523 | mpz_init (s); |

524 | mpz_init (p); |

525 | mpz_init (pp); |

526 | mpz_init (q); |

527 | mpz_init (u); |

528 | mpz_init (v); |

529 | |

530 | /* Implementation inits. */ |

531 | if (T > GMP_SPECT_MAXT) |

532 | T = GMP_SPECT_MAXT; /* FIXME: Lazy. */ |

533 | |

534 | /* S1 [Initialize.] */ |

535 | ui_t = 2 - 1; /* NOTE: `t' in description == ui_t + 1 |

536 | for easy indexing */ |

537 | mpz_set (h, a); |

538 | mpz_set (hp, m); |

539 | mpz_set_ui (p, 1); |

540 | mpz_set_ui (pp, 0); |

541 | mpz_set (r, a); |

542 | mpz_pow_ui (s, a, 2); |

543 | mpz_add_ui (s, s, 1); /* s = 1 + a^2 */ |

544 | |

545 | /* S2 [Euclidean step.] */ |

546 | while (1) |

547 | { |

548 | if (g_debug > DEBUG_1) |

549 | { |

550 | mpz_mul (tmp1, h, pp); |

551 | mpz_mul (tmp2, hp, p); |

552 | mpz_sub (tmp1, tmp1, tmp2); |

553 | if (mpz_cmpabs (m, tmp1)) |

554 | { |

555 | printf ("***BUG***: h*pp - hp*p = "); |

556 | mpz_out_str (stdout, 10, tmp1); |

557 | printf ("\n"); |

558 | } |

559 | } |

560 | if (g_debug > DEBUG_2) |

561 | { |

562 | printf ("hp = "); |

563 | mpz_out_str (stdout, 10, hp); |

564 | printf ("\nh = "); |

565 | mpz_out_str (stdout, 10, h); |

566 | printf ("\n"); |

567 | fflush (stdout); |

568 | } |

569 | |

570 | if (mpz_sgn (h)) |

571 | mpz_tdiv_q (q, hp, h); /* q = floor(hp/h) */ |

572 | else |

573 | mpz_set_ui (q, 1); |

574 | |

575 | if (g_debug > DEBUG_2) |

576 | { |

577 | printf ("q = "); |

578 | mpz_out_str (stdout, 10, q); |

579 | printf ("\n"); |

580 | fflush (stdout); |

581 | } |

582 | |

583 | mpz_mul (tmp1, q, h); |

584 | mpz_sub (u, hp, tmp1); /* u = hp - q*h */ |

585 | |

586 | mpz_mul (tmp1, q, p); |

587 | mpz_sub (v, pp, tmp1); /* v = pp - q*p */ |

588 | |

589 | mpz_pow_ui (tmp1, u, 2); |

590 | mpz_pow_ui (tmp2, v, 2); |

591 | mpz_add (tmp1, tmp1, tmp2); |

592 | if (mpz_cmp (tmp1, s) < 0) |

593 | { |

594 | mpz_set (s, tmp1); /* s = u^2 + v^2 */ |

595 | mpz_set (hp, h); /* hp = h */ |

596 | mpz_set (h, u); /* h = u */ |

597 | mpz_set (pp, p); /* pp = p */ |

598 | mpz_set (p, v); /* p = v */ |

599 | } |

600 | else |

601 | break; |

602 | } |

603 | |

604 | /* S3 [Compute v2.] */ |

605 | mpz_sub (u, u, h); |

606 | mpz_sub (v, v, p); |

607 | |

608 | mpz_pow_ui (tmp1, u, 2); |

609 | mpz_pow_ui (tmp2, v, 2); |

610 | mpz_add (tmp1, tmp1, tmp2); |

611 | if (mpz_cmp (tmp1, s) < 0) |

612 | { |

613 | mpz_set (s, tmp1); /* s = u^2 + v^2 */ |

614 | mpz_set (hp, u); |

615 | mpz_set (pp, v); |

616 | } |

617 | mpf_set_z (f_tmp1, s); |

618 | mpf_sqrt (rop[ui_t - 1], f_tmp1); |

619 | |

620 | /* S4 [Advance t.] */ |

621 | mpz_neg (U[0][0], h); |

622 | mpz_set (U[0][1], p); |

623 | mpz_neg (U[1][0], hp); |

624 | mpz_set (U[1][1], pp); |

625 | |

626 | mpz_set (V[0][0], pp); |

627 | mpz_set (V[0][1], hp); |

628 | mpz_neg (V[1][0], p); |

629 | mpz_neg (V[1][1], h); |

630 | if (mpz_cmp_ui (pp, 0) > 0) |

631 | { |

632 | mpz_neg (V[0][0], V[0][0]); |

633 | mpz_neg (V[0][1], V[0][1]); |

634 | mpz_neg (V[1][0], V[1][0]); |

635 | mpz_neg (V[1][1], V[1][1]); |

636 | } |

637 | |

638 | while (ui_t + 1 != T) /* S4 loop */ |

639 | { |

640 | ui_t++; |

641 | mpz_mul (r, a, r); |

642 | mpz_mod (r, r, m); |

643 | |

644 | /* Add new row and column to U and V. They are initialized with |

645 | all elements set to zero, so clearing is not necessary. */ |

646 | |

647 | mpz_neg (U[ui_t][0], r); /* U: First col in new row. */ |

648 | mpz_set_ui (U[ui_t][ui_t], 1); /* U: Last col in new row. */ |

649 | |

650 | mpz_set (V[ui_t][ui_t], m); /* V: Last col in new row. */ |

651 | |

652 | /* "Finally, for 1 <= i < t, |

653 | set q = round (vi1 * r / m), |

654 | vit = vi1*r - q*m, |

655 | and Ut=Ut+q*Ui */ |

656 | |

657 | for (ui_i = 0; ui_i < ui_t; ui_i++) |

658 | { |

659 | mpz_mul (tmp1, V[ui_i][0], r); /* tmp1=vi1*r */ |

660 | zdiv_round (q, tmp1, m); /* q=round(vi1*r/m) */ |

661 | mpz_mul (tmp2, q, m); /* tmp2=q*m */ |

662 | mpz_sub (V[ui_i][ui_t], tmp1, tmp2); |

663 | |

664 | for (ui_j = 0; ui_j <= ui_t; ui_j++) /* U[t] = U[t] + q*U[i] */ |

665 | { |

666 | mpz_mul (tmp1, q, U[ui_i][ui_j]); /* tmp=q*uij */ |

667 | mpz_add (U[ui_t][ui_j], U[ui_t][ui_j], tmp1); /* utj = utj + q*uij */ |

668 | } |

669 | } |

670 | |

671 | /* s = min (s, zdot (U[t], U[t]) */ |

672 | vz_dot (tmp1, U[ui_t], U[ui_t], ui_t + 1); |

673 | if (mpz_cmp (tmp1, s) < 0) |

674 | mpz_set (s, tmp1); |

675 | |

676 | ui_k = ui_t; |

677 | ui_j = 0; /* WARNING: ui_j no longer a temp. */ |

678 | |

679 | /* S5 [Transform.] */ |

680 | if (g_debug > DEBUG_2) |

681 | printf ("(t, k, j, q1, q2, ...)\n"); |

682 | do |

683 | { |

684 | if (g_debug > DEBUG_2) |

685 | printf ("(%u, %u, %u", ui_t + 1, ui_k + 1, ui_j + 1); |

686 | |

687 | for (ui_i = 0; ui_i <= ui_t; ui_i++) |

688 | { |

689 | if (ui_i != ui_j) |

690 | { |

691 | vz_dot (tmp1, V[ui_i], V[ui_j], ui_t + 1); /* tmp1=dot(Vi,Vj). */ |

692 | mpz_abs (tmp2, tmp1); |

693 | mpz_mul_ui (tmp2, tmp2, 2); /* tmp2 = 2*abs(dot(Vi,Vj) */ |

694 | vz_dot (tmp3, V[ui_j], V[ui_j], ui_t + 1); /* tmp3=dot(Vj,Vj). */ |

695 | |

696 | if (mpz_cmp (tmp2, tmp3) > 0) |

697 | { |

698 | zdiv_round (q, tmp1, tmp3); /* q=round(Vi.Vj/Vj.Vj) */ |

699 | if (g_debug > DEBUG_2) |

700 | { |

701 | printf (", "); |

702 | mpz_out_str (stdout, 10, q); |

703 | } |

704 | |

705 | for (ui_l = 0; ui_l <= ui_t; ui_l++) |

706 | { |

707 | mpz_mul (tmp1, q, V[ui_j][ui_l]); |

708 | mpz_sub (V[ui_i][ui_l], V[ui_i][ui_l], tmp1); /* Vi=Vi-q*Vj */ |

709 | mpz_mul (tmp1, q, U[ui_i][ui_l]); |

710 | mpz_add (U[ui_j][ui_l], U[ui_j][ui_l], tmp1); /* Uj=Uj+q*Ui */ |

711 | } |

712 | |

713 | vz_dot (tmp1, U[ui_j], U[ui_j], ui_t + 1); /* tmp1=dot(Uj,Uj) */ |

714 | if (mpz_cmp (tmp1, s) < 0) /* s = min(s,dot(Uj,Uj)) */ |

715 | mpz_set (s, tmp1); |

716 | ui_k = ui_j; |

717 | } |

718 | else if (g_debug > DEBUG_2) |

719 | printf (", #"); /* 2|Vi.Vj| <= Vj.Vj */ |

720 | } |

721 | else if (g_debug > DEBUG_2) |

722 | printf (", *"); /* i == j */ |

723 | } |

724 | |

725 | if (g_debug > DEBUG_2) |

726 | printf (")\n"); |

727 | |

728 | /* S6 [Advance j.] */ |

729 | if (ui_j == ui_t) |

730 | ui_j = 0; |

731 | else |

732 | ui_j++; |

733 | } |

734 | while (ui_j != ui_k); /* S5 */ |

735 | |

736 | /* From Knuth p. 104: "The exhaustive search in steps S8-S10 |

737 | reduces the value of s only rarely." */ |

738 | #ifdef DO_SEARCH |

739 | /* S7 [Prepare for search.] */ |

740 | /* Find minimum in (x[1], ..., x[t]) satisfying condition |

741 | x[k]^2 <= f(y[1], ...,y[t]) * dot(V[k],V[k]) */ |

742 | |

743 | ui_k = ui_t; |

744 | if (g_debug > DEBUG_2) |

745 | { |

746 | printf ("searching..."); |

747 | /*for (f = 0; f < ui_t*/ |

748 | fflush (stdout); |

749 | } |

750 | |

751 | /* Z[i] = floor (sqrt (floor (dot(V[i],V[i]) * s / m^2))); */ |

752 | mpz_pow_ui (tmp1, m, 2); |

753 | mpf_set_z (f_tmp1, tmp1); |

754 | mpf_set_z (f_tmp2, s); |

755 | mpf_div (f_tmp1, f_tmp2, f_tmp1); /* f_tmp1 = s/m^2 */ |

756 | for (ui_i = 0; ui_i <= ui_t; ui_i++) |

757 | { |

758 | vz_dot (tmp1, V[ui_i], V[ui_i], ui_t + 1); |

759 | mpf_set_z (f_tmp2, tmp1); |

760 | mpf_mul (f_tmp2, f_tmp2, f_tmp1); |

761 | f_floor (f_tmp2, f_tmp2); |

762 | mpf_sqrt (f_tmp2, f_tmp2); |

763 | mpz_set_f (Z[ui_i], f_tmp2); |

764 | } |

765 | |

766 | /* S8 [Advance X[k].] */ |

767 | do |

768 | { |

769 | if (g_debug > DEBUG_2) |

770 | { |

771 | printf ("X[%u] = ", ui_k); |

772 | mpz_out_str (stdout, 10, X[ui_k]); |

773 | printf ("\tZ[%u] = ", ui_k); |

774 | mpz_out_str (stdout, 10, Z[ui_k]); |

775 | printf ("\n"); |

776 | fflush (stdout); |

777 | } |

778 | |

779 | if (mpz_cmp (X[ui_k], Z[ui_k])) |

780 | { |

781 | mpz_add_ui (X[ui_k], X[ui_k], 1); |

782 | for (ui_i = 0; ui_i <= ui_t; ui_i++) |

783 | mpz_add (Y[ui_i], Y[ui_i], U[ui_k][ui_i]); |

784 | |

785 | /* S9 [Advance k.] */ |

786 | while (++ui_k <= ui_t) |

787 | { |

788 | mpz_neg (X[ui_k], Z[ui_k]); |

789 | mpz_mul_ui (tmp1, Z[ui_k], 2); |

790 | for (ui_i = 0; ui_i <= ui_t; ui_i++) |

791 | { |

792 | mpz_mul (tmp2, tmp1, U[ui_k][ui_i]); |

793 | mpz_sub (Y[ui_i], Y[ui_i], tmp2); |

794 | } |

795 | } |

796 | vz_dot (tmp1, Y, Y, ui_t + 1); |

797 | if (mpz_cmp (tmp1, s) < 0) |

798 | mpz_set (s, tmp1); |

799 | } |

800 | } |

801 | while (--ui_k); |

802 | #endif /* DO_SEARCH */ |

803 | mpf_set_z (f_tmp1, s); |

804 | mpf_sqrt (rop[ui_t - 1], f_tmp1); |

805 | #ifdef DO_SEARCH |

806 | if (g_debug > DEBUG_2) |

807 | printf ("done.\n"); |

808 | #endif /* DO_SEARCH */ |

809 | } /* S4 loop */ |

810 | |

811 | /* Clear GMP variables. */ |

812 | |

813 | mpf_clear (f_tmp1); |

814 | mpf_clear (f_tmp2); |

815 | for (ui_i = 0; ui_i < GMP_SPECT_MAXT; ui_i++) |

816 | { |

817 | for (ui_j = 0; ui_j < GMP_SPECT_MAXT; ui_j++) |

818 | { |

819 | mpz_clear (U[ui_i][ui_j]); |

820 | mpz_clear (V[ui_i][ui_j]); |

821 | } |

822 | mpz_clear (X[ui_i]); |

823 | mpz_clear (Y[ui_i]); |

824 | mpz_clear (Z[ui_i]); |

825 | } |

826 | mpz_clear (tmp1); |

827 | mpz_clear (tmp2); |

828 | mpz_clear (tmp3); |

829 | mpz_clear (h); |

830 | mpz_clear (hp); |

831 | mpz_clear (r); |

832 | mpz_clear (s); |

833 | mpz_clear (p); |

834 | mpz_clear (pp); |

835 | mpz_clear (q); |

836 | mpz_clear (u); |

837 | mpz_clear (v); |

838 | |

839 | return; |

840 | } |

841 |

**Note:**See TracBrowser for help on using the repository browser.