qcn.c
上传用户:qaz666999
上传日期:2022-08-06
资源大小:2570k
文件大小:4k
- /* Use mpz_kronecker_ui() to calculate an estimate for the quadratic
- class number h(d), for a given negative fundamental discriminant, using
- Dirichlet's analytic formula.
- Copyright 1999, 2000, 2001, 2002 Free Software Foundation, Inc.
- This file is part of the GNU MP Library.
- This program is free software; you can redistribute it and/or modify it
- under the terms of the GNU General Public License as published by the Free
- Software Foundation; either version 3 of the License, or (at your option)
- any later version.
- This program is distributed in the hope that it will be useful, but WITHOUT
- ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
- FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for
- more details.
- You should have received a copy of the GNU General Public License along with
- this program. If not, see http://www.gnu.org/licenses/. */
- /* Usage: qcn [-p limit] <discriminant>...
- A fundamental discriminant means one of the form D or 4*D with D
- square-free. Each argument is checked to see it's congruent to 0 or 1
- mod 4 (as all discriminants must be), and that it's negative, but there's
- no check on D being square-free.
- This program is a bit of a toy, there are better methods for calculating
- the class number and class group structure.
- Reference:
- Daniel Shanks, "Class Number, A Theory of Factorization, and Genera",
- Proc. Symp. Pure Math., vol 20, 1970, pages 415-440.
- */
- #include <math.h>
- #include <stdio.h>
- #include <stdlib.h>
- #include <string.h>
- #include "gmp.h"
- #ifndef M_PI
- #define M_PI 3.14159265358979323846
- #endif
- /* A simple but slow primality test. */
- int
- prime_p (unsigned long n)
- {
- unsigned long i, limit;
- if (n == 2)
- return 1;
- if (n < 2 || !(n&1))
- return 0;
- limit = (unsigned long) floor (sqrt ((double) n));
- for (i = 3; i <= limit; i+=2)
- if ((n % i) == 0)
- return 0;
- return 1;
- }
- /* The formula is as follows, with d < 0.
- w * sqrt(-d) inf p
- h(d) = ------------ * product --------
- 2 * pi p=2 p - (d/p)
- (d/p) is the Kronecker symbol and the product is over primes p. w is 6
- when d=-3, 4 when d=-4, or 2 otherwise.
- Calculating the product up to p=infinity would take a long time, so for
- the estimate primes up to 132,000 are used. Shanks found this giving an
- accuracy of about 1 part in 1000, in normal cases. */
- unsigned long p_limit = 132000;
- double
- qcn_estimate (mpz_t d)
- {
- double h;
- unsigned long p;
- /* p=2 */
- h = sqrt (-mpz_get_d (d)) / M_PI
- * 2.0 / (2.0 - mpz_kronecker_ui (d, 2));
- if (mpz_cmp_si (d, -3) == 0) h *= 3;
- else if (mpz_cmp_si (d, -4) == 0) h *= 2;
- for (p = 3; p <= p_limit; p += 2)
- if (prime_p (p))
- h *= (double) p / (double) (p - mpz_kronecker_ui (d, p));
- return h;
- }
- void
- qcn_str (char *num)
- {
- mpz_t z;
- mpz_init_set_str (z, num, 0);
- if (mpz_sgn (z) >= 0)
- {
- mpz_out_str (stdout, 0, z);
- printf (" is not supported (negatives only)n");
- }
- else if (mpz_fdiv_ui (z, 4) != 0 && mpz_fdiv_ui (z, 4) != 1)
- {
- mpz_out_str (stdout, 0, z);
- printf (" is not a discriminant (must == 0 or 1 mod 4)n");
- }
- else
- {
- printf ("h(");
- mpz_out_str (stdout, 0, z);
- printf (") approx %.1fn", qcn_estimate (z));
- }
- mpz_clear (z);
- }
- int
- main (int argc, char *argv[])
- {
- int i;
- int saw_number = 0;
- for (i = 1; i < argc; i++)
- {
- if (strcmp (argv[i], "-p") == 0)
- {
- i++;
- if (i >= argc)
- {
- fprintf (stderr, "Missing argument to -pn");
- exit (1);
- }
- p_limit = atoi (argv[i]);
- }
- else
- {
- qcn_str (argv[i]);
- saw_number = 1;
- }
- }
- if (! saw_number)
- {
- /* some default output */
- qcn_str ("-85702502803"); /* is 16259 */
- qcn_str ("-328878692999"); /* is 1499699 */
- qcn_str ("-928185925902146563"); /* is 52739552 */
- qcn_str ("-84148631888752647283"); /* is 496652272 */
- return 0;
- }
- return 0;
- }