diff options
Diffstat (limited to 'vendor/gmp-6.3.0/tests/rand/statlib.c')
-rw-r--r-- | vendor/gmp-6.3.0/tests/rand/statlib.c | 836 |
1 files changed, 836 insertions, 0 deletions
diff --git a/vendor/gmp-6.3.0/tests/rand/statlib.c b/vendor/gmp-6.3.0/tests/rand/statlib.c new file mode 100644 index 0000000..db05380 --- /dev/null +++ b/vendor/gmp-6.3.0/tests/rand/statlib.c @@ -0,0 +1,836 @@ +/* statlib.c -- Statistical functions for testing the randomness of + number sequences. */ + +/* +Copyright 1999, 2000 Free Software Foundation, Inc. + +This file is part of the GNU MP Library test suite. + +The GNU MP Library test suite is free software; you can redistribute it +and/or modify it under the terms of the GNU General Public License as +published by the Free Software Foundation; either version 3 of the License, +or (at your option) any later version. + +The GNU MP Library test suite is distributed in the hope that it will be +useful, but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General +Public License for more details. + +You should have received a copy of the GNU General Public License along with +the GNU MP Library test suite. If not, see https://www.gnu.org/licenses/. */ + +/* The theories for these functions are taken from D. Knuth's "The Art +of Computer Programming: Volume 2, Seminumerical Algorithms", Third +Edition, Addison Wesley, 1998. */ + +/* Implementation notes. + +The Kolmogorov-Smirnov test. + +Eq. (13) in Knuth, p. 50, says that if X1, X2, ..., Xn are independent +observations arranged into ascending order + + Kp = sqr(n) * max(j/n - F(Xj)) for all 1<=j<=n + Km = sqr(n) * max(F(Xj) - (j-1)/n)) for all 1<=j<=n + +where F(x) = Pr(X <= x) = probability that (X <= x), which for a +uniformly distributed random real number between zero and one is +exactly the number itself (x). + + +The answer to exercise 23 gives the following implementation, which +doesn't need the observations to be sorted in ascending order: + +for (k = 0; k < m; k++) + a[k] = 1.0 + b[k] = 0.0 + c[k] = 0 + +for (each observation Xj) + Y = F(Xj) + k = floor (m * Y) + a[k] = min (a[k], Y) + b[k] = max (b[k], Y) + c[k] += 1 + + j = 0 + rp = rm = 0 + for (k = 0; k < m; k++) + if (c[k] > 0) + rm = max (rm, a[k] - j/n) + j += c[k] + rp = max (rp, j/n - b[k]) + +Kp = sqr (n) * rp +Km = sqr (n) * rm + +*/ + +#include <stdio.h> +#include <stdlib.h> +#include <math.h> + +#include "gmpstat.h" + +/* ks (Kp, Km, X, P, n) -- Perform a Kolmogorov-Smirnov test on the N + real numbers between zero and one in vector X. P is the + distribution function, called for each entry in X, which should + calculate the probability of X being greater than or equal to any + number in the sequence. (For a uniformly distributed sequence of + real numbers between zero and one, this is simply equal to X.) The + result is put in Kp and Km. */ + +void +ks (mpf_t Kp, + mpf_t Km, + mpf_t X[], + void (P) (mpf_t, mpf_t), + unsigned long int n) +{ + mpf_t Kt; /* temp */ + mpf_t f_x; + mpf_t f_j; /* j */ + mpf_t f_jnq; /* j/n or (j-1)/n */ + unsigned long int j; + + /* Sort the vector in ascending order. */ + qsort (X, n, sizeof (__mpf_struct), mpf_cmp); + + /* K-S test. */ + /* Kp = sqr(n) * max(j/n - F(Xj)) for all 1<=j<=n + Km = sqr(n) * max(F(Xj) - (j-1)/n)) for all 1<=j<=n + */ + + mpf_init (Kt); mpf_init (f_x); mpf_init (f_j); mpf_init (f_jnq); + mpf_set_ui (Kp, 0); mpf_set_ui (Km, 0); + for (j = 1; j <= n; j++) + { + P (f_x, X[j-1]); + mpf_set_ui (f_j, j); + + mpf_div_ui (f_jnq, f_j, n); + mpf_sub (Kt, f_jnq, f_x); + if (mpf_cmp (Kt, Kp) > 0) + mpf_set (Kp, Kt); + if (g_debug > DEBUG_2) + { + printf ("j=%lu ", j); + printf ("P()="); mpf_out_str (stdout, 10, 2, f_x); printf ("\t"); + + printf ("jnq="); mpf_out_str (stdout, 10, 2, f_jnq); printf (" "); + printf ("diff="); mpf_out_str (stdout, 10, 2, Kt); printf (" "); + printf ("Kp="); mpf_out_str (stdout, 10, 2, Kp); printf ("\t"); + } + mpf_sub_ui (f_j, f_j, 1); + mpf_div_ui (f_jnq, f_j, n); + mpf_sub (Kt, f_x, f_jnq); + if (mpf_cmp (Kt, Km) > 0) + mpf_set (Km, Kt); + + if (g_debug > DEBUG_2) + { + printf ("jnq="); mpf_out_str (stdout, 10, 2, f_jnq); printf (" "); + printf ("diff="); mpf_out_str (stdout, 10, 2, Kt); printf (" "); + printf ("Km="); mpf_out_str (stdout, 10, 2, Km); printf (" "); + printf ("\n"); + } + } + mpf_sqrt_ui (Kt, n); + mpf_mul (Kp, Kp, Kt); + mpf_mul (Km, Km, Kt); + + mpf_clear (Kt); mpf_clear (f_x); mpf_clear (f_j); mpf_clear (f_jnq); +} + +/* ks_table(val, n) -- calculate probability for Kp/Km less than or + equal to VAL with N observations. See [Knuth section 3.3.1] */ + +void +ks_table (mpf_t p, mpf_t val, const unsigned int n) +{ + /* We use Eq. (27), Knuth p.58, skipping O(1/n) for simplicity. + This shortcut will result in too high probabilities, especially + when n is small. + + Pr(Kp(n) <= s) = 1 - pow(e, -2*s^2) * (1 - 2/3*s/sqrt(n) + O(1/n)) */ + + /* We have 's' in variable VAL and store the result in P. */ + + mpf_t t1, t2; + + mpf_init (t1); mpf_init (t2); + + /* t1 = 1 - 2/3 * s/sqrt(n) */ + mpf_sqrt_ui (t1, n); + mpf_div (t1, val, t1); + mpf_mul_ui (t1, t1, 2); + mpf_div_ui (t1, t1, 3); + mpf_ui_sub (t1, 1, t1); + + /* t2 = pow(e, -2*s^2) */ +#ifndef OLDGMP + mpf_pow_ui (t2, val, 2); /* t2 = s^2 */ + mpf_set_d (t2, exp (-(2.0 * mpf_get_d (t2)))); +#else + /* hmmm, gmp doesn't have pow() for floats. use doubles. */ + mpf_set_d (t2, pow (M_E, -(2 * pow (mpf_get_d (val), 2)))); +#endif + + /* p = 1 - t1 * t2 */ + mpf_mul (t1, t1, t2); + mpf_ui_sub (p, 1, t1); + + mpf_clear (t1); mpf_clear (t2); +} + +static double x2_table_X[][7] = { + { -2.33, -1.64, -.674, 0.0, 0.674, 1.64, 2.33 }, /* x */ + { 5.4289, 2.6896, .454276, 0.0, .454276, 2.6896, 5.4289} /* x^2 */ +}; + +#define _2D3 ((double) .6666666666) + +/* x2_table (t, v, n) -- return chi-square table row for V in T[]. */ +void +x2_table (double t[], + unsigned int v) +{ + int f; + + + /* FIXME: Do a table lookup for v <= 30 since the following formula + [Knuth, vol 2, 3.3.1] is only good for v > 30. */ + + /* value = v + sqrt(2*v) * X[p] + (2/3) * X[p]^2 - 2/3 + O(1/sqrt(t) */ + /* NOTE: The O() term is ignored for simplicity. */ + + for (f = 0; f < 7; f++) + t[f] = + v + + sqrt (2 * v) * x2_table_X[0][f] + + _2D3 * x2_table_X[1][f] - _2D3; +} + + +/* P(p, x) -- Distribution function. Calculate the probability of X +being greater than or equal to any number in the sequence. For a +random real number between zero and one given by a uniformly +distributed random number generator, this is simply equal to X. */ + +static void +P (mpf_t p, mpf_t x) +{ + mpf_set (p, x); +} + +/* mpf_freqt() -- Frequency test using KS on N real numbers between zero + and one. See [Knuth vol 2, p.61]. */ +void +mpf_freqt (mpf_t Kp, + mpf_t Km, + mpf_t X[], + const unsigned long int n) +{ + ks (Kp, Km, X, P, n); +} + + +/* The Chi-square test. Eq. (8) in Knuth vol. 2 says that if Y[] + holds the observations and p[] is the probability for.. (to be + continued!) + + V = 1/n * sum((s=1 to k) Y[s]^2 / p[s]) - n */ + +void +x2 (mpf_t V, /* result */ + unsigned long int X[], /* data */ + unsigned int k, /* #of categories */ + void (P) (mpf_t, unsigned long int, void *), /* probability func */ + void *x, /* extra user data passed to P() */ + unsigned long int n) /* #of samples */ +{ + unsigned int f; + mpf_t f_t, f_t2; /* temp floats */ + + mpf_init (f_t); mpf_init (f_t2); + + + mpf_set_ui (V, 0); + for (f = 0; f < k; f++) + { + if (g_debug > DEBUG_2) + fprintf (stderr, "%u: P()=", f); + mpf_set_ui (f_t, X[f]); + mpf_mul (f_t, f_t, f_t); /* f_t = X[f]^2 */ + P (f_t2, f, x); /* f_t2 = Pr(f) */ + if (g_debug > DEBUG_2) + mpf_out_str (stderr, 10, 2, f_t2); + mpf_div (f_t, f_t, f_t2); + mpf_add (V, V, f_t); + if (g_debug > DEBUG_2) + { + fprintf (stderr, "\tV="); + mpf_out_str (stderr, 10, 2, V); + fprintf (stderr, "\t"); + } + } + if (g_debug > DEBUG_2) + fprintf (stderr, "\n"); + mpf_div_ui (V, V, n); + mpf_sub_ui (V, V, n); + + mpf_clear (f_t); mpf_clear (f_t2); +} + +/* Pzf(p, s, x) -- Probability for category S in mpz_freqt(). It's + 1/d for all S. X is a pointer to an unsigned int holding 'd'. */ +static void +Pzf (mpf_t p, unsigned long int s, void *x) +{ + mpf_set_ui (p, 1); + mpf_div_ui (p, p, *((unsigned int *) x)); +} + +/* mpz_freqt(V, X, imax, n) -- Frequency test on integers. [Knuth, + vol 2, 3.3.2]. Keep IMAX low on this one, since we loop from 0 to + IMAX. 128 or 256 could be nice. + + X[] must not contain numbers outside the range 0 <= X <= IMAX. + + Return value is number of observations actually used, after + discarding entries out of range. + + Since X[] contains integers between zero and IMAX, inclusive, we + have IMAX+1 categories. + + Note that N should be at least 5*IMAX. Result is put in V and can + be compared to output from x2_table (v=IMAX). */ + +unsigned long int +mpz_freqt (mpf_t V, + mpz_t X[], + unsigned int imax, + const unsigned long int n) +{ + unsigned long int *v; /* result */ + unsigned int f; + unsigned int d; /* number of categories = imax+1 */ + unsigned int uitemp; + unsigned long int usedn; + + + d = imax + 1; + + v = (unsigned long int *) calloc (imax + 1, sizeof (unsigned long int)); + if (NULL == v) + { + fprintf (stderr, "mpz_freqt(): out of memory\n"); + exit (1); + } + + /* count */ + usedn = n; /* actual number of observations */ + for (f = 0; f < n; f++) + { + uitemp = mpz_get_ui(X[f]); + if (uitemp > imax) /* sanity check */ + { + if (g_debug) + fprintf (stderr, "mpz_freqt(): warning: input insanity: %u, "\ + "ignored.\n", uitemp); + usedn--; + continue; + } + v[uitemp]++; + } + + if (g_debug > DEBUG_2) + { + fprintf (stderr, "counts:\n"); + for (f = 0; f <= imax; f++) + fprintf (stderr, "%u:\t%lu\n", f, v[f]); + } + + /* chi-square with k=imax+1 and P(x)=1/(imax+1) for all x.*/ + x2 (V, v, d, Pzf, (void *) &d, usedn); + + free (v); + return (usedn); +} + +/* debug dummy to drag in dump funcs */ +void +foo_debug () +{ + if (0) + { + mpf_dump (0); +#ifndef OLDGMP + mpz_dump (0); +#endif + } +} + +/* merit (rop, t, v, m) -- calculate merit for spectral test result in + dimension T, see Knuth p. 105. BUGS: Only valid for 2 <= T <= + 6. */ +void +merit (mpf_t rop, unsigned int t, mpf_t v, mpz_t m) +{ + int f; + mpf_t f_m, f_const, f_pi; + + mpf_init (f_m); + mpf_set_z (f_m, m); + mpf_init_set_d (f_const, M_PI); + mpf_init_set_d (f_pi, M_PI); + + switch (t) + { + case 2: /* PI */ + break; + case 3: /* PI * 4/3 */ + mpf_mul_ui (f_const, f_const, 4); + mpf_div_ui (f_const, f_const, 3); + break; + case 4: /* PI^2 * 1/2 */ + mpf_mul (f_const, f_const, f_pi); + mpf_div_ui (f_const, f_const, 2); + break; + case 5: /* PI^2 * 8/15 */ + mpf_mul (f_const, f_const, f_pi); + mpf_mul_ui (f_const, f_const, 8); + mpf_div_ui (f_const, f_const, 15); + break; + case 6: /* PI^3 * 1/6 */ + mpf_mul (f_const, f_const, f_pi); + mpf_mul (f_const, f_const, f_pi); + mpf_div_ui (f_const, f_const, 6); + break; + default: + fprintf (stderr, + "spect (merit): can't calculate merit for dimensions > 6\n"); + mpf_set_ui (f_const, 0); + break; + } + + /* rop = v^t */ + mpf_set (rop, v); + for (f = 1; f < t; f++) + mpf_mul (rop, rop, v); + mpf_mul (rop, rop, f_const); + mpf_div (rop, rop, f_m); + + mpf_clear (f_m); + mpf_clear (f_const); + mpf_clear (f_pi); +} + +double +merit_u (unsigned int t, mpf_t v, mpz_t m) +{ + mpf_t rop; + double res; + + mpf_init (rop); + merit (rop, t, v, m); + res = mpf_get_d (rop); + mpf_clear (rop); + return res; +} + +/* f_floor (rop, op) -- Set rop = floor (op). */ +void +f_floor (mpf_t rop, mpf_t op) +{ + mpz_t z; + + mpz_init (z); + + /* No mpf_floor(). Convert to mpz and back. */ + mpz_set_f (z, op); + mpf_set_z (rop, z); + + mpz_clear (z); +} + + +/* vz_dot (rop, v1, v2, nelem) -- compute dot product of z-vectors V1, + V2. N is number of elements in vectors V1 and V2. */ + +void +vz_dot (mpz_t rop, mpz_t V1[], mpz_t V2[], unsigned int n) +{ + mpz_t t; + + mpz_init (t); + mpz_set_ui (rop, 0); + while (n--) + { + mpz_mul (t, V1[n], V2[n]); + mpz_add (rop, rop, t); + } + + mpz_clear (t); +} + +void +spectral_test (mpf_t rop[], unsigned int T, mpz_t a, mpz_t m) +{ + /* Knuth "Seminumerical Algorithms, Third Edition", section 3.3.4 + (pp. 101-103). */ + + /* v[t] = min { sqrt (x[1]^2 + ... + x[t]^2) | + x[1] + a*x[2] + ... + pow (a, t-1) * x[t] is congruent to 0 (mod m) } */ + + + /* Variables. */ + unsigned int ui_t; + unsigned int ui_i, ui_j, ui_k, ui_l; + mpf_t f_tmp1, f_tmp2; + mpz_t tmp1, tmp2, tmp3; + mpz_t U[GMP_SPECT_MAXT][GMP_SPECT_MAXT], + V[GMP_SPECT_MAXT][GMP_SPECT_MAXT], + X[GMP_SPECT_MAXT], + Y[GMP_SPECT_MAXT], + Z[GMP_SPECT_MAXT]; + mpz_t h, hp, r, s, p, pp, q, u, v; + + /* GMP inits. */ + mpf_init (f_tmp1); + mpf_init (f_tmp2); + for (ui_i = 0; ui_i < GMP_SPECT_MAXT; ui_i++) + { + for (ui_j = 0; ui_j < GMP_SPECT_MAXT; ui_j++) + { + mpz_init_set_ui (U[ui_i][ui_j], 0); + mpz_init_set_ui (V[ui_i][ui_j], 0); + } + mpz_init_set_ui (X[ui_i], 0); + mpz_init_set_ui (Y[ui_i], 0); + mpz_init (Z[ui_i]); + } + mpz_init (tmp1); + mpz_init (tmp2); + mpz_init (tmp3); + mpz_init (h); + mpz_init (hp); + mpz_init (r); + mpz_init (s); + mpz_init (p); + mpz_init (pp); + mpz_init (q); + mpz_init (u); + mpz_init (v); + + /* Implementation inits. */ + if (T > GMP_SPECT_MAXT) + T = GMP_SPECT_MAXT; /* FIXME: Lazy. */ + + /* S1 [Initialize.] */ + ui_t = 2 - 1; /* NOTE: `t' in description == ui_t + 1 + for easy indexing */ + mpz_set (h, a); + mpz_set (hp, m); + mpz_set_ui (p, 1); + mpz_set_ui (pp, 0); + mpz_set (r, a); + mpz_pow_ui (s, a, 2); + mpz_add_ui (s, s, 1); /* s = 1 + a^2 */ + + /* S2 [Euclidean step.] */ + while (1) + { + if (g_debug > DEBUG_1) + { + mpz_mul (tmp1, h, pp); + mpz_mul (tmp2, hp, p); + mpz_sub (tmp1, tmp1, tmp2); + if (mpz_cmpabs (m, tmp1)) + { + printf ("***BUG***: h*pp - hp*p = "); + mpz_out_str (stdout, 10, tmp1); + printf ("\n"); + } + } + if (g_debug > DEBUG_2) + { + printf ("hp = "); + mpz_out_str (stdout, 10, hp); + printf ("\nh = "); + mpz_out_str (stdout, 10, h); + printf ("\n"); + fflush (stdout); + } + + if (mpz_sgn (h)) + mpz_tdiv_q (q, hp, h); /* q = floor(hp/h) */ + else + mpz_set_ui (q, 1); + + if (g_debug > DEBUG_2) + { + printf ("q = "); + mpz_out_str (stdout, 10, q); + printf ("\n"); + fflush (stdout); + } + + mpz_mul (tmp1, q, h); + mpz_sub (u, hp, tmp1); /* u = hp - q*h */ + + mpz_mul (tmp1, q, p); + mpz_sub (v, pp, tmp1); /* v = pp - q*p */ + + mpz_pow_ui (tmp1, u, 2); + mpz_pow_ui (tmp2, v, 2); + mpz_add (tmp1, tmp1, tmp2); + if (mpz_cmp (tmp1, s) < 0) + { + mpz_set (s, tmp1); /* s = u^2 + v^2 */ + mpz_set (hp, h); /* hp = h */ + mpz_set (h, u); /* h = u */ + mpz_set (pp, p); /* pp = p */ + mpz_set (p, v); /* p = v */ + } + else + break; + } + + /* S3 [Compute v2.] */ + mpz_sub (u, u, h); + mpz_sub (v, v, p); + + mpz_pow_ui (tmp1, u, 2); + mpz_pow_ui (tmp2, v, 2); + mpz_add (tmp1, tmp1, tmp2); + if (mpz_cmp (tmp1, s) < 0) + { + mpz_set (s, tmp1); /* s = u^2 + v^2 */ + mpz_set (hp, u); + mpz_set (pp, v); + } + mpf_set_z (f_tmp1, s); + mpf_sqrt (rop[ui_t - 1], f_tmp1); + + /* S4 [Advance t.] */ + mpz_neg (U[0][0], h); + mpz_set (U[0][1], p); + mpz_neg (U[1][0], hp); + mpz_set (U[1][1], pp); + + mpz_set (V[0][0], pp); + mpz_set (V[0][1], hp); + mpz_neg (V[1][0], p); + mpz_neg (V[1][1], h); + if (mpz_cmp_ui (pp, 0) > 0) + { + mpz_neg (V[0][0], V[0][0]); + mpz_neg (V[0][1], V[0][1]); + mpz_neg (V[1][0], V[1][0]); + mpz_neg (V[1][1], V[1][1]); + } + + while (ui_t + 1 != T) /* S4 loop */ + { + ui_t++; + mpz_mul (r, a, r); + mpz_mod (r, r, m); + + /* Add new row and column to U and V. They are initialized with + all elements set to zero, so clearing is not necessary. */ + + mpz_neg (U[ui_t][0], r); /* U: First col in new row. */ + mpz_set_ui (U[ui_t][ui_t], 1); /* U: Last col in new row. */ + + mpz_set (V[ui_t][ui_t], m); /* V: Last col in new row. */ + + /* "Finally, for 1 <= i < t, + set q = round (vi1 * r / m), + vit = vi1*r - q*m, + and Ut=Ut+q*Ui */ + + for (ui_i = 0; ui_i < ui_t; ui_i++) + { + mpz_mul (tmp1, V[ui_i][0], r); /* tmp1=vi1*r */ + zdiv_round (q, tmp1, m); /* q=round(vi1*r/m) */ + mpz_mul (tmp2, q, m); /* tmp2=q*m */ + mpz_sub (V[ui_i][ui_t], tmp1, tmp2); + + for (ui_j = 0; ui_j <= ui_t; ui_j++) /* U[t] = U[t] + q*U[i] */ + { + mpz_mul (tmp1, q, U[ui_i][ui_j]); /* tmp=q*uij */ + mpz_add (U[ui_t][ui_j], U[ui_t][ui_j], tmp1); /* utj = utj + q*uij */ + } + } + + /* s = min (s, zdot (U[t], U[t]) */ + vz_dot (tmp1, U[ui_t], U[ui_t], ui_t + 1); + if (mpz_cmp (tmp1, s) < 0) + mpz_set (s, tmp1); + + ui_k = ui_t; + ui_j = 0; /* WARNING: ui_j no longer a temp. */ + + /* S5 [Transform.] */ + if (g_debug > DEBUG_2) + printf ("(t, k, j, q1, q2, ...)\n"); + do + { + if (g_debug > DEBUG_2) + printf ("(%u, %u, %u", ui_t + 1, ui_k + 1, ui_j + 1); + + for (ui_i = 0; ui_i <= ui_t; ui_i++) + { + if (ui_i != ui_j) + { + vz_dot (tmp1, V[ui_i], V[ui_j], ui_t + 1); /* tmp1=dot(Vi,Vj). */ + mpz_abs (tmp2, tmp1); + mpz_mul_ui (tmp2, tmp2, 2); /* tmp2 = 2*abs(dot(Vi,Vj) */ + vz_dot (tmp3, V[ui_j], V[ui_j], ui_t + 1); /* tmp3=dot(Vj,Vj). */ + + if (mpz_cmp (tmp2, tmp3) > 0) + { + zdiv_round (q, tmp1, tmp3); /* q=round(Vi.Vj/Vj.Vj) */ + if (g_debug > DEBUG_2) + { + printf (", "); + mpz_out_str (stdout, 10, q); + } + + for (ui_l = 0; ui_l <= ui_t; ui_l++) + { + mpz_mul (tmp1, q, V[ui_j][ui_l]); + mpz_sub (V[ui_i][ui_l], V[ui_i][ui_l], tmp1); /* Vi=Vi-q*Vj */ + mpz_mul (tmp1, q, U[ui_i][ui_l]); + mpz_add (U[ui_j][ui_l], U[ui_j][ui_l], tmp1); /* Uj=Uj+q*Ui */ + } + + vz_dot (tmp1, U[ui_j], U[ui_j], ui_t + 1); /* tmp1=dot(Uj,Uj) */ + if (mpz_cmp (tmp1, s) < 0) /* s = min(s,dot(Uj,Uj)) */ + mpz_set (s, tmp1); + ui_k = ui_j; + } + else if (g_debug > DEBUG_2) + printf (", #"); /* 2|Vi.Vj| <= Vj.Vj */ + } + else if (g_debug > DEBUG_2) + printf (", *"); /* i == j */ + } + + if (g_debug > DEBUG_2) + printf (")\n"); + + /* S6 [Advance j.] */ + if (ui_j == ui_t) + ui_j = 0; + else + ui_j++; + } + while (ui_j != ui_k); /* S5 */ + + /* From Knuth p. 104: "The exhaustive search in steps S8-S10 + reduces the value of s only rarely." */ +#ifdef DO_SEARCH + /* S7 [Prepare for search.] */ + /* Find minimum in (x[1], ..., x[t]) satisfying condition + x[k]^2 <= f(y[1], ...,y[t]) * dot(V[k],V[k]) */ + + ui_k = ui_t; + if (g_debug > DEBUG_2) + { + printf ("searching..."); + /*for (f = 0; f < ui_t*/ + fflush (stdout); + } + + /* Z[i] = floor (sqrt (floor (dot(V[i],V[i]) * s / m^2))); */ + mpz_pow_ui (tmp1, m, 2); + mpf_set_z (f_tmp1, tmp1); + mpf_set_z (f_tmp2, s); + mpf_div (f_tmp1, f_tmp2, f_tmp1); /* f_tmp1 = s/m^2 */ + for (ui_i = 0; ui_i <= ui_t; ui_i++) + { + vz_dot (tmp1, V[ui_i], V[ui_i], ui_t + 1); + mpf_set_z (f_tmp2, tmp1); + mpf_mul (f_tmp2, f_tmp2, f_tmp1); + f_floor (f_tmp2, f_tmp2); + mpf_sqrt (f_tmp2, f_tmp2); + mpz_set_f (Z[ui_i], f_tmp2); + } + + /* S8 [Advance X[k].] */ + do + { + if (g_debug > DEBUG_2) + { + printf ("X[%u] = ", ui_k); + mpz_out_str (stdout, 10, X[ui_k]); + printf ("\tZ[%u] = ", ui_k); + mpz_out_str (stdout, 10, Z[ui_k]); + printf ("\n"); + fflush (stdout); + } + + if (mpz_cmp (X[ui_k], Z[ui_k])) + { + mpz_add_ui (X[ui_k], X[ui_k], 1); + for (ui_i = 0; ui_i <= ui_t; ui_i++) + mpz_add (Y[ui_i], Y[ui_i], U[ui_k][ui_i]); + + /* S9 [Advance k.] */ + while (++ui_k <= ui_t) + { + mpz_neg (X[ui_k], Z[ui_k]); + mpz_mul_ui (tmp1, Z[ui_k], 2); + for (ui_i = 0; ui_i <= ui_t; ui_i++) + { + mpz_mul (tmp2, tmp1, U[ui_k][ui_i]); + mpz_sub (Y[ui_i], Y[ui_i], tmp2); + } + } + vz_dot (tmp1, Y, Y, ui_t + 1); + if (mpz_cmp (tmp1, s) < 0) + mpz_set (s, tmp1); + } + } + while (--ui_k); +#endif /* DO_SEARCH */ + mpf_set_z (f_tmp1, s); + mpf_sqrt (rop[ui_t - 1], f_tmp1); +#ifdef DO_SEARCH + if (g_debug > DEBUG_2) + printf ("done.\n"); +#endif /* DO_SEARCH */ + } /* S4 loop */ + + /* Clear GMP variables. */ + + mpf_clear (f_tmp1); + mpf_clear (f_tmp2); + for (ui_i = 0; ui_i < GMP_SPECT_MAXT; ui_i++) + { + for (ui_j = 0; ui_j < GMP_SPECT_MAXT; ui_j++) + { + mpz_clear (U[ui_i][ui_j]); + mpz_clear (V[ui_i][ui_j]); + } + mpz_clear (X[ui_i]); + mpz_clear (Y[ui_i]); + mpz_clear (Z[ui_i]); + } + mpz_clear (tmp1); + mpz_clear (tmp2); + mpz_clear (tmp3); + mpz_clear (h); + mpz_clear (hp); + mpz_clear (r); + mpz_clear (s); + mpz_clear (p); + mpz_clear (pp); + mpz_clear (q); + mpz_clear (u); + mpz_clear (v); + + return; +} |