github.com/aergoio/aergo@v1.3.1/libtool/src/gmp-6.1.2/demos/qcn.c (about)

     1  /* Use mpz_kronecker_ui() to calculate an estimate for the quadratic
     2     class number h(d), for a given negative fundamental discriminant, using
     3     Dirichlet's analytic formula.
     4  
     5  Copyright 1999-2002 Free Software Foundation, Inc.
     6  
     7  This file is part of the GNU MP Library.
     8  
     9  This program is free software; you can redistribute it and/or modify it
    10  under the terms of the GNU General Public License as published by the Free
    11  Software Foundation; either version 3 of the License, or (at your option)
    12  any later version.
    13  
    14  This program is distributed in the hope that it will be useful, but WITHOUT
    15  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
    16  FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License for
    17  more details.
    18  
    19  You should have received a copy of the GNU General Public License along with
    20  this program.  If not, see https://www.gnu.org/licenses/.  */
    21  
    22  
    23  /* Usage: qcn [-p limit] <discriminant>...
    24  
    25     A fundamental discriminant means one of the form D or 4*D with D
    26     square-free.  Each argument is checked to see it's congruent to 0 or 1
    27     mod 4 (as all discriminants must be), and that it's negative, but there's
    28     no check on D being square-free.
    29  
    30     This program is a bit of a toy, there are better methods for calculating
    31     the class number and class group structure.
    32  
    33     Reference:
    34  
    35     Daniel Shanks, "Class Number, A Theory of Factorization, and Genera",
    36     Proc. Symp. Pure Math., vol 20, 1970, pages 415-440.
    37  
    38  */
    39  
    40  #include <math.h>
    41  #include <stdio.h>
    42  #include <stdlib.h>
    43  #include <string.h>
    44  
    45  #include "gmp.h"
    46  
    47  #ifndef M_PI
    48  #define M_PI  3.14159265358979323846
    49  #endif
    50  
    51  
    52  /* A simple but slow primality test.  */
    53  int
    54  prime_p (unsigned long n)
    55  {
    56    unsigned long  i, limit;
    57  
    58    if (n == 2)
    59      return 1;
    60    if (n < 2 || !(n&1))
    61      return 0;
    62  
    63    limit = (unsigned long) floor (sqrt ((double) n));
    64    for (i = 3; i <= limit; i+=2)
    65      if ((n % i) == 0)
    66        return 0;
    67  
    68    return 1;
    69  }
    70  
    71  
    72  /* The formula is as follows, with d < 0.
    73  
    74  	       w * sqrt(-d)      inf      p
    75  	h(d) = ------------ *  product --------
    76  		  2 * pi         p=2   p - (d/p)
    77  
    78  
    79     (d/p) is the Kronecker symbol and the product is over primes p.  w is 6
    80     when d=-3, 4 when d=-4, or 2 otherwise.
    81  
    82     Calculating the product up to p=infinity would take a long time, so for
    83     the estimate primes up to 132,000 are used.  Shanks found this giving an
    84     accuracy of about 1 part in 1000, in normal cases.  */
    85  
    86  unsigned long  p_limit = 132000;
    87  
    88  double
    89  qcn_estimate (mpz_t d)
    90  {
    91    double  h;
    92    unsigned long  p;
    93  
    94    /* p=2 */
    95    h = sqrt (-mpz_get_d (d)) / M_PI
    96      * 2.0 / (2.0 - mpz_kronecker_ui (d, 2));
    97  
    98    if (mpz_cmp_si (d, -3) == 0)       h *= 3;
    99    else if (mpz_cmp_si (d, -4) == 0)  h *= 2;
   100  
   101    for (p = 3; p <= p_limit; p += 2)
   102      if (prime_p (p))
   103        h *= (double) p / (double) (p - mpz_kronecker_ui (d, p));
   104  
   105    return h;
   106  }
   107  
   108  
   109  void
   110  qcn_str (char *num)
   111  {
   112    mpz_t  z;
   113  
   114    mpz_init_set_str (z, num, 0);
   115  
   116    if (mpz_sgn (z) >= 0)
   117      {
   118        mpz_out_str (stdout, 0, z);
   119        printf (" is not supported (negatives only)\n");
   120      }
   121    else if (mpz_fdiv_ui (z, 4) != 0 && mpz_fdiv_ui (z, 4) != 1)
   122      {
   123        mpz_out_str (stdout, 0, z);
   124        printf (" is not a discriminant (must == 0 or 1 mod 4)\n");
   125      }
   126    else
   127      {
   128        printf ("h(");
   129        mpz_out_str (stdout, 0, z);
   130        printf (") approx %.1f\n", qcn_estimate (z));
   131      }
   132    mpz_clear (z);
   133  }
   134  
   135  
   136  int
   137  main (int argc, char *argv[])
   138  {
   139    int  i;
   140    int  saw_number = 0;
   141  
   142    for (i = 1; i < argc; i++)
   143      {
   144        if (strcmp (argv[i], "-p") == 0)
   145  	{
   146  	  i++;
   147  	  if (i >= argc)
   148  	    {
   149  	      fprintf (stderr, "Missing argument to -p\n");
   150  	      exit (1);
   151  	    }
   152  	  p_limit = atoi (argv[i]);
   153  	}
   154        else
   155  	{
   156  	  qcn_str (argv[i]);
   157  	  saw_number = 1;
   158  	}
   159      }
   160  
   161    if (! saw_number)
   162      {
   163        /* some default output */
   164        qcn_str ("-85702502803");           /* is 16259   */
   165        qcn_str ("-328878692999");          /* is 1499699 */
   166        qcn_str ("-928185925902146563");    /* is 52739552 */
   167        qcn_str ("-84148631888752647283");  /* is 496652272 */
   168        return 0;
   169      }
   170  
   171    return 0;
   172  }