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