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