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