demos/qcn.c

   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, 2000, 2001, 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 http://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 }