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 test suite.
8
9 The GNU MP Library test suite is free software; you can redistribute it
10 and/or modify it under the terms of the GNU General Public License as
11 published by the Free Software Foundation; either version 3 of the License,
12 or (at your option) any later version.
13
14 The GNU MP Library test suite is distributed in the hope that it will be
15 useful, but WITHOUT ANY WARRANTY; without even the implied warranty of
16 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
17 Public License for more details.
18
19 You should have received a copy of the GNU General Public License along with
20 the GNU MP Library test suite. If not, see https://www.gnu.org/licenses/. */
21
22 /* The theories for these functions are taken from D. Knuth's "The Art
23 of Computer Programming: Volume 2, Seminumerical Algorithms", Third
24 Edition, Addison Wesley, 1998. */
25
26 /* Implementation notes.
27
28 The Kolmogorov-Smirnov test.
29
30 Eq. (13) in Knuth, p. 50, says that if X1, X2, ..., Xn are independent
31 observations arranged into ascending order
32
33 Kp = sqr(n) * max(j/n - F(Xj)) for all 1<=j<=n
34 Km = sqr(n) * max(F(Xj) - (j-1)/n)) for all 1<=j<=n
35
36 where F(x) = Pr(X <= x) = probability that (X <= x), which for a
37 uniformly distributed random real number between zero and one is
38 exactly the number itself (x).
39
40
41 The answer to exercise 23 gives the following implementation, which
42 doesn't need the observations to be sorted in ascending order:
43
44 for (k = 0; k < m; k++)
45 a[k] = 1.0
46 b[k] = 0.0
47 c[k] = 0
48
49 for (each observation Xj)
50 Y = F(Xj)
51 k = floor (m * Y)
52 a[k] = min (a[k], Y)
53 b[k] = max (b[k], Y)
54 c[k] += 1
55
56 j = 0
57 rp = rm = 0
58 for (k = 0; k < m; k++)
59 if (c[k] > 0)
60 rm = max (rm, a[k] - j/n)
61 j += c[k]
62 rp = max (rp, j/n - b[k])
63
64 Kp = sqr (n) * rp
65 Km = sqr (n) * rm
66
67 */
68
69 #include <stdio.h>
70 #include <stdlib.h>
71 #include <math.h>
72
73 #include "gmpstat.h"
74
75 /* ks (Kp, Km, X, P, n) -- Perform a Kolmogorov-Smirnov test on the N
76 real numbers between zero and one in vector X. P is the
77 distribution function, called for each entry in X, which should
78 calculate the probability of X being greater than or equal to any
79 number in the sequence. (For a uniformly distributed sequence of
80 real numbers between zero and one, this is simply equal to X.) The
81 result is put in Kp and Km. */
82
83 void
84 ks (mpf_t Kp,
85 mpf_t Km,
86 mpf_t X[],
87 void (P) (mpf_t, mpf_t),
88 unsigned long int n)
89 {
90 mpf_t Kt; /* temp */
91 mpf_t f_x;
92 mpf_t f_j; /* j */
93 mpf_t f_jnq; /* j/n or (j-1)/n */
94 unsigned long int j;
95
96 /* Sort the vector in ascending order. */
97 qsort (X, n, sizeof (__mpf_struct), mpf_cmp);
98
99 /* K-S test. */
100 /* Kp = sqr(n) * max(j/n - F(Xj)) for all 1<=j<=n
101 Km = sqr(n) * max(F(Xj) - (j-1)/n)) for all 1<=j<=n
102 */
103
104 mpf_init (Kt); mpf_init (f_x); mpf_init (f_j); mpf_init (f_jnq);
105 mpf_set_ui (Kp, 0); mpf_set_ui (Km, 0);
106 for (j = 1; j <= n; j++)
107 {
108 P (f_x, X[j-1]);
109 mpf_set_ui (f_j, j);
110
111 mpf_div_ui (f_jnq, f_j, n);
112 mpf_sub (Kt, f_jnq, f_x);
113 if (mpf_cmp (Kt, Kp) > 0)
114 mpf_set (Kp, Kt);
115 if (g_debug > DEBUG_2)
116 {
117 printf ("j=%lu ", j);
118 printf ("P()="); mpf_out_str (stdout, 10, 2, f_x); printf ("\t");
119
120 printf ("jnq="); mpf_out_str (stdout, 10, 2, f_jnq); printf (" ");
121 printf ("diff="); mpf_out_str (stdout, 10, 2, Kt); printf (" ");
122 printf ("Kp="); mpf_out_str (stdout, 10, 2, Kp); printf ("\t");
123 }
124 mpf_sub_ui (f_j, f_j, 1);
125 mpf_div_ui (f_jnq, f_j, n);
126 mpf_sub (Kt, f_x, f_jnq);
127 if (mpf_cmp (Kt, Km) > 0)
128 mpf_set (Km, Kt);
129
130 if (g_debug > DEBUG_2)
131 {
132 printf ("jnq="); mpf_out_str (stdout, 10, 2, f_jnq); printf (" ");
133 printf ("diff="); mpf_out_str (stdout, 10, 2, Kt); printf (" ");
134 printf ("Km="); mpf_out_str (stdout, 10, 2, Km); printf (" ");
135 printf ("\n");
136 }
137 }
138 mpf_sqrt_ui (Kt, n);
139 mpf_mul (Kp, Kp, Kt);
140 mpf_mul (Km, Km, Kt);
141
142 mpf_clear (Kt); mpf_clear (f_x); mpf_clear (f_j); mpf_clear (f_jnq);
143 }
144
145 /* ks_table(val, n) -- calculate probability for Kp/Km less than or
146 equal to VAL with N observations. See [Knuth section 3.3.1] */
147
148 void
149 ks_table (mpf_t p, mpf_t val, const unsigned int n)
150 {
151 /* We use Eq. (27), Knuth p.58, skipping O(1/n) for simplicity.
152 This shortcut will result in too high probabilities, especially
153 when n is small.
154
155 Pr(Kp(n) <= s) = 1 - pow(e, -2*s^2) * (1 - 2/3*s/sqrt(n) + O(1/n)) */
156
157 /* We have 's' in variable VAL and store the result in P. */
158
159 mpf_t t1, t2;
160
161 mpf_init (t1); mpf_init (t2);
162
163 /* t1 = 1 - 2/3 * s/sqrt(n) */
164 mpf_sqrt_ui (t1, n);
165 mpf_div (t1, val, t1);
166 mpf_mul_ui (t1, t1, 2);
167 mpf_div_ui (t1, t1, 3);
168 mpf_ui_sub (t1, 1, t1);
169
170 /* t2 = pow(e, -2*s^2) */
171 #ifndef OLDGMP
172 mpf_pow_ui (t2, val, 2); /* t2 = s^2 */
173 mpf_set_d (t2, exp (-(2.0 * mpf_get_d (t2))));
174 #else
175 /* hmmm, gmp doesn't have pow() for floats. use doubles. */
176 mpf_set_d (t2, pow (M_E, -(2 * pow (mpf_get_d (val), 2))));
177 #endif
178
179 /* p = 1 - t1 * t2 */
180 mpf_mul (t1, t1, t2);
181 mpf_ui_sub (p, 1, t1);
182
183 mpf_clear (t1); mpf_clear (t2);
184 }
185
186 static double x2_table_X[][7] = {
187 { -2.33, -1.64, -.674, 0.0, 0.674, 1.64, 2.33 }, /* x */
188 { 5.4289, 2.6896, .454276, 0.0, .454276, 2.6896, 5.4289} /* x^2 */
189 };
190
191 #define _2D3 ((double) .6666666666)
192
193 /* x2_table (t, v, n) -- return chi-square table row for V in T[]. */
194 void
195 x2_table (double t[],
196 unsigned int v)
197 {
198 int f;
199
200
201 /* FIXME: Do a table lookup for v <= 30 since the following formula
202 [Knuth, vol 2, 3.3.1] is only good for v > 30. */
203
204 /* value = v + sqrt(2*v) * X[p] + (2/3) * X[p]^2 - 2/3 + O(1/sqrt(t) */
205 /* NOTE: The O() term is ignored for simplicity. */
206
207 for (f = 0; f < 7; f++)
208 t[f] =
209 v +
210 sqrt (2 * v) * x2_table_X[0][f] +
211 _2D3 * x2_table_X[1][f] - _2D3;
212 }
213
214
215 /* P(p, x) -- Distribution function. Calculate the probability of X
216 being greater than or equal to any number in the sequence. For a
217 random real number between zero and one given by a uniformly
218 distributed random number generator, this is simply equal to X. */
219
220 static void
221 P (mpf_t p, mpf_t x)
222 {
223 mpf_set (p, x);
224 }
225
226 /* mpf_freqt() -- Frequency test using KS on N real numbers between zero
227 and one. See [Knuth vol 2, p.61]. */
228 void
229 mpf_freqt (mpf_t Kp,
230 mpf_t Km,
231 mpf_t X[],
232 const unsigned long int n)
233 {
234 ks (Kp, Km, X, P, n);
235 }
236
237
238 /* The Chi-square test. Eq. (8) in Knuth vol. 2 says that if Y[]
239 holds the observations and p[] is the probability for.. (to be
240 continued!)
241
242 V = 1/n * sum((s=1 to k) Y[s]^2 / p[s]) - n */
243
244 void
245 x2 (mpf_t V, /* result */
246 unsigned long int X[], /* data */
247 unsigned int k, /* #of categories */
248 void (P) (mpf_t, unsigned long int, void *), /* probability func */
249 void *x, /* extra user data passed to P() */
250 unsigned long int n) /* #of samples */
251 {
252 unsigned int f;
253 mpf_t f_t, f_t2; /* temp floats */
254
255 mpf_init (f_t); mpf_init (f_t2);
256
257
258 mpf_set_ui (V, 0);
259 for (f = 0; f < k; f++)
260 {
261 if (g_debug > DEBUG_2)
262 fprintf (stderr, "%u: P()=", f);
263 mpf_set_ui (f_t, X[f]);
264 mpf_mul (f_t, f_t, f_t); /* f_t = X[f]^2 */
265 P (f_t2, f, x); /* f_t2 = Pr(f) */
266 if (g_debug > DEBUG_2)
267 mpf_out_str (stderr, 10, 2, f_t2);
268 mpf_div (f_t, f_t, f_t2);
269 mpf_add (V, V, f_t);
270 if (g_debug > DEBUG_2)
271 {
272 fprintf (stderr, "\tV=");
273 mpf_out_str (stderr, 10, 2, V);
274 fprintf (stderr, "\t");
275 }
276 }
277 if (g_debug > DEBUG_2)
278 fprintf (stderr, "\n");
279 mpf_div_ui (V, V, n);
280 mpf_sub_ui (V, V, n);
281
282 mpf_clear (f_t); mpf_clear (f_t2);
283 }
284
285 /* Pzf(p, s, x) -- Probability for category S in mpz_freqt(). It's
286 1/d for all S. X is a pointer to an unsigned int holding 'd'. */
287 static void
288 Pzf (mpf_t p, unsigned long int s, void *x)
289 {
290 mpf_set_ui (p, 1);
291 mpf_div_ui (p, p, *((unsigned int *) x));
292 }
293
294 /* mpz_freqt(V, X, imax, n) -- Frequency test on integers. [Knuth,
295 vol 2, 3.3.2]. Keep IMAX low on this one, since we loop from 0 to
296 IMAX. 128 or 256 could be nice.
297
298 X[] must not contain numbers outside the range 0 <= X <= IMAX.
299
300 Return value is number of observations actually used, after
301 discarding entries out of range.
302
303 Since X[] contains integers between zero and IMAX, inclusive, we
304 have IMAX+1 categories.
305
306 Note that N should be at least 5*IMAX. Result is put in V and can
307 be compared to output from x2_table (v=IMAX). */
308
309 unsigned long int
310 mpz_freqt (mpf_t V,
311 mpz_t X[],
312 unsigned int imax,
313 const unsigned long int n)
314 {
315 unsigned long int *v; /* result */
316 unsigned int f;
317 unsigned int d; /* number of categories = imax+1 */
318 unsigned int uitemp;
319 unsigned long int usedn;
320
321
322 d = imax + 1;
323
324 v = (unsigned long int *) calloc (imax + 1, sizeof (unsigned long int));
325 if (NULL == v)
326 {
327 fprintf (stderr, "mpz_freqt(): out of memory\n");
328 exit (1);
329 }
330
331 /* count */
332 usedn = n; /* actual number of observations */
333 for (f = 0; f < n; f++)
334 {
335 uitemp = mpz_get_ui(X[f]);
336 if (uitemp > imax) /* sanity check */
337 {
338 if (g_debug)
339 fprintf (stderr, "mpz_freqt(): warning: input insanity: %u, "\
340 "ignored.\n", uitemp);
341 usedn--;
342 continue;
343 }
344 v[uitemp]++;
345 }
346
347 if (g_debug > DEBUG_2)
348 {
349 fprintf (stderr, "counts:\n");
350 for (f = 0; f <= imax; f++)
351 fprintf (stderr, "%u:\t%lu\n", f, v[f]);
352 }
353
354 /* chi-square with k=imax+1 and P(x)=1/(imax+1) for all x.*/
355 x2 (V, v, d, Pzf, (void *) &d, usedn);
356
357 free (v);
358 return (usedn);
359 }
360
361 /* debug dummy to drag in dump funcs */
362 void
363 foo_debug ()
364 {
365 if (0)
366 {
367 mpf_dump (0);
368 #ifndef OLDGMP
369 mpz_dump (0);
370 #endif
371 }
372 }
373
374 /* merit (rop, t, v, m) -- calculate merit for spectral test result in
375 dimension T, see Knuth p. 105. BUGS: Only valid for 2 <= T <=
376 6. */
377 void
378 merit (mpf_t rop, unsigned int t, mpf_t v, mpz_t m)
379 {
380 int f;
381 mpf_t f_m, f_const, f_pi;
382
383 mpf_init (f_m);
384 mpf_set_z (f_m, m);
385 mpf_init_set_d (f_const, M_PI);
386 mpf_init_set_d (f_pi, M_PI);
387
388 switch (t)
389 {
390 case 2: /* PI */
391 break;
392 case 3: /* PI * 4/3 */
393 mpf_mul_ui (f_const, f_const, 4);
394 mpf_div_ui (f_const, f_const, 3);
395 break;
396 case 4: /* PI^2 * 1/2 */
397 mpf_mul (f_const, f_const, f_pi);
398 mpf_div_ui (f_const, f_const, 2);
399 break;
400 case 5: /* PI^2 * 8/15 */
401 mpf_mul (f_const, f_const, f_pi);
402 mpf_mul_ui (f_const, f_const, 8);
403 mpf_div_ui (f_const, f_const, 15);
404 break;
405 case 6: /* PI^3 * 1/6 */
406 mpf_mul (f_const, f_const, f_pi);
407 mpf_mul (f_const, f_const, f_pi);
408 mpf_div_ui (f_const, f_const, 6);
409 break;
410 default:
411 fprintf (stderr,
412 "spect (merit): can't calculate merit for dimensions > 6\n");
413 mpf_set_ui (f_const, 0);
414 break;
415 }
416
417 /* rop = v^t */
418 mpf_set (rop, v);
419 for (f = 1; f < t; f++)
420 mpf_mul (rop, rop, v);
421 mpf_mul (rop, rop, f_const);
422 mpf_div (rop, rop, f_m);
423
424 mpf_clear (f_m);
425 mpf_clear (f_const);
426 mpf_clear (f_pi);
427 }
428
429 double
430 merit_u (unsigned int t, mpf_t v, mpz_t m)
431 {
432 mpf_t rop;
433 double res;
434
435 mpf_init (rop);
436 merit (rop, t, v, m);
437 res = mpf_get_d (rop);
438 mpf_clear (rop);
439 return res;
440 }
441
442 /* f_floor (rop, op) -- Set rop = floor (op). */
443 void
444 f_floor (mpf_t rop, mpf_t op)
445 {
446 mpz_t z;
447
448 mpz_init (z);
449
450 /* No mpf_floor(). Convert to mpz and back. */
451 mpz_set_f (z, op);
452 mpf_set_z (rop, z);
453
454 mpz_clear (z);
455 }
456
457
458 /* vz_dot (rop, v1, v2, nelem) -- compute dot product of z-vectors V1,
459 V2. N is number of elements in vectors V1 and V2. */
460
461 void
462 vz_dot (mpz_t rop, mpz_t V1[], mpz_t V2[], unsigned int n)
463 {
464 mpz_t t;
465
466 mpz_init (t);
467 mpz_set_ui (rop, 0);
468 while (n--)
469 {
470 mpz_mul (t, V1[n], V2[n]);
471 mpz_add (rop, rop, t);
472 }
473
474 mpz_clear (t);
475 }
476
477 void
478 spectral_test (mpf_t rop[], unsigned int T, mpz_t a, mpz_t m)
479 {
480 /* Knuth "Seminumerical Algorithms, Third Edition", section 3.3.4
481 (pp. 101-103). */
482
483 /* v[t] = min { sqrt (x[1]^2 + ... + x[t]^2) |
484 x[1] + a*x[2] + ... + pow (a, t-1) * x[t] is congruent to 0 (mod m) } */
485
486
487 /* Variables. */
488 unsigned int ui_t;
489 unsigned int ui_i, ui_j, ui_k, ui_l;
490 mpf_t f_tmp1, f_tmp2;
491 mpz_t tmp1, tmp2, tmp3;
492 mpz_t U[GMP_SPECT_MAXT][GMP_SPECT_MAXT],
493 V[GMP_SPECT_MAXT][GMP_SPECT_MAXT],
494 X[GMP_SPECT_MAXT],
495 Y[GMP_SPECT_MAXT],
496 Z[GMP_SPECT_MAXT];
497 mpz_t h, hp, r, s, p, pp, q, u, v;
498
499 /* GMP inits. */
500 mpf_init (f_tmp1);
501 mpf_init (f_tmp2);
502 for (ui_i = 0; ui_i < GMP_SPECT_MAXT; ui_i++)
503 {
504 for (ui_j = 0; ui_j < GMP_SPECT_MAXT; ui_j++)
505 {
506 mpz_init_set_ui (U[ui_i][ui_j], 0);
507 mpz_init_set_ui (V[ui_i][ui_j], 0);
508 }
509 mpz_init_set_ui (X[ui_i], 0);
510 mpz_init_set_ui (Y[ui_i], 0);
511 mpz_init (Z[ui_i]);
512 }
513 mpz_init (tmp1);
514 mpz_init (tmp2);
515 mpz_init (tmp3);
516 mpz_init (h);
517 mpz_init (hp);
518 mpz_init (r);
519 mpz_init (s);
520 mpz_init (p);
521 mpz_init (pp);
522 mpz_init (q);
523 mpz_init (u);
524 mpz_init (v);
525
526 /* Implementation inits. */
527 if (T > GMP_SPECT_MAXT)
528 T = GMP_SPECT_MAXT; /* FIXME: Lazy. */
529
530 /* S1 [Initialize.] */
531 ui_t = 2 - 1; /* NOTE: `t' in description == ui_t + 1
532 for easy indexing */
533 mpz_set (h, a);
534 mpz_set (hp, m);
535 mpz_set_ui (p, 1);
536 mpz_set_ui (pp, 0);
537 mpz_set (r, a);
538 mpz_pow_ui (s, a, 2);
539 mpz_add_ui (s, s, 1); /* s = 1 + a^2 */
540
541 /* S2 [Euclidean step.] */
542 while (1)
543 {
544 if (g_debug > DEBUG_1)
545 {
546 mpz_mul (tmp1, h, pp);
547 mpz_mul (tmp2, hp, p);
548 mpz_sub (tmp1, tmp1, tmp2);
549 if (mpz_cmpabs (m, tmp1))
550 {
551 printf ("***BUG***: h*pp - hp*p = ");
552 mpz_out_str (stdout, 10, tmp1);
553 printf ("\n");
554 }
555 }
556 if (g_debug > DEBUG_2)
557 {
558 printf ("hp = ");
559 mpz_out_str (stdout, 10, hp);
560 printf ("\nh = ");
561 mpz_out_str (stdout, 10, h);
562 printf ("\n");
563 fflush (stdout);
564 }
565
566 if (mpz_sgn (h))
567 mpz_tdiv_q (q, hp, h); /* q = floor(hp/h) */
568 else
569 mpz_set_ui (q, 1);
570
571 if (g_debug > DEBUG_2)
572 {
573 printf ("q = ");
574 mpz_out_str (stdout, 10, q);
575 printf ("\n");
576 fflush (stdout);
577 }
578
579 mpz_mul (tmp1, q, h);
580 mpz_sub (u, hp, tmp1); /* u = hp - q*h */
581
582 mpz_mul (tmp1, q, p);
583 mpz_sub (v, pp, tmp1); /* v = pp - q*p */
584
585 mpz_pow_ui (tmp1, u, 2);
586 mpz_pow_ui (tmp2, v, 2);
587 mpz_add (tmp1, tmp1, tmp2);
588 if (mpz_cmp (tmp1, s) < 0)
589 {
590 mpz_set (s, tmp1); /* s = u^2 + v^2 */
591 mpz_set (hp, h); /* hp = h */
592 mpz_set (h, u); /* h = u */
593 mpz_set (pp, p); /* pp = p */
594 mpz_set (p, v); /* p = v */
595 }
596 else
597 break;
598 }
599
600 /* S3 [Compute v2.] */
601 mpz_sub (u, u, h);
602 mpz_sub (v, v, p);
603
604 mpz_pow_ui (tmp1, u, 2);
605 mpz_pow_ui (tmp2, v, 2);
606 mpz_add (tmp1, tmp1, tmp2);
607 if (mpz_cmp (tmp1, s) < 0)
608 {
609 mpz_set (s, tmp1); /* s = u^2 + v^2 */
610 mpz_set (hp, u);
611 mpz_set (pp, v);
612 }
613 mpf_set_z (f_tmp1, s);
614 mpf_sqrt (rop[ui_t - 1], f_tmp1);
615
616 /* S4 [Advance t.] */
617 mpz_neg (U[0][0], h);
618 mpz_set (U[0][1], p);
619 mpz_neg (U[1][0], hp);
620 mpz_set (U[1][1], pp);
621
622 mpz_set (V[0][0], pp);
623 mpz_set (V[0][1], hp);
624 mpz_neg (V[1][0], p);
625 mpz_neg (V[1][1], h);
626 if (mpz_cmp_ui (pp, 0) > 0)
627 {
628 mpz_neg (V[0][0], V[0][0]);
629 mpz_neg (V[0][1], V[0][1]);
630 mpz_neg (V[1][0], V[1][0]);
631 mpz_neg (V[1][1], V[1][1]);
632 }
633
634 while (ui_t + 1 != T) /* S4 loop */
635 {
636 ui_t++;
637 mpz_mul (r, a, r);
638 mpz_mod (r, r, m);
639
640 /* Add new row and column to U and V. They are initialized with
641 all elements set to zero, so clearing is not necessary. */
642
643 mpz_neg (U[ui_t][0], r); /* U: First col in new row. */
644 mpz_set_ui (U[ui_t][ui_t], 1); /* U: Last col in new row. */
645
646 mpz_set (V[ui_t][ui_t], m); /* V: Last col in new row. */
647
648 /* "Finally, for 1 <= i < t,
649 set q = round (vi1 * r / m),
650 vit = vi1*r - q*m,
651 and Ut=Ut+q*Ui */
652
653 for (ui_i = 0; ui_i < ui_t; ui_i++)
654 {
655 mpz_mul (tmp1, V[ui_i][0], r); /* tmp1=vi1*r */
656 zdiv_round (q, tmp1, m); /* q=round(vi1*r/m) */
657 mpz_mul (tmp2, q, m); /* tmp2=q*m */
658 mpz_sub (V[ui_i][ui_t], tmp1, tmp2);
659
660 for (ui_j = 0; ui_j <= ui_t; ui_j++) /* U[t] = U[t] + q*U[i] */
661 {
662 mpz_mul (tmp1, q, U[ui_i][ui_j]); /* tmp=q*uij */
663 mpz_add (U[ui_t][ui_j], U[ui_t][ui_j], tmp1); /* utj = utj + q*uij */
664 }
665 }
666
667 /* s = min (s, zdot (U[t], U[t]) */
668 vz_dot (tmp1, U[ui_t], U[ui_t], ui_t + 1);
669 if (mpz_cmp (tmp1, s) < 0)
670 mpz_set (s, tmp1);
671
672 ui_k = ui_t;
673 ui_j = 0; /* WARNING: ui_j no longer a temp. */
674
675 /* S5 [Transform.] */
676 if (g_debug > DEBUG_2)
677 printf ("(t, k, j, q1, q2, ...)\n");
678 do
679 {
680 if (g_debug > DEBUG_2)
681 printf ("(%u, %u, %u", ui_t + 1, ui_k + 1, ui_j + 1);
682
683 for (ui_i = 0; ui_i <= ui_t; ui_i++)
684 {
685 if (ui_i != ui_j)
686 {
687 vz_dot (tmp1, V[ui_i], V[ui_j], ui_t + 1); /* tmp1=dot(Vi,Vj). */
688 mpz_abs (tmp2, tmp1);
689 mpz_mul_ui (tmp2, tmp2, 2); /* tmp2 = 2*abs(dot(Vi,Vj) */
690 vz_dot (tmp3, V[ui_j], V[ui_j], ui_t + 1); /* tmp3=dot(Vj,Vj). */
691
692 if (mpz_cmp (tmp2, tmp3) > 0)
693 {
694 zdiv_round (q, tmp1, tmp3); /* q=round(Vi.Vj/Vj.Vj) */
695 if (g_debug > DEBUG_2)
696 {
697 printf (", ");
698 mpz_out_str (stdout, 10, q);
699 }
700
701 for (ui_l = 0; ui_l <= ui_t; ui_l++)
702 {
703 mpz_mul (tmp1, q, V[ui_j][ui_l]);
704 mpz_sub (V[ui_i][ui_l], V[ui_i][ui_l], tmp1); /* Vi=Vi-q*Vj */
705 mpz_mul (tmp1, q, U[ui_i][ui_l]);
706 mpz_add (U[ui_j][ui_l], U[ui_j][ui_l], tmp1); /* Uj=Uj+q*Ui */
707 }
708
709 vz_dot (tmp1, U[ui_j], U[ui_j], ui_t + 1); /* tmp1=dot(Uj,Uj) */
710 if (mpz_cmp (tmp1, s) < 0) /* s = min(s,dot(Uj,Uj)) */
711 mpz_set (s, tmp1);
712 ui_k = ui_j;
713 }
714 else if (g_debug > DEBUG_2)
715 printf (", #"); /* 2|Vi.Vj| <= Vj.Vj */
716 }
717 else if (g_debug > DEBUG_2)
718 printf (", *"); /* i == j */
719 }
720
721 if (g_debug > DEBUG_2)
722 printf (")\n");
723
724 /* S6 [Advance j.] */
725 if (ui_j == ui_t)
726 ui_j = 0;
727 else
728 ui_j++;
729 }
730 while (ui_j != ui_k); /* S5 */
731
732 /* From Knuth p. 104: "The exhaustive search in steps S8-S10
733 reduces the value of s only rarely." */
734 #ifdef DO_SEARCH
735 /* S7 [Prepare for search.] */
736 /* Find minimum in (x[1], ..., x[t]) satisfying condition
737 x[k]^2 <= f(y[1], ...,y[t]) * dot(V[k],V[k]) */
738
739 ui_k = ui_t;
740 if (g_debug > DEBUG_2)
741 {
742 printf ("searching...");
743 /*for (f = 0; f < ui_t*/
744 fflush (stdout);
745 }
746
747 /* Z[i] = floor (sqrt (floor (dot(V[i],V[i]) * s / m^2))); */
748 mpz_pow_ui (tmp1, m, 2);
749 mpf_set_z (f_tmp1, tmp1);
750 mpf_set_z (f_tmp2, s);
751 mpf_div (f_tmp1, f_tmp2, f_tmp1); /* f_tmp1 = s/m^2 */
752 for (ui_i = 0; ui_i <= ui_t; ui_i++)
753 {
754 vz_dot (tmp1, V[ui_i], V[ui_i], ui_t + 1);
755 mpf_set_z (f_tmp2, tmp1);
756 mpf_mul (f_tmp2, f_tmp2, f_tmp1);
757 f_floor (f_tmp2, f_tmp2);
758 mpf_sqrt (f_tmp2, f_tmp2);
759 mpz_set_f (Z[ui_i], f_tmp2);
760 }
761
762 /* S8 [Advance X[k].] */
763 do
764 {
765 if (g_debug > DEBUG_2)
766 {
767 printf ("X[%u] = ", ui_k);
768 mpz_out_str (stdout, 10, X[ui_k]);
769 printf ("\tZ[%u] = ", ui_k);
770 mpz_out_str (stdout, 10, Z[ui_k]);
771 printf ("\n");
772 fflush (stdout);
773 }
774
775 if (mpz_cmp (X[ui_k], Z[ui_k]))
776 {
777 mpz_add_ui (X[ui_k], X[ui_k], 1);
778 for (ui_i = 0; ui_i <= ui_t; ui_i++)
779 mpz_add (Y[ui_i], Y[ui_i], U[ui_k][ui_i]);
780
781 /* S9 [Advance k.] */
782 while (++ui_k <= ui_t)
783 {
784 mpz_neg (X[ui_k], Z[ui_k]);
785 mpz_mul_ui (tmp1, Z[ui_k], 2);
786 for (ui_i = 0; ui_i <= ui_t; ui_i++)
787 {
788 mpz_mul (tmp2, tmp1, U[ui_k][ui_i]);
789 mpz_sub (Y[ui_i], Y[ui_i], tmp2);
790 }
791 }
792 vz_dot (tmp1, Y, Y, ui_t + 1);
793 if (mpz_cmp (tmp1, s) < 0)
794 mpz_set (s, tmp1);
795 }
796 }
797 while (--ui_k);
798 #endif /* DO_SEARCH */
799 mpf_set_z (f_tmp1, s);
800 mpf_sqrt (rop[ui_t - 1], f_tmp1);
801 #ifdef DO_SEARCH
802 if (g_debug > DEBUG_2)
803 printf ("done.\n");
804 #endif /* DO_SEARCH */
805 } /* S4 loop */
806
807 /* Clear GMP variables. */
808
809 mpf_clear (f_tmp1);
810 mpf_clear (f_tmp2);
811 for (ui_i = 0; ui_i < GMP_SPECT_MAXT; ui_i++)
812 {
813 for (ui_j = 0; ui_j < GMP_SPECT_MAXT; ui_j++)
814 {
815 mpz_clear (U[ui_i][ui_j]);
816 mpz_clear (V[ui_i][ui_j]);
817 }
818 mpz_clear (X[ui_i]);
819 mpz_clear (Y[ui_i]);
820 mpz_clear (Z[ui_i]);
821 }
822 mpz_clear (tmp1);
823 mpz_clear (tmp2);
824 mpz_clear (tmp3);
825 mpz_clear (h);
826 mpz_clear (hp);
827 mpz_clear (r);
828 mpz_clear (s);
829 mpz_clear (p);
830 mpz_clear (pp);
831 mpz_clear (q);
832 mpz_clear (u);
833 mpz_clear (v);
834
835 return;
836 }