Annotation of OpenXM_contrib/gmp/tests/rand/findlc.c, Revision 1.1.1.2
1.1 maekawa 1: /*
1.1.1.2 ! ohara 2: Copyright 2000 Free Software Foundation, Inc.
1.1 maekawa 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 2.1 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; see the file COPYING.LIB. If not, write to
18: the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
19: MA 02111-1307, USA.
20: */
21:
22: #include <stdio.h>
23: #include <stdlib.h>
24: #include <unistd.h>
25: #include <signal.h>
26: #include <math.h>
27: #include "gmp.h"
28: #include "gmpstat.h"
29:
30: #define RCSID(msg) \
31: static /**/const char *const rcsid[] = { (char *)rcsid, "\100(#)" msg }
32:
1.1.1.2 ! ohara 33: RCSID("$Id: findlc.c,v 1.4 2000/10/14 23:18:06 kevin Exp $");
1.1 maekawa 34:
35: int g_debug = 0;
36:
37: static mpz_t a;
38:
39: static void
40: sh_status (int sig)
41: {
42: printf ("sh_status: signal %d caught. dumping status.\n", sig);
43:
44: printf (" a = ");
45: mpz_out_str (stdout, 10, a);
46: printf ("\n");
47: fflush (stdout);
48:
49: if (SIGSEGV == sig) /* remove SEGV handler */
50: signal (SIGSEGV, SIG_DFL);
51: }
52:
53: /* Input is a modulus (m). We shall find multiplyer (a) and adder (c)
54: conforming to the rules found in the first comment block in file
55: mpz/urandom.c.
56:
57: Then run a spectral test on the generator and discard any
58: multipliers not passing. */
59:
60: /* TODO:
61:
62: . find a better algorithm than a+=8; bigger jumps perhaps?
63:
64: */
65:
66: void
67: mpz_true_random (mpz_t s, unsigned long int nbits)
68: {
69: #if __FreeBSD__
70: FILE *fs;
71: char c[1];
72: int i;
73:
74: mpz_set_ui (s, 0);
75: for (i = 0; i < nbits; i += 8)
76: {
77: for (;;)
78: {
79: int nread;
80: fs = fopen ("/dev/random", "r");
81: nread = fread (c, 1, 1, fs);
82: fclose (fs);
83: if (nread != 0)
84: break;
85: sleep (1);
86: }
87: mpz_mul_2exp (s, s, 8);
88: mpz_add_ui (s, s, ((unsigned long int) c[0]) & 0xff);
89: printf ("%d random bits\n", i + 8);
90: }
91: if (nbits % 8 != 0)
92: mpz_mod_2exp (s, s, nbits);
93: #endif
94: }
95:
96: int
97: main (int argc, char *argv[])
98: {
99: const char usage[] = "usage: findlc [-dv] m2exp [low_merit [high_merit]]\n";
100: int f;
101: int v_lose, m_lose, v_best, m_best;
102: int c;
103: int debug = 1;
104: int cnt_high_merit;
105: mpz_t m;
106: unsigned long int m2exp;
107: #define DIMS 6 /* dimensions run in spectral test */
108: mpf_t v[DIMS-1]; /* spectral test result (there's no v
109: for 1st dimension */
110: mpf_t f_merit, low_merit, high_merit;
111: mpz_t acc, minus8;
112: mpz_t min, max;
113: mpz_t s;
114:
115:
116: mpz_init (m);
117: mpz_init (a);
118: for (f = 0; f < DIMS-1; f++)
119: mpf_init (v[f]);
120: mpf_init (f_merit);
121: mpf_init_set_d (low_merit, .1);
122: mpf_init_set_d (high_merit, .1);
123:
124: while ((c = getopt (argc, argv, "a:di:hv")) != -1)
125: switch (c)
126: {
127: case 'd': /* debug */
128: g_debug++;
129: break;
130:
131: case 'v': /* print version */
132: puts (rcsid[1]);
133: exit (0);
134:
135: case 'h':
136: case '?':
137: default:
138: fputs (usage, stderr);
139: exit (1);
140: }
141:
142: argc -= optind;
143: argv += optind;
144:
145: if (argc < 1)
146: {
147: fputs (usage, stderr);
148: exit (1);
149: }
150:
151: /* Install signal handler. */
152: if (SIG_ERR == signal (SIGSEGV, sh_status))
153: {
154: perror ("signal (SIGSEGV)");
155: exit (1);
156: }
157: if (SIG_ERR == signal (SIGHUP, sh_status))
158: {
159: perror ("signal (SIGHUP)");
160: exit (1);
161: }
162:
163: printf ("findlc: version: %s\n", rcsid[1]);
164: m2exp = atol (argv[0]);
165: mpz_init_set_ui (m, 1);
166: mpz_mul_2exp (m, m, m2exp);
167: printf ("m = 0x");
168: mpz_out_str (stdout, 16, m);
169: puts ("");
170:
171: if (argc > 1) /* have low_merit */
172: mpf_set_str (low_merit, argv[1], 0);
173: if (argc > 2) /* have high_merit */
174: mpf_set_str (high_merit, argv[2], 0);
175:
176: if (debug)
177: {
178: fprintf (stderr, "low_merit = ");
179: mpf_out_str (stderr, 10, 2, low_merit);
180: fprintf (stderr, "; high_merit = ");
181: mpf_out_str (stderr, 10, 2, high_merit);
182: fputs ("\n", stderr);
183: }
184:
185: mpz_init (minus8);
186: mpz_set_si (minus8, -8L);
187: mpz_init_set_ui (acc, 0);
188: mpz_init (s);
189: mpz_init_set_d (min, 0.01 * pow (2.0, (double) m2exp));
190: mpz_init_set_d (max, 0.99 * pow (2.0, (double) m2exp));
191:
192: mpz_true_random (s, m2exp); /* Start. */
193: mpz_setbit (s, 0); /* Make it odd. */
194:
195: v_best = m_best = 2*(DIMS-1);
196: for (;;)
197: {
198: mpz_add (acc, acc, s);
199: mpz_mod_2exp (acc, acc, m2exp);
200: #if later
201: mpz_and_si (a, acc, -8L);
202: #else
203: mpz_and (a, acc, minus8);
204: #endif
205: mpz_add_ui (a, a, 5);
206: if (mpz_cmp (a, min) <= 0 || mpz_cmp (a, max) >= 0)
207: continue;
208:
209: spectral_test (v, DIMS, a, m);
210: for (f = 0, v_lose = m_lose = 0, cnt_high_merit = DIMS-1;
211: f < DIMS-1; f++)
212: {
213: merit (f_merit, f + 2, v[f], m);
214:
215: if (mpf_cmp_ui (v[f], 1 << (30 / (f + 2) + (f == 2))) < 0)
216: v_lose++;
217:
218: if (mpf_cmp (f_merit, low_merit) < 0)
219: m_lose++;
220:
221: if (mpf_cmp (f_merit, high_merit) >= 0)
222: cnt_high_merit--;
223: }
224:
225: if (0 == v_lose && 0 == m_lose)
226: {
227: mpz_out_str (stdout, 10, a); puts (""); fflush (stdout);
228: if (0 == cnt_high_merit)
229: break; /* leave loop */
230: }
231: if (v_lose < v_best)
232: {
233: v_best = v_lose;
234: printf ("best (v_lose=%d; m_lose=%d): ", v_lose, m_lose);
235: mpz_out_str (stdout, 10, a); puts (""); fflush (stdout);
236: }
237: if (m_lose < m_best)
238: {
239: m_best = m_lose;
240: printf ("best (v_lose=%d; m_lose=%d): ", v_lose, m_lose);
241: mpz_out_str (stdout, 10, a); puts (""); fflush (stdout);
242: }
243: }
244:
245: mpz_clear (m);
246: mpz_clear (a);
247: for (f = 0; f < DIMS-1; f++)
248: mpf_clear (v[f]);
249: mpf_clear (f_merit);
250: mpf_clear (low_merit);
251: mpf_clear (high_merit);
252:
253: printf ("done.\n");
254: return 0;
255: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to findlc.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / gmp / tests / rand