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
5Copyright 1999, 2000, 2001, 2002 Free Software Foundation, Inc.
6
7This file is part of the GNU MP Library.
8
9This program is free software; you can redistribute it and/or modify it
10under the terms of the GNU General Public License as published by the Free
11Software Foundation; either version 3 of the License, or (at your option)
12any later version.
13
14This program is distributed in the hope that it will be useful, but WITHOUT
15ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
16FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License for
17more details.
18
19You should have received a copy of the GNU General Public License along with
20this 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.  */
53int
54prime_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
86unsigned long  p_limit = 132000;
87
88double
89qcn_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
109void
110qcn_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
136int
137main (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}
173