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     ! 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: }