findlc.c
上传用户:qaz666999
上传日期:2022-08-06
资源大小:2570k
文件大小:6k
源码类别:

数学计算

开发平台:

Unix_Linux

  1. /*
  2. Copyright 2000 Free Software Foundation, Inc.
  3. This file is part of the GNU MP Library.
  4. The GNU MP Library is free software; you can redistribute it and/or modify
  5. it under the terms of the GNU Lesser General Public License as published by
  6. the Free Software Foundation; either version 3 of the License, or (at your
  7. option) any later version.
  8. The GNU MP Library is distributed in the hope that it will be useful, but
  9. WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  10. or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
  11. License for more details.
  12. You should have received a copy of the GNU Lesser General Public License
  13. along with the GNU MP Library.  If not, see http://www.gnu.org/licenses/.  */
  14. #include <stdio.h>
  15. #include <stdlib.h>
  16. #include <unistd.h>
  17. #include <signal.h>
  18. #include <math.h>
  19. #include "gmp.h"
  20. #include "gmpstat.h"
  21. #define RCSID(msg) 
  22. static /**/const char *const rcsid[] = { (char *)rcsid, "100(#)" msg }
  23. RCSID("$Id$");
  24. int g_debug = 0;
  25. static mpz_t a;
  26. static void
  27. sh_status (int sig)
  28. {
  29.   printf ("sh_status: signal %d caught. dumping status.n", sig);
  30.   printf ("  a = ");
  31.   mpz_out_str (stdout, 10, a);
  32.   printf ("n");
  33.   fflush (stdout);
  34.   if (SIGSEGV == sig) /* remove SEGV handler */
  35.     signal (SIGSEGV, SIG_DFL);
  36. }
  37. /* Input is a modulus (m).  We shall find multiplier (a) and adder (c)
  38.    conforming to the rules found in the first comment block in file
  39.    mpz/urandom.c.
  40.    Then run a spectral test on the generator and discard any
  41.    multipliers not passing.  */
  42. /* TODO:
  43.    . find a better algorithm than a+=8; bigger jumps perhaps?
  44. */
  45. void
  46. mpz_true_random (mpz_t s, unsigned long int nbits)
  47. {
  48. #if __FreeBSD__
  49.   FILE *fs;
  50.   char c[1];
  51.   int i;
  52.   mpz_set_ui (s, 0);
  53.   for (i = 0; i < nbits; i += 8)
  54.     {
  55.       for (;;)
  56. {
  57.   int nread;
  58.   fs = fopen ("/dev/random", "r");
  59.   nread = fread (c, 1, 1, fs);
  60.   fclose (fs);
  61.   if (nread != 0)
  62.     break;
  63.   sleep (1);
  64. }
  65.       mpz_mul_2exp (s, s, 8);
  66.       mpz_add_ui (s, s, ((unsigned long int) c[0]) & 0xff);
  67.       printf ("%d random bitsn", i + 8);
  68.     }
  69.   if (nbits % 8 != 0)
  70.     mpz_mod_2exp (s, s, nbits);
  71. #endif
  72. }
  73. int
  74. main (int argc, char *argv[])
  75. {
  76.   const char usage[] = "usage: findlc [-dv] m2exp [low_merit [high_merit]]n";
  77.   int f;
  78.   int v_lose, m_lose, v_best, m_best;
  79.   int c;
  80.   int debug = 1;
  81.   int cnt_high_merit;
  82.   mpz_t m;
  83.   unsigned long int m2exp;
  84. #define DIMS 6 /* dimensions run in spectral test */
  85.   mpf_t v[DIMS-1]; /* spectral test result (there's no v
  86.    for 1st dimension */
  87.   mpf_t f_merit, low_merit, high_merit;
  88.   mpz_t acc, minus8;
  89.   mpz_t min, max;
  90.   mpz_t s;
  91.   mpz_init (m);
  92.   mpz_init (a);
  93.   for (f = 0; f < DIMS-1; f++)
  94.     mpf_init (v[f]);
  95.   mpf_init (f_merit);
  96.   mpf_init_set_d (low_merit, .1);
  97.   mpf_init_set_d (high_merit, .1);
  98.   while ((c = getopt (argc, argv, "a:di:hv")) != -1)
  99.     switch (c)
  100.       {
  101.       case 'd': /* debug */
  102. g_debug++;
  103. break;
  104.       case 'v': /* print version */
  105. puts (rcsid[1]);
  106. exit (0);
  107.       case 'h':
  108.       case '?':
  109.       default:
  110. fputs (usage, stderr);
  111. exit (1);
  112.       }
  113.   argc -= optind;
  114.   argv += optind;
  115.   if (argc < 1)
  116.     {
  117.       fputs (usage, stderr);
  118.       exit (1);
  119.     }
  120.   /* Install signal handler. */
  121.   if (SIG_ERR == signal (SIGSEGV, sh_status))
  122.     {
  123.       perror ("signal (SIGSEGV)");
  124.       exit (1);
  125.     }
  126.   if (SIG_ERR == signal (SIGHUP, sh_status))
  127.     {
  128.       perror ("signal (SIGHUP)");
  129.       exit (1);
  130.     }
  131.   printf ("findlc: version: %sn", rcsid[1]);
  132.   m2exp = atol (argv[0]);
  133.   mpz_init_set_ui (m, 1);
  134.   mpz_mul_2exp (m, m, m2exp);
  135.   printf ("m = 0x");
  136.   mpz_out_str (stdout, 16, m);
  137.   puts ("");
  138.   if (argc > 1) /* have low_merit */
  139.     mpf_set_str (low_merit, argv[1], 0);
  140.   if (argc > 2) /* have high_merit */
  141.     mpf_set_str (high_merit, argv[2], 0);
  142.   if (debug)
  143.     {
  144.       fprintf (stderr, "low_merit = ");
  145.       mpf_out_str (stderr, 10, 2, low_merit);
  146.       fprintf (stderr, "; high_merit = ");
  147.       mpf_out_str (stderr, 10, 2, high_merit);
  148.       fputs ("n", stderr);
  149.     }
  150.   mpz_init (minus8);
  151.   mpz_set_si (minus8, -8L);
  152.   mpz_init_set_ui (acc, 0);
  153.   mpz_init (s);
  154.   mpz_init_set_d (min, 0.01 * pow (2.0, (double) m2exp));
  155.   mpz_init_set_d (max, 0.99 * pow (2.0, (double) m2exp));
  156.   mpz_true_random (s, m2exp); /* Start.  */
  157.   mpz_setbit (s, 0); /* Make it odd.  */
  158.   v_best = m_best = 2*(DIMS-1);
  159.   for (;;)
  160.     {
  161.       mpz_add (acc, acc, s);
  162.       mpz_mod_2exp (acc, acc, m2exp);
  163. #if later
  164.       mpz_and_si (a, acc, -8L);
  165. #else
  166.       mpz_and (a, acc, minus8);
  167. #endif
  168.       mpz_add_ui (a, a, 5);
  169.       if (mpz_cmp (a, min) <= 0 || mpz_cmp (a, max) >= 0)
  170. continue;
  171.       spectral_test (v, DIMS, a, m);
  172.       for (f = 0, v_lose = m_lose = 0, cnt_high_merit = DIMS-1;
  173.    f < DIMS-1; f++)
  174. {
  175.   merit (f_merit, f + 2, v[f], m);
  176.   if (mpf_cmp_ui (v[f], 1 << (30 / (f + 2) + (f == 2))) < 0)
  177.     v_lose++;
  178.   if (mpf_cmp (f_merit, low_merit) < 0)
  179.     m_lose++;
  180.   if (mpf_cmp (f_merit, high_merit) >= 0)
  181.     cnt_high_merit--;
  182. }
  183.       if (0 == v_lose && 0 == m_lose)
  184. {
  185.   mpz_out_str (stdout, 10, a); puts (""); fflush (stdout);
  186.   if (0 == cnt_high_merit)
  187.     break; /* leave loop */
  188. }
  189.       if (v_lose < v_best)
  190. {
  191.   v_best = v_lose;
  192.   printf ("best (v_lose=%d; m_lose=%d): ", v_lose, m_lose);
  193.   mpz_out_str (stdout, 10, a); puts (""); fflush (stdout);
  194. }
  195.       if (m_lose < m_best)
  196. {
  197.   m_best = m_lose;
  198.   printf ("best (v_lose=%d; m_lose=%d): ", v_lose, m_lose);
  199.   mpz_out_str (stdout, 10, a); puts (""); fflush (stdout);
  200. }
  201.     }
  202.   mpz_clear (m);
  203.   mpz_clear (a);
  204.   for (f = 0; f < DIMS-1; f++)
  205.     mpf_clear (v[f]);
  206.   mpf_clear (f_merit);
  207.   mpf_clear (low_merit);
  208.   mpf_clear (high_merit);
  209.   printf ("done.n");
  210.   return 0;
  211. }