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

数学计算

开发平台:

Unix_Linux

  1. /* spect.c -- the spectral test */
  2. /*
  3. Copyright 1999 Free Software Foundation, Inc.
  4. This file is part of the GNU MP Library.
  5. The GNU MP Library is free software; you can redistribute it and/or modify
  6. it under the terms of the GNU Lesser General Public License as published by
  7. the Free Software Foundation; either version 3 of the License, or (at your
  8. option) any later version.
  9. The GNU MP Library is distributed in the hope that it will be useful, but
  10. WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  11. or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
  12. License for more details.
  13. You should have received a copy of the GNU Lesser General Public License
  14. along with the GNU MP Library.  If not, see http://www.gnu.org/licenses/.  */
  15. /* T is upper dimension.  Z_A is the LC multiplier, which is
  16.    relatively prime to Z_M, the LC modulus.  The result is put in
  17.    rop[] with v[t] in rop[t-2]. */
  18. /* BUGS: Due to lazy allocation scheme, maximum T is hard coded to MAXT. */
  19. #include <stdio.h>
  20. #include <stdlib.h>
  21. #include <unistd.h>
  22. #include <math.h>
  23. #include "gmp.h"
  24. #include "gmpstat.h"
  25. int g_debug = 0;
  26. int
  27. main (int argc, char *argv[])
  28. {
  29.   const char usage[] = "usage: spect [-d] a m nn";
  30.   int c;
  31.   unsigned int n;
  32.   mpz_t a, m;
  33.   mpf_t res[GMP_SPECT_MAXT], res_min[GMP_SPECT_MAXT], f_tmp;
  34.   register int f;
  35.   mpz_init (a);
  36.   mpz_init (m);
  37.   for (f = 0; f < GMP_SPECT_MAXT; f++)
  38.     {
  39.       mpf_init (res[f]);
  40.       mpf_init (res_min[f]);
  41.     }
  42.   mpf_init (f_tmp);
  43.   mpf_set_ui (res_min[0], 32768); /* 2^15 */
  44.   mpf_set_ui (res_min[1], 1024); /* 2^10 */
  45.   mpf_set_ui (res_min[2], 256); /* 2^8 */
  46.   mpf_set_ui (res_min[3], 64); /* 2^6 */
  47.   mpf_set_ui (res_min[4], 32); /* 2^5 */
  48.   while ((c = getopt (argc, argv, "dh")) != -1)
  49.     switch (c)
  50.       {
  51.       case 'd': /* debug */
  52. g_debug++;
  53. break;
  54.       case 'h':
  55.       default:
  56. fputs (usage, stderr);
  57. exit (1);
  58.       }
  59.   argc -= optind;
  60.   argv += optind;
  61.   if (argc < 3)
  62.     {
  63.       fputs (usage, stderr);
  64.       exit (1);
  65.     }
  66.   mpz_set_str (a, argv[0], 0);
  67.   mpz_set_str (m, argv[1], 0);
  68.   n = (unsigned int) atoi (argv[2]);
  69.   if (n + 1 > GMP_SPECT_MAXT)
  70.     n = GMP_SPECT_MAXT + 1;
  71.   spectral_test (res, n, a, m);
  72.   for (f = 0; f < n - 1; f++)
  73.     {
  74.       /* print v */
  75.       printf ("%d: v = ", f + 2);
  76.       mpf_out_str (stdout, 10, 4, res[f]);
  77. #ifdef PRINT_RAISED_BY_TWO_AS_WELL
  78.       printf (" (^2 = ");
  79.       mpf_mul (f_tmp, res[f], res[f]);
  80.       mpf_out_str (stdout, 10, 4, f_tmp);
  81.       printf (")");
  82. #endif /* PRINT_RAISED_BY_TWO_AS_WELL */
  83.       /* print merit */
  84.       printf (" m = ");
  85.       merit (f_tmp, f + 2, res[f], m);
  86.       mpf_out_str (stdout, 10, 4, f_tmp);
  87.       if (mpf_cmp (res[f], res_min[f]) < 0)
  88. printf ("t*** v too low ***");
  89.       if (mpf_get_d (f_tmp) < .1)
  90. printf ("t*** merit too low ***");
  91.       puts ("");
  92.     }
  93.   mpz_clear (a);
  94.   mpz_clear (m);
  95.   for (f = 0; f < GMP_SPECT_MAXT; f++)
  96.     {
  97.       mpf_clear (res[f]);
  98.       mpf_clear (res_min[f]);
  99.     }
  100.   mpf_clear (f_tmp);
  101.   return 0;
  102. }
  103. void
  104. debug_foo()
  105. {
  106.   if (0)
  107.     {
  108.       mpz_dump (0);
  109.       mpf_dump (0);
  110.     }
  111. }