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.2

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