/*
  Author:  Pate Williams (c) 1997

  Algorithm 7.4.5 (p(z) and p'(z)). 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 p(struct complex omega1, struct complex q,
                 struct complex u)
{
  struct complex a, qn, sum, sum1, sum2, sum3, term;

  qn = q;
  sum1.x = sum1.y = 0.0;
  do {
    a = cmul(qn, u);
    a.x = 1.0 - a.x;
    a.y = - a.y;
    a = cmul(a, a);
    term = cdiv(qn, a);
    sum1 = cadd(sum1, term);
    qn = cmul(q, qn);
  } while (cabs(term) > 1.0e-20);
  sum1 = cmul(sum1, u);
  qn = q;
  sum2.x = sum2.y = 0.0;
  do {
    a = csub(qn, u);
    term = cdiv(qn, cmul(a, a));
    sum2 = cadd(sum2, term);
    qn = cmul(q, qn);
  } while (cabs(term) > 1.0e-20);
  sum2 = cmul(sum2, u);
  qn = q;
  sum3.x = sum3.y = 0.0;
  do {
    a.x = 1.0 - qn.x;
    a.y = - qn.x;
    term = cdiv(qn, cmul(a, a));
    sum3 = cadd(sum3, term);
    qn = cmul(q, qn);
  } while (cabs(term) > 1.0e-20);
  sum3.x *= 2.0;
  sum3.y *= 2.0;
  a.x = 1.0 - u.x;
  a.y = - u.y;
  sum = cadd(cdiv(u, cmul(a, a)), sum1);
  sum = cadd(sum, sum2);
  sum = csub(sum, sum3);
  sum.x += 1.0 / 12.0;
  a.x = 0.0;
  a.y = 2.0 * M_PI;
  a = cdiv(a, omega1);
  return cmul(cmul(a, a), sum);
}

struct complex pp(struct complex omega1, struct complex q,
                  struct complex u)
{
  struct complex a, b, c, d, qn, sum, sum1, sum2, term;

  qn = q;
  sum1.x = sum1.y = 0.0;
  do {
    a = cmul(qn, u);
    b.x = 1.0 + a.x;
    b.y = a.y;
    c.x = 1.0 - a.x;
    c.y = - a.y;
    d = cmul(c, cmul(c, c));
    term = cdiv(b, d);
    term = cmul(qn, term);
    qn = cmul(q, qn);
  } while (cabs(term) > 1.0e-20);
  qn = q;
  sum2.x = sum2.y = 0.0;
  do {
    a = cadd(qn, u);
    b = csub(qn, u);
    d = cmul(b, cmul(b, b));
    term = cdiv(a, d);
    term = cmul(qn, term);
    qn = cmul(q, qn);
  } while (cabs(term) > 1.0e-20);
  a.x = 1.0 + u.x;
  a.y = u.y;
  b.x = 1.0 - u.x;
  b.y = - u.y;
  d = cmul(b, cmul(b, b));
  sum = cadd(cdiv(a, d), sum1);
  sum = cadd(sum, sum2);
  a.x = 0.0;
  a.y = 2.0 * M_PI;
  b = cdiv(a, omega1);
  c = cmul(b, cmul(b, b));
  d = cmul(c, u);
  return cmul(d, sum);
}

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;
  }
}

int pandpp(struct complex omega1, struct complex omega2,
           struct complex z, struct complex *P,
           struct complex *PP)
{
   long A[2][2], n;
   struct complex a, b, q, ratio, tau, temp, u;

   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);
   omega2.x *= A[1][0];
   omega2.y *= A[1][0];
   omega1.x *= A[1][1];
   omega1.y *= A[1][1];
   omega1 = cadd(omega2, omega1);
   z = cdiv(z, omega1);
   a.x = 0.0;
   a.y = 2.0 * M_PI;
   n = floor(0.5 + z.y / tau.y);
   b.x = n;
   b.y = 0.0;
   z = csub(z, b);
   b.x = floor(z.x + 0.5);
   b.y = 0.0;
   z = csub(z, b);
   if (cabs(z) == 0.0) return 0;
   q = cexp(cmul(a, tau));
   u = cexp(cmul(a, z));
   *P  =  p(omega1, q, u);
   *PP = pp(omega1, q, u);
   return 1;
}

int main(void)
{
  int flag;
  struct complex P, PP;
  struct complex omega1 = {1.0, 4.0};
  struct complex omega2 = {1.5, 2.0};
  struct complex z = {0.5, 0.5};

  flag = pandpp(omega1, omega2, z, &P, &PP);
  printf("omega1 = "); cprint(omega1);
  printf("omega2 = "); cprint(omega2);
  printf("z      = "); cprint(z);
  if (flag == 0)
    printf("z is in the lattice\n");
  else {
    printf("p(z)  = "); cprint(P);
    printf("p'(z) = "); cprint(PP);
  }
  return 0;
}