tests/rand/findlc.c

   1 /*
   2 Copyright 2000 Free Software Foundation, Inc.
   3
   4 This file is part of the GNU MP Library.
   5
   6 The GNU MP Library is free software; you can redistribute it and/or modify
   7 it under the terms of the GNU Lesser General Public License as published by
   8 the Free Software Foundation; either version 3 of the License, or (at your
   9 option) any later version.
  10
  11 The GNU MP Library is distributed in the hope that it will be useful, but
  12 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  13 or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
  14 License for more details.
  15
  16 You should have received a copy of the GNU Lesser General Public License
  17 along 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) \
  28 static /**/const char *const rcsid[] = { (char *)rcsid, "\100(#)" msg }
  29
  30 RCSID("$Id$");
  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 }