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 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) \ 28 static /**/const char *const rcsid[] = { (char *)rcsid, "\100(#)" msg } 29 30 RCSID("$Id: findlc.c,v 1.1.1.2 2013/11/29 07:49:48 mrg Exp $"); 31 32 int g_debug = 0; 33 34 static mpz_t a; 35 36 static void 37 sh_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 63 void 64 mpz_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 93 int 94 main (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