1/* statlib.c -- Statistical functions for testing the randomness of 2 number sequences. */ 3 4/* 5Copyright 1999, 2000 Free Software Foundation, Inc. 6 7This file is part of the GNU MP Library test suite. 8 9The GNU MP Library test suite is free software; you can redistribute it 10and/or modify it under the terms of the GNU General Public License as 11published by the Free Software Foundation; either version 3 of the License, 12or (at your option) any later version. 13 14The GNU MP Library test suite is distributed in the hope that it will be 15useful, but WITHOUT ANY WARRANTY; without even the implied warranty of 16MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General 17Public License for more details. 18 19You should have received a copy of the GNU General Public License along with 20the 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 23of Computer Programming: Volume 2, Seminumerical Algorithms", Third 24Edition, Addison Wesley, 1998. */ 25 26/* Implementation notes. 27 28The Kolmogorov-Smirnov test. 29 30Eq. (13) in Knuth, p. 50, says that if X1, X2, ..., Xn are independent 31observations 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 36where F(x) = Pr(X <= x) = probability that (X <= x), which for a 37uniformly distributed random real number between zero and one is 38exactly the number itself (x). 39 40 41The answer to exercise 23 gives the following implementation, which 42doesn't need the observations to be sorted in ascending order: 43 44for (k = 0; k < m; k++) 45 a[k] = 1.0 46 b[k] = 0.0 47 c[k] = 0 48 49for (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 64Kp = sqr (n) * rp 65Km = 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 83void 84ks (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 148void 149ks_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 186static 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[]. */ 194void 195x2_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 216being greater than or equal to any number in the sequence. For a 217random real number between zero and one given by a uniformly 218distributed random number generator, this is simply equal to X. */ 219 220static void 221P (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]. */ 228void 229mpf_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 244void 245x2 (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'. */ 287static void 288Pzf (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 309unsigned long int 310mpz_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 */ 362void 363foo_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. */ 377void 378merit (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 429double 430merit_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). */ 443void 444f_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 461void 462vz_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 477void 478spectral_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} 837