/*
  Author:  Pate Williams (c) 1997

  Algorithm 7.4.7 (Periods of an Elliptic Curve
  over R). See "A Course in Computational Alge-
  braic Number Theory" by Henri Cohen page 398.
*/

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

double AGM(double a, double b)
/* calculates the Arithmetic-Geometric Mean of a and b */
{
  double a0 = a, a1, b0 = b, b1;

  do {
    a1 = 0.5 * (a0 + b0);
    b1 = sqrt(a0 * b0);
    a0 = a1, b0 = b1;
  } while (fabs(a1 - b1) > 1.0e-20);
  return a1;
}

void solve(double a0, double a1, double a2,
           struct complex *z1, struct complex *z2,
           struct complex *z3)
/* solves the cubic equation
   z ^ 3 + a2 * z ^ 2 + a1 * z + a0 = 0 */
{
  double a, d, q, r, t1, t2;
  struct complex expon, s, t, s1, s2, u, v;

  q = a1 / 3.0 - a2 * a2 / 9.0;
  r = (a1 * a2 - 3.0 * a0) / 6.0 - a2 * a2 * a2 / 27.0;
  d = sqrt(q * q * q + r * r);
  t1 = r + d;
  t2 = r - d;
  expon.x = 1.0 / 3.0;
  expon.y = 0.0;
  s.x = t1;
  s.y = 0.0;
  if (t1 != 0.0) s1 = cpow(s, expon);
  else s1.x = s1.y = 0.0;
  s.x = t2;
  if (t2 != 0.0 )s2 = cpow(s, expon);
  else s2.x = s2.y = 0.0;
  s = cadd(s1, s2);
  a = a2 / 3.0;
  t.x = - a;
  t.y = 0.0;
  *z1 = cadd(s, t);
  s.x *= - 0.5;
  s.y *= - 0.5;
  u = csub(s1, s2);
  v.x = 0.0;
  v.y = 0.5 * sqrt(3.0);
  *z2 = cadd(cadd(s, t), cmul(u, v));
  *z3 = cconj(*z2);
}

void periods(double a1, double a2, double a3,
             double a4, double a6,
             struct complex *omega1,
             struct complex *omega2)
{
  double a, b, c, t;
  double b2, b4, b6, b8, c0, c1, c2, delta, e[3], e1;
  int i, j;
  struct complex z1, z2, z3;

  b2 = a1 * a1 + 4.0 * a2;
  b4 = a1 * a3 + 2.0 * a4;
  b6 = a3 * a3 + 4.0 * a6;
  b8 = a1 * a1 * a6 + 4.0 * a2 * a6
     - a1 * a3 * a4 + a2 * a3 * a3
     - a4 * a4;
  delta = - b2 * b2 * b8 - 8.0 * b4 * b4 * b4
        - 27.0 * b6 * b6 + 9.0 * b2 * b4 * b6;
  c0 = b6 / 4.0;
  c1 = b4 / 2.0;
  c2 = b2 / 4.0;
  solve(c0, c1, c2, &z1, &z2, &z3);
  if (delta >= 0.0) {
    e[0] = z1.x;
    e[1] = z2.x;
    e[2] = z3.x;
    /* sort the real roots into descending order */
    for (i = 0; i < 2; i++)
      for (j = i + 1; j < 3; j++)
        if (e[i] < e[j])
          t = e[i], e[i] = e[j], e[j] = t;
    omega1->x = M_PI / AGM(sqrt(e[0] - e[2]), sqrt(e[0] - e[1]));
    omega1->y = 0.0;
    omega2->x = 0.0;
    omega2->y = M_PI / AGM(sqrt(e[0] - e[2]), sqrt(e[1] - e[2]));
  }
  else {
    e1 = z1.x;
    a = 3.0 * e1 + b2 / 4.0;
    b = sqrt(3.0 * e1 * e1 + 0.5 * b2 * e1 + 0.25 * b4);
    c = 2.0 * sqrt(b);
    omega1->x = 2.0 * M_PI / AGM(c, sqrt(2.0 * b + a));
    omega1->y = 0.0;
    omega2->x = - 0.5 * omega1->x;
    omega2->y = M_PI / AGM(c, sqrt(2.0 * b - a));
  }
}

int main(void)
{
  double a1 = 1.0, a2 = 2.0, a3 = 3.0;
  double a4 = 4.0, a6 = 6.0;
  double b0 = - 2.0, b1 = 0.0, b2 = 0.0;
  struct complex omega1, omega2;
  struct complex z1, z2, z3;

  solve(b0, b1, b2, &z1, &z2, &z3);
  printf("the roots of z ^ 3 - 2 = 0 are:\n");
  printf("z1 = "); cprint(z1);
  printf("z2 = "); cprint(z2);
  printf("z3 = "); cprint(z3);
  periods(a1, a2, a3, a4, a6, &omega1, &omega2);
  printf("omega1 = "); cprint(omega1);
  printf("omega2 = "); cprint(omega2);
  omega2 = cdiv(omega2, omega1);
  printf("2 / 1  = "); cprint(omega2);
  return 0;
}