Return to qcn.c CVS log

Up to [local] / OpenXM_contrib / gmp / demos

Annotation of OpenXM_contrib/gmp/demos/qcn.c, Revision 1.1.1.1

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