source: trunk/third/gmp/tests/rand/statlib.c @ 18191

Revision 18191, 20.3 KB checked in by ghudson, 22 years ago (diff)
This commit was generated by cvs2svn to compensate for changes in r18190, which included commits to RCS files with non-trunk default branches.
Line 
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.
8
9The GNU MP Library is free software; you can redistribute it and/or modify
10it under the terms of the GNU Lesser General Public License as published by
11the Free Software Foundation; either version 2.1 of the License, or (at your
12option) any later version.
13
14The GNU MP Library is distributed in the hope that it will be useful, but
15WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
16or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
17License for more details.
18
19You should have received a copy of the GNU Lesser General Public License
20along with the GNU MP Library; see the file COPYING.LIB.  If not, write to
21the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
22MA 02111-1307, USA.
23*/
24
25/* The theories for these functions are taken from D. Knuth's "The Art
26of Computer Programming: Volume 2, Seminumerical Algorithms", Third
27Edition, Addison Wesley, 1998. */
28
29/* Implementation notes.
30
31The Kolmogorov-Smirnov test.
32
33Eq. (13) in Knuth, p. 50, says that if X1, X2, ..., Xn are independent
34observations 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
39where F(x) = Pr(X <= x) = probability that (X <= x), which for a
40uniformly distributed random real number between zero and one is
41exactly the number itself (x).
42
43
44The answer to exercise 23 gives the following implementation, which
45doesn't need the observations to be sorted in ascending order:
46
47for (k = 0; k < m; k++)
48        a[k] = 1.0
49        b[k] = 0.0
50        c[k] = 0
51
52for (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
67Kp = sqr (n) * rp
68Km = 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
87void
88ks (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
152void
153ks_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
190static 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[]. */
198void
199x2_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
220being greater than or equal to any number in the sequence.  For a
221random real number between zero and one given by a uniformly
222distributed random number generator, this is simply equal to X. */
223
224static void
225P (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]. */
232void
233mpf_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
248void
249x2 (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'. */
291static void
292Pzf (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
313unsigned long int
314mpz_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 */
366void
367foo_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. */
381void
382merit (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
433double
434merit_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). */
447void
448f_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
465void
466vz_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
481void
482spectral_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
Note: See TracBrowser for help on using the repository browser.