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 }