Annotation of OpenXM_contrib/gmp/tests/rand/statlib.c, Revision 1.1.1.2
1.1 maekawa 1: /* statlib.c -- Statistical functions for testing the randomness of
2: number sequences. */
3:
4: /*
1.1.1.2 ! ohara 5: Copyright 1999, 2000 Free Software Foundation, Inc.
1.1 maekawa 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])
1.1.1.2 ! ohara 66:
1.1 maekawa 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);
1.1.1.2 ! ohara 542: mpz_pow_ui (s, a, 2);
1.1 maekawa 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." */
1.1.1.2 ! ohara 738: #ifdef DO_SEARCH
1.1 maekawa 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);
1.1.1.2 ! ohara 805: #ifdef DO_SEARCH
1.1 maekawa 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