/*
  Author:  Pate Williams (c) 1997

  Heuristic Algorithm 5.6.13 (Class Number for D > 0).
  See "A Course in Computational Algebraic Number
  Theory" by Henri Cohen pages 268 - 269.
*/

#include <assert.h>
#include <float.h>
#include <math.h>
#include <stdio.h>
#include "lip.h"

double ztoreal(verylong za)
{
  double x;
  char digit[8192];
  long i, length = z2log(za);
  verylong zb = 0, zc = 0;

  assert(length <= 8192l);
  zcopy(za, &zb);
  for (i = 0; i < length; i++) {
    zlowbits(zb, 1l, &zc);
    digit[i] = (char) zc[1];
    zrshift(zb, 1l, &zc);
    zcopy(zc, &zb);
  }
  x = digit[length - 1];
  for (i = length - 2; i >= 0; i--)
    x = 2 * x + digit[i];
  return x;
}

double regulator(verylong zD, verylong za,
                 verylong zb)
{
  double L, f = sqrt(ztoreal(zD)), R = 1.0;
  long e = 0, l;
  verylong zA = 0, zd = 0, zp = 0, zq = 0, zr = 0;
  verylong zs = 0, zt = 0, zx = 0, zy = 0;
  verylong za2 = 0, zb2 = 0, zp1 = 0, zq1 = 0;

  L = DBL_MAX / f;
  l = log(L) / log(2.0);
  zsqrt(zD, &zd, &zr);
  zcopy(zb, &zp);
  zlshift(za, 1l, &za2);
  zcopy(za2, &zq);
  zmod(zb, za2, &zb2);
  zsq(zp, &zr);
  zsub(zD, zr, &zs);
  zdiv(zs, zq, &zq1, &zr);
  L2:
  zadd(zp, zd, &zt);
  zdiv(zt, zq, &zA, &zr);
  if (zscompare(zr, 0l) < 0) {
    zcopy(zq, &zt);
    zabs(&zt);
    zadd(zr, zt, &zs);
    zcopy(zs, &zr);
  }
  zcopy(zp, &zp1);
  zsub(zd, zr, &zp);
  zcopy(zq, &zt);
  zsub(zp, zp1, &zs);
  zmul(zA, zs, &zx);
  zsub(zq1, zx, &zy);
  zcopy(zy, &zq);
  zcopy(zt, &zq1);
  R *= (ztoreal(zp) + f) / ztoreal(zq);
  if (R >= L) R /= L, e++;
  zmod(zp, za2, &zx);
  if (zcompare(zq, za2) != 0 || zcompare(zx, zb2) != 0)
    goto L2;
  R = log(R) + e * l * log(2.0);
  zfree(&zA);
  zfree(&zd);
  zfree(&zp);
  zfree(&zq);
  zfree(&zr);
  zfree(&zs);
  zfree(&zt);
  zfree(&zx);
  zfree(&zy);
  zfree(&za2);
  zfree(&zb2);
  zfree(&zp1);
  zfree(&zq1);
  return R;
}

int class_number(verylong zD, verylong za,
                 verylong zb, double *r)
{
  int j;
  double D, h1, product, R = regulator(zD, za, zb);
  long c, h, i, m, p;
  verylong zp = 0;

  zpstart2();
  D = ztoreal(zD);
  h = 0;
  if (R >= pow(D, 0.25)) {
    h1 = sqrt(ztoreal(zD)) / (2.0 * R);
    c = 0;
    L3:
    /* compute the Euler product */
    product = 1.0;
    for (i = 0; i < 500; i++) {
      p = zpnext();
      zintoz(p, &zp);
      j = zjacobi(zD, zp);
      if (j != 0) product /= 1.0 - j / (double) p;
    }
    h1 *= product;
    m = h1 + 0.5;
    if (fabs(m - h1) > 0.1) {
      c = 0;
      goto L3;
    }
    if (m != h) {
      h = m;
      c = 1;
      goto L3;
    }
    else {
      c++;
      if (c <= 5) goto L3;
    }
  }
  zfree(&zp);
  *r = R;
  return h;
}

int main(void)
{
  double R;
  long a[5] = {1, 2, 3, 4, 5};
  long b[5] = {4, 6, 9, 13, 15};
  long D[5] = {8, 12, 33, 89, 105}, i;
  verylong zD = 0, za = 0, zb = 0;

  printf("--------------\n");
  printf(" D  h   R(D)\n");
  printf("--------------\n");
  for (i = 0; i <  5; i++) {
    zintoz(D[i], &zD);
    zintoz(a[i], &za);
    zintoz(b[i], &zb);
    printf("%3ld ", D[i]);
    printf("%d ", class_number(zD, za, zb, &R));
    printf("%lf\n", R);
  }
  printf("--------------\n");
  zfree(&zD);
  zfree(&za);
  zfree(&zb);
  return 0;
}