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 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 | |
---|