1139826Simp/* 253541SshinCopyright 2000 Free Software Foundation, Inc. 353541Sshin 453541SshinThis file is part of the GNU MP Library test suite. 553541Sshin 653541SshinThe GNU MP Library test suite is free software; you can redistribute it 753541Sshinand/or modify it under the terms of the GNU General Public License as 853541Sshinpublished by the Free Software Foundation; either version 3 of the License, 953541Sshinor (at your option) any later version. 1053541Sshin 1153541SshinThe GNU MP Library test suite is distributed in the hope that it will be 1253541Sshinuseful, but WITHOUT ANY WARRANTY; without even the implied warranty of 1353541SshinMERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General 1453541SshinPublic License for more details. 1553541Sshin 1653541SshinYou should have received a copy of the GNU General Public License along with 1753541Sshinthe GNU MP Library test suite. If not, see https://www.gnu.org/licenses/. */ 1853541Sshin 1953541Sshin#include <stdio.h> 2053541Sshin#include <stdlib.h> 2153541Sshin#include <unistd.h> 2253541Sshin#include <signal.h> 2353541Sshin#include <math.h> 2453541Sshin#include "gmpstat.h" 2553541Sshin 2653541Sshin#define RCSID(msg) \ 2753541Sshinstatic /**/const char *const rcsid[] = { (char *)rcsid, "\100(#)" msg } 2853541Sshin 2953541SshinRCSID("$Id: findlc.c,v 1.1.1.4 2020/09/27 00:27:04 mrg Exp $"); 3053541Sshin 3153541Sshinint g_debug = 0; 32139826Simp 3353541Sshinstatic mpz_t a; 34171328Srwatson 35171328Srwatsonstatic void 3653541Sshinsh_status (int sig) 3753541Sshin{ 3853541Sshin printf ("sh_status: signal %d caught. dumping status.\n", sig); 3953541Sshin 4053541Sshin printf (" a = "); 4153541Sshin mpz_out_str (stdout, 10, a); 4253541Sshin printf ("\n"); 4353541Sshin fflush (stdout); 4453541Sshin 4553541Sshin if (SIGSEGV == sig) /* remove SEGV handler */ 4653541Sshin signal (SIGSEGV, SIG_DFL); 4753541Sshin} 4853541Sshin 4953541Sshin/* Input is a modulus (m). We shall find multiplier (a) and adder (c) 5053541Sshin conforming to the rules found in the first comment block in file 5153541Sshin mpz/urandom.c. 5253541Sshin 5353541Sshin Then run a spectral test on the generator and discard any 5453541Sshin multipliers not passing. */ 5553541Sshin 5653541Sshin/* TODO: 5753541Sshin 5853541Sshin . find a better algorithm than a+=8; bigger jumps perhaps? 5953541Sshin 6053541Sshin*/ 6153541Sshin 6253541Sshinvoid 6353541Sshinmpz_true_random (mpz_t s, unsigned long int nbits) 6453541Sshin{ 6562587Sitojun#if __FreeBSD__ 6653541Sshin FILE *fs; 6755205Speter char c[1]; 6854263Sshin int i; 6954263Sshin 70171328Srwatson mpz_set_ui (s, 0); 7153541Sshin for (i = 0; i < nbits; i += 8) 72171328Srwatson { 73171328Srwatson for (;;) 74171328Srwatson { 75171328Srwatson int nread; 7655205Speter fs = fopen ("/dev/random", "r"); 7753541Sshin nread = fread (c, 1, 1, fs); 7853541Sshin 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