1 /*
2 Copyright 2000 Free Software Foundation, Inc.
3
4 This file is part of the GNU MP Library test suite.
5
6 The GNU MP Library test suite is free software; you can redistribute it
7 and/or modify it under the terms of the GNU General Public License as
8 published by the Free Software Foundation; either version 3 of the License,
9 or (at your option) any later version.
10
11 The GNU MP Library test suite is distributed in the hope that it will be
12 useful, but WITHOUT ANY WARRANTY; without even the implied warranty of
13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
14 Public License for more details.
15
16 You should have received a copy of the GNU General Public License along with
17 the GNU MP Library test suite. If not, see https://www.gnu.org/licenses/. */
18
19 #include <stdio.h>
20 #include <stdlib.h>
21 #include <unistd.h>
22 #include <signal.h>
23 #include <math.h>
24 #include "gmpstat.h"
25
26 #define RCSID(msg) \
27 static /**/const char *const rcsid[] = { (char *)rcsid, "\100(#)" msg }
28
29 RCSID("$Id: findlc.c,v 1.1.1.4 2020/09/27 00:27:04 mrg Exp $");
30
31 int g_debug = 0;
32
33 static mpz_t a;
34
35 static void
sh_status(int sig)36 sh_status (int sig)
37 {
38 printf ("sh_status: signal %d caught. dumping status.\n", sig);
39
40 printf (" a = ");
41 mpz_out_str (stdout, 10, a);
42 printf ("\n");
43 fflush (stdout);
44
45 if (SIGSEGV == sig) /* remove SEGV handler */
46 signal (SIGSEGV, SIG_DFL);
47 }
48
49 /* Input is a modulus (m). We shall find multiplier (a) and adder (c)
50 conforming to the rules found in the first comment block in file
51 mpz/urandom.c.
52
53 Then run a spectral test on the generator and discard any
54 multipliers not passing. */
55
56 /* TODO:
57
58 . find a better algorithm than a+=8; bigger jumps perhaps?
59
60 */
61
62 void
mpz_true_random(mpz_t s,unsigned long int nbits)63 mpz_true_random (mpz_t s, unsigned long int nbits)
64 {
65 #if __FreeBSD__
66 FILE *fs;
67 char c[1];
68 int i;
69
70 mpz_set_ui (s, 0);
71 for (i = 0; i < nbits; i += 8)
72 {
73 for (;;)
74 {
75 int nread;
76 fs = fopen ("/dev/random", "r");
77 nread = fread (c, 1, 1, fs);
78 fclose (fs);
79 if (nread != 0)
80 break;
81 sleep (1);
82 }
83 mpz_mul_2exp (s, s, 8);
84 mpz_add_ui (s, s, ((unsigned long int) c[0]) & 0xff);
85 printf ("%d random bits\n", i + 8);
86 }
87 if (nbits % 8 != 0)
88 mpz_mod_2exp (s, s, nbits);
89 #endif
90 }
91
92 int
main(int argc,char * argv[])93 main (int argc, char *argv[])
94 {
95 const char usage[] = "usage: findlc [-dv] m2exp [low_merit [high_merit]]\n";
96 int f;
97 int v_lose, m_lose, v_best, m_best;
98 int c;
99 int debug = 1;
100 int cnt_high_merit;
101 mpz_t m;
102 unsigned long int m2exp;
103 #define DIMS 6 /* dimensions run in spectral test */
104 mpf_t v[DIMS-1]; /* spectral test result (there's no v
105 for 1st dimension */
106 mpf_t f_merit, low_merit, high_merit;
107 mpz_t acc, minus8;
108 mpz_t min, max;
109 mpz_t s;
110
111
112 mpz_init (m);
113 mpz_init (a);
114 for (f = 0; f < DIMS-1; f++)
115 mpf_init (v[f]);
116 mpf_init (f_merit);
117 mpf_init_set_d (low_merit, .1);
118 mpf_init_set_d (high_merit, .1);
119
120 while ((c = getopt (argc, argv, "a:di:hv")) != -1)
121 switch (c)
122 {
123 case 'd': /* debug */
124 g_debug++;
125 break;
126
127 case 'v': /* print version */
128 puts (rcsid[1]);
129 exit (0);
130
131 case 'h':
132 case '?':
133 default:
134 fputs (usage, stderr);
135 exit (1);
136 }
137
138 argc -= optind;
139 argv += optind;
140
141 if (argc < 1)
142 {
143 fputs (usage, stderr);
144 exit (1);
145 }
146
147 /* Install signal handler. */
148 if (SIG_ERR == signal (SIGSEGV, sh_status))
149 {
150 perror ("signal (SIGSEGV)");
151 exit (1);
152 }
153 if (SIG_ERR == signal (SIGHUP, sh_status))
154 {
155 perror ("signal (SIGHUP)");
156 exit (1);
157 }
158
159 printf ("findlc: version: %s\n", rcsid[1]);
160 m2exp = atol (argv[0]);
161 mpz_init_set_ui (m, 1);
162 mpz_mul_2exp (m, m, m2exp);
163 printf ("m = 0x");
164 mpz_out_str (stdout, 16, m);
165 puts ("");
166
167 if (argc > 1) /* have low_merit */
168 mpf_set_str (low_merit, argv[1], 0);
169 if (argc > 2) /* have high_merit */
170 mpf_set_str (high_merit, argv[2], 0);
171
172 if (debug)
173 {
174 fprintf (stderr, "low_merit = ");
175 mpf_out_str (stderr, 10, 2, low_merit);
176 fprintf (stderr, "; high_merit = ");
177 mpf_out_str (stderr, 10, 2, high_merit);
178 fputs ("\n", stderr);
179 }
180
181 mpz_init (minus8);
182 mpz_set_si (minus8, -8L);
183 mpz_init_set_ui (acc, 0);
184 mpz_init (s);
185 mpz_init_set_d (min, 0.01 * pow (2.0, (double) m2exp));
186 mpz_init_set_d (max, 0.99 * pow (2.0, (double) m2exp));
187
188 mpz_true_random (s, m2exp); /* Start. */
189 mpz_setbit (s, 0); /* Make it odd. */
190
191 v_best = m_best = 2*(DIMS-1);
192 for (;;)
193 {
194 mpz_add (acc, acc, s);
195 mpz_mod_2exp (acc, acc, m2exp);
196 #if later
197 mpz_and_si (a, acc, -8L);
198 #else
199 mpz_and (a, acc, minus8);
200 #endif
201 mpz_add_ui (a, a, 5);
202 if (mpz_cmp (a, min) <= 0 || mpz_cmp (a, max) >= 0)
203 continue;
204
205 spectral_test (v, DIMS, a, m);
206 for (f = 0, v_lose = m_lose = 0, cnt_high_merit = DIMS-1;
207 f < DIMS-1; f++)
208 {
209 merit (f_merit, f + 2, v[f], m);
210
211 if (mpf_cmp_ui (v[f], 1 << (30 / (f + 2) + (f == 2))) < 0)
212 v_lose++;
213
214 if (mpf_cmp (f_merit, low_merit) < 0)
215 m_lose++;
216
217 if (mpf_cmp (f_merit, high_merit) >= 0)
218 cnt_high_merit--;
219 }
220
221 if (0 == v_lose && 0 == m_lose)
222 {
223 mpz_out_str (stdout, 10, a); puts (""); fflush (stdout);
224 if (0 == cnt_high_merit)
225 break; /* leave loop */
226 }
227 if (v_lose < v_best)
228 {
229 v_best = v_lose;
230 printf ("best (v_lose=%d; m_lose=%d): ", v_lose, m_lose);
231 mpz_out_str (stdout, 10, a); puts (""); fflush (stdout);
232 }
233 if (m_lose < m_best)
234 {
235 m_best = m_lose;
236 printf ("best (v_lose=%d; m_lose=%d): ", v_lose, m_lose);
237 mpz_out_str (stdout, 10, a); puts (""); fflush (stdout);
238 }
239 }
240
241 mpz_clear (m);
242 mpz_clear (a);
243 for (f = 0; f < DIMS-1; f++)
244 mpf_clear (v[f]);
245 mpf_clear (f_merit);
246 mpf_clear (low_merit);
247 mpf_clear (high_merit);
248
249 printf ("done.\n");
250 return 0;
251 }
252