/*
  Author:  Pate Williams (c) 1997

  Complex arithmetic. See "Numerical Computation
  Using C" by Robert Glassey Appendix I pages 263
  through 266.
  "Algorithm 7.6.1 (Hilbert Class Polynomial).
  Given a negative discriminant D, this algorithm
  computes the monic polynomial of degree h(D) in
  Z[X] of which j((D + sqrt(D)) / 2) is a root."
  -Henri Cohen- See "A Course in Computational
  Algebraic Number Theory" by Henri Cohen page
  415.
*/

#include <math.h>
#include <stdio.h>
#include <stdlib.h>

#define POLY_SIZE 8192l

struct complex cadd(struct complex u, struct complex v)
{
  struct complex w;

  w.x = u.x + v.x;
  w.y = u.y + v.y;
  return w;
}

struct complex csub(struct complex u, struct complex v)
{
  struct complex w;

  w.x = u.x - v.x;
  w.y = u.y - v.y;
  return w;
}

struct complex cdiv(struct complex u, struct complex v)
{
  double temp1, temp2;
  struct complex w;

  if (cabs(v) < 5.0e-20 * cabs(u)) {
    fprintf(stderr, "fatal error\ndivision by zero in cdiv");
    exit(1);
  }
  if (fabs(v.x) <= fabs(v.y)) {
    temp1 = v.x / v.y;
    temp2 = v.y + temp1 * v.x;
    w.x = (temp1 * u.x + u.y) / temp2;
    w.y = (temp1 * u.y - u.x) / temp2;
  }
  else {
    temp1 = v.y / v.x;
    temp2 = v.x + temp1 * v.y;
    w.x = (u.x + temp1 * u.y) / temp2;
    w.y = (u.y - temp1 * u.x) / temp2;
  }
  return w;
}

struct complex cmul(struct complex u, struct complex v)
{
  struct complex w;

  w.x = u.x * v.x - u.y * v.y;
  w.y = u.x * v.y + u.y * v.x;
  return w;
}

struct complex cconj(struct complex u)
{
  struct complex w;

  w.x = u.x;
  w.y = - u.y;
  return w;
}

struct complex csqrt(struct complex u)
{
  double r = sqrt(cabs(u));
  double theta = atan2(u.y, u.x), phi = 0.5 * theta;
  struct complex w;

  w.x = r * cos(phi);
  w.y = r * sin(phi);
  return w;
}

struct complex cexp(struct complex u)
{
  double e = exp(u.x);
  struct complex w;

  w.x = e * cos(u.y);
  w.y = e * sin(u.y);
  return w;
}

struct complex csin(struct complex u)
{
  struct complex w;

  w.x = sin(u.x) * cosh(u.y);
  w.y = cos(u.x) * sinh(u.y);
  return w;
}

struct complex ccos(struct complex u)
{
  struct complex w;

  w.x = cos(u.x) * cosh(u.y);
  w.y = - sin(u.x) * sinh(u.y);
  return w;
}

struct complex clog(struct complex u)
{
  double r = cabs(u);
  struct complex w;

  w.x = log(r);
  w.y = atan2(u.y, u.x);
  return w;
}

struct complex cpow(struct complex u, struct complex v)
{
  return cexp(cmul(v, clog(u)));
}

void cprint(struct complex u)
{
  printf("(%lf, %lf)\n", u.x, u.y);
}

long round(double x)
{
  if (x >= 0.0) return x + 0.5;
  return x - 0.5;
}

void cpoly_mul(long m, long n, struct complex *a,
               struct complex *b, struct complex *c,
               long *p)
{
  long i, j, k;
  struct complex ai, bj, sum, term;

  *p = m + n;
  for (k = 0; k <= *p; k++) {
    sum.x = sum.y = 0.0;
    for (i = 0; i <= k; i++) {
      j = k - i;
      if (i > m) ai.x = ai.y = 0.0;
      else ai = a[i];
      if (j > n) bj.x = bj.y = 0.0;
      else bj = b[j];
      term = cmul(ai, bj);
      sum = cadd(sum, term);
    }
    c[k] = sum;
  }
}

