/*
  Author:  Pate Williams (c) 1997

  Algorithm 7.4.2 (Reduction). See "A Course in
  Computational Algebraic Number Theory" by Henri
  Cohen page 395.
*/

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

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 main(void)
{
  long A[2][2];
  struct complex tau = {2.0, 0.5};

  printf("tau = "); cprint(tau);
  reduction(A, &tau);
  printf("%2ld %2ld\n", A[0][0], A[0][1]);
  printf("%2ld %2ld\n", A[1][0], A[1][1]);
  printf("tau = "); cprint(tau);
  return 0;
}