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