/*
  Author:  Pate Williams (c) 1997

  Algorithm 7.4.3 (g2 and g3). See "A Course in
  Computational Algebraic Number Theory" by Henri
  Cohen page 396.
*/

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

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, "division 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);
}

struct complex g2(struct complex omega1, struct complex q)
{
  double rn;
  long n = 1;
  struct complex denom, expon, g2, qn, sum, term;

  sum.x = sum.y = 0.0;
  qn = q;
  do {
    rn = n;
    rn = rn * rn * rn;
    denom.x = 1.0 - qn.x;
    denom.y = - qn.y;
    term = cdiv(qn, denom);
    term.x *= rn;
    term.y *= rn;
    sum = cadd(sum, term);
    qn = cmul(q, qn);
    n++;
  } while (cabs(term) > 1.0e-20);
  sum.x *= 240.0;
  sum.y *= 240.0;
  sum.x += 1.0;
  expon.x = 2.0 * M_PI;
  expon.y = 0.0;
  term = cdiv(expon, omega1);
  expon.x = 4.0;
  expon.y = 0.0;
  term = cpow(term, expon);
  g2 = cmul(sum, term);
  g2.x /= 12.0;
  g2.y /= 12.0;
  return g2;
}

struct complex g3(struct complex omega1, struct complex q)
{
  double rn;
  long n = 1;
  struct complex denom, expon, g3, qn, sum, term;

  sum.x = sum.y = 0.0;
  qn = q;
  do {
    rn = n;
    rn = pow(rn, 5.0);
    denom.x = 1.0 - qn.x;
    denom.y = - qn.y;
    term = cdiv(qn, denom);
    term.x *= rn;
    term.y *= rn;
    sum = cadd(sum, term);
    qn = cmul(q, qn);
    n++;
  } while (cabs(term) > 1.0e-20);
  sum.x *= 504.0;
  sum.y *= 504.0;
  sum.x = 1.0 - sum.x;
  sum.y = - sum.y;
  expon.x = 2.0 * M_PI;
  expon.y = 0.0;
  term = cdiv(expon, omega1);
  expon.x = 6.0;
  expon.y = 0.0;
  term = cpow(term, expon);
  g3 = cmul(sum, term);
  g3.x /= 216.0;
  g3.y /= 216.0;
  return g3;
}

void reduction(long A[2][2], struct complex *tau)
{
  long a, b, n;
  struct complex ctau, m;

  A[0][0] = A[1][1] = 1, A[0][1] = A[1][0] = 0;
  L2:
  n = tau->x + 0.5;
  tau->x -= n;
  A[0][0] -= n * A[1][0];
  A[0][1] -= n * A[1][1];
  ctau = cconj(*tau);
  m = cmul(*tau, ctau);
  if (fabs(m.x) < 1.0) {
    ctau.x = - ctau.x;
    ctau.y = - ctau.y;
    *tau = cdiv(ctau, m);
    a = A[0][0], b = A[0][1];
    A[0][0] = - A[1][0], A[0][1] = - A[1][1];
    A[1][0] = a, A[1][1] = b;
    goto L2;
  }
}

void g2andg3(struct complex omega1, struct complex omega2,
             struct complex *G2, struct complex *G3)
{
   long A[2][2];
   struct complex omega, q, ratio, tau, temp;

   ratio = cdiv(omega2, omega1);
   if (ratio.y < 0.0)
     temp = omega1, omega1 = omega2, omega2 = temp;
   tau = cdiv(omega2, omega1);
   printf("tau  = "); cprint(tau);
   reduction(A, &tau);
   printf("tau' = "); cprint(tau);
   ratio.x = 0.0;
   ratio.y = 2.0 * M_PI;
   q = cexp(cmul(ratio, tau));
   omega2.x *= A[1][0];
   omega2.y *= A[1][0];
   omega1.x *= A[1][1];
   omega1.y *= A[1][1];
   omega = cadd(omega2, omega1);
   *G2 = g2(omega, q);
   *G3 = g3(omega, q);
}

int main(void)
{
  struct complex G2, G3;
  struct complex omega1 = {1.0, 4.0};
  struct complex omega2 = {1.5, 2.0};

  g2andg3(omega1, omega2, &G2, &G3);
  printf("omega1 = "); cprint(omega1);
  printf("omega2 = "); cprint(omega2);
  printf("g2 = "); cprint(G2);
  printf("g3 = "); cprint(G3);
  return 0;
}