struct complex Delta(struct complex q)
{
  double sign = - 1.0;
  long e1, e2, n = 1;
  struct complex expon1 = {0.0, 0.0}, expon2 = {0.0, 0.0};
  struct complex sum = {0.0, 0.0}, term1, term2;

  do {
    e1 = n * (3 * n - 1) / 2;
    e2 = n * (3 * n + 1) / 2;
    expon1.x = e1;
    expon2.x = e2;
    term1 = cpow(q, expon1);
    term2 = cpow(q, expon2);
    term1 = cadd(term1, term2);
    term1.x *= sign;
    term1.y *= sign;
    sum = cadd(sum, term1);
    sign = - sign;
    n++;
  } while (cabs(term1) > 1.0e-20);
  sum.x += 1.0;
  expon1.x = 24;
  return cmul(q, cpow(sum, expon1));
}

struct complex j(struct complex tau)
{
  struct complex c = {0.0, 2.0 * M_PI}, f, f1, q, q2;

  q = cexp(cmul(c, tau));
  q2 = cmul(q, q);
  f = cdiv(Delta(q2), Delta(q));
  c.x = 256.0;
  c.y = 0.0;
  f1 = cmul(c, f);
  f1.x += 1.0;
  c.x = 3.0;
  f1 = cpow(f1, c);
  return cdiv(f1, f);
}

void Hilbert(long D, long *P, long *dP)
{
  long B, a, b, dQ, dR, i, t;
  double x;
  struct complex J, ca = {0.0, 0.0}, cb = {0.0, 0.0};
  struct complex cD = {0.0, 0.0}, sqrtD, tau;
  struct complex cP[POLY_SIZE], cQ[2], cR[POLY_SIZE];

  cD.x = D;
  sqrtD = csqrt(cD);
  cP[0].x = 1.0;
  cP[0].y = 0.0;
  *dP = 0;
  b = D % 2;
  if (b < 0) b += 2;
  B = sqrt(labs(D) / 3.0);
  L2:
  t = (b * b - D) / 4;
  a = max(b, 1);
  L3:
  if (t % a != 0) goto L4;
  ca.x = 2 * a;
  ca.y = 0.0;
  cb.x = - b;
  cb.y = 0.0;
  tau = cdiv(cadd(cb, sqrtD), ca);
  J = j(tau);
  if (a == b || a * a == t || b == 0) {
    cQ[0].x = - J.x;
    cQ[0].y = - J.y;
    cQ[1].x = 1.0;
    cQ[1].y = 0.0;
    dQ = 1;
  }
  else {
    x = cabs(J);
    cQ[0].x = x * x;
    cQ[0].y = 0.0;
    cQ[1].x = - 2.0 * J.x;
    cQ[1].y = 0.0;
    cQ[2].x = 1.0;
    cQ[2].y = 0.0;
    dQ = 2;
  }
  cpoly_mul(*dP, dQ, cP, cQ, cR, &dR);
  *dP = dR;
  for (i = 0; i <= dR; i++) cP[i] = cR[i];
  L4:
  a++;
  if (a * a <= t) goto L3;
  b += 2;
  if (b <= B) goto L2;
  for (i = 0; i <= *dP; i++) P[i] = round(cP[i].x);
}

int main(void)
{
  long D[17] = {- 4, - 7, - 8, - 11,
                - 15, - 20, - 24, -35, - 40, - 23,
                - 32, - 31, - 39, - 55, - 56, - 47,
                - 79};
  long dP, i, j, P[POLY_SIZE];

  printf(" D  h(D) Hilbert Class Polynomial\n");
  for (i = 0; i < 17; i++) {
    Hilbert(D[i], P, &dP);
    printf("%3ld %ld ", D[i], dP);
    for (j = dP; j >= 0; j--) printf("%11ld ", P[j]);
    printf("\n");
  }
  return 0;
}