1/*
2Copyright 2000 Free Software Foundation, Inc.
3
4This file is part of the GNU MP Library test suite.
5
6The GNU MP Library test suite is free software; you can redistribute it
7and/or modify it under the terms of the GNU General Public License as
8published by the Free Software Foundation; either version 3 of the License,
9or (at your option) any later version.
10
11The GNU MP Library test suite is distributed in the hope that it will be
12useful, but WITHOUT ANY WARRANTY; without even the implied warranty of
13MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General
14Public License for more details.
15
16You should have received a copy of the GNU General Public License along with
17the 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) \
27static /**/const char *const rcsid[] = { (char *)rcsid, "\100(#)" msg }
28
29RCSID("$Id: findlc.c,v 1.1.1.4 2020/09/27 00:27:04 mrg Exp $");
30
31int g_debug = 0;
32
33static mpz_t a;
34
35static void
36sh_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
62void
63mpz_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
92int
93main (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