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