xref: /netbsd-src/external/lgpl3/gmp/dist/tests/rand/spect.c (revision 6cd39ddb8550f6fa1bff3fed32053d7f19fd0453)
1 /* spect.c -- the spectral test */
2 
3 /*
4 Copyright 1999 Free Software Foundation, Inc.
5 
6 This file is part of the GNU MP Library test suite.
7 
8 The GNU MP Library test suite is free software; you can redistribute it
9 and/or modify it under the terms of the GNU General Public License as
10 published by the Free Software Foundation; either version 3 of the License,
11 or (at your option) any later version.
12 
13 The GNU MP Library test suite is distributed in the hope that it will be
14 useful, but WITHOUT ANY WARRANTY; without even the implied warranty of
15 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General
16 Public License for more details.
17 
18 You should have received a copy of the GNU General Public License along with
19 the GNU MP Library test suite.  If not, see http://www.gnu.org/licenses/.  */
20 
21 /* T is upper dimension.  Z_A is the LC multiplier, which is
22    relatively prime to Z_M, the LC modulus.  The result is put in
23    rop[] with v[t] in rop[t-2]. */
24 
25 /* BUGS: Due to lazy allocation scheme, maximum T is hard coded to MAXT. */
26 
27 #include <stdio.h>
28 #include <stdlib.h>
29 #include <unistd.h>
30 #include <math.h>
31 #include "gmp.h"
32 
33 #include "gmpstat.h"
34 
35 int g_debug = 0;
36 
37 int
38 main (int argc, char *argv[])
39 {
40   const char usage[] = "usage: spect [-d] a m n\n";
41   int c;
42   unsigned int n;
43   mpz_t a, m;
44   mpf_t res[GMP_SPECT_MAXT], res_min[GMP_SPECT_MAXT], f_tmp;
45   register int f;
46 
47 
48   mpz_init (a);
49   mpz_init (m);
50   for (f = 0; f < GMP_SPECT_MAXT; f++)
51     {
52       mpf_init (res[f]);
53       mpf_init (res_min[f]);
54     }
55   mpf_init (f_tmp);
56   mpf_set_ui (res_min[0], 32768); /* 2^15 */
57   mpf_set_ui (res_min[1], 1024); /* 2^10 */
58   mpf_set_ui (res_min[2], 256); /* 2^8 */
59   mpf_set_ui (res_min[3], 64); /* 2^6 */
60   mpf_set_ui (res_min[4], 32); /* 2^5 */
61 
62   while ((c = getopt (argc, argv, "dh")) != -1)
63     switch (c)
64       {
65       case 'd':			/* debug */
66 	g_debug++;
67 	break;
68       case 'h':
69       default:
70 	fputs (usage, stderr);
71 	exit (1);
72       }
73   argc -= optind;
74   argv += optind;
75 
76   if (argc < 3)
77     {
78       fputs (usage, stderr);
79       exit (1);
80     }
81 
82   mpz_set_str (a, argv[0], 0);
83   mpz_set_str (m, argv[1], 0);
84   n = (unsigned int) atoi (argv[2]);
85   if (n + 1 > GMP_SPECT_MAXT)
86     n = GMP_SPECT_MAXT + 1;
87 
88   spectral_test (res, n, a, m);
89 
90   for (f = 0; f < n - 1; f++)
91     {
92       /* print v */
93       printf ("%d: v = ", f + 2);
94       mpf_out_str (stdout, 10, 4, res[f]);
95 
96 #ifdef PRINT_RAISED_BY_TWO_AS_WELL
97       printf (" (^2 = ");
98       mpf_mul (f_tmp, res[f], res[f]);
99       mpf_out_str (stdout, 10, 4, f_tmp);
100       printf (")");
101 #endif /* PRINT_RAISED_BY_TWO_AS_WELL */
102 
103       /* print merit */
104       printf (" m = ");
105       merit (f_tmp, f + 2, res[f], m);
106       mpf_out_str (stdout, 10, 4, f_tmp);
107 
108       if (mpf_cmp (res[f], res_min[f]) < 0)
109 	printf ("\t*** v too low ***");
110       if (mpf_get_d (f_tmp) < .1)
111 	printf ("\t*** merit too low ***");
112 
113       puts ("");
114     }
115 
116   mpz_clear (a);
117   mpz_clear (m);
118   for (f = 0; f < GMP_SPECT_MAXT; f++)
119     {
120       mpf_clear (res[f]);
121       mpf_clear (res_min[f]);
122     }
123   mpf_clear (f_tmp);
124 
125   return 0;
126 }
127 
128 
129 void
130 debug_foo()
131 {
132   if (0)
133     {
134       mpz_dump (0);
135       mpf_dump (0);
136     }
137 }
138