tests/rand/statlib.c

   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 3 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.  If not, see http://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
  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
 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
 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
 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
 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
 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
 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
 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
 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
 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
 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
 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
 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
 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 }