/*
  Author:  Pate Williams (c) 1997

  Algorithm 5.7.1 (Fundamental Unit Using
  Continued Fractions). See "A Course in
  Computational Algebraic Number Theory"
  by Henri Cohen page 270.
*/

#include <stdio.h>
#include "lip.h"

void fundamental(verylong za, verylong zb,
                 verylong zc, verylong *ze,
                 long *iterations)
{
  verylong zA = 0, zD = 0, zQ = 0, zd = 0, zf = 0;
  verylong zg = 0, zm = 0, zn = 0, zp = 0, zq = 0;
  verylong zr = 0, zt = 0, zu = 0, zv = 0;
  verylong zu1 = 0, zu2 = 0, zv1 = 0, zv2 = 0;

  *iterations = 0;
  zlshift(za, 1l, &zm);
  zmod(zb, zm, &zn);
  zsq(zb, &zd);
  zmul(za, zc, &zp);
  zlshift(zp, 2l, &zq);
  zsub(zd, zq, &zD);
  zsqrt(zD, &zd, &zp);
  zcopy(zb, &zu1);
  znegate(&zu1);
  zcopy(zm, &zu2);
  zone(&zv1);
  zzero(&zv2);
  zcopy(zb, &zp);
  zcopy(zm, &zq);
  L2:
  *iterations = *iterations + 1;
  zadd(zp, zd, &zf);
  zdiv(zf, zq, &zA, &zr);
  zmul(zA, zq, &zf);
  zsub(zf, zp, &zr);
  zcopy(zr, &zp);
  zsq(zp, &zf);
  zsub(zD, zf, &zg);
  zdiv(zg, zq, &zQ, &zr);
  zcopy(zQ, &zq);
  zmul(zA, zu2, &zf);
  zadd(zf, zu1, &zt);
  zcopy(zu2, &zu1);
  zcopy(zt, &zu2);
  zmul(zA, zv2, &zf);
  zadd(zf, zv1, &zt);
  zcopy(zv2, &zv1);
  zcopy(zt, &zv2);
  zmod(zp, zm, &zr);
  if (zcompare(zq, zm) == 0 && zcompare(zn, zr) == 0) {
    zdiv(zu2, za, &zu, &zr);
    zabs(&zu);
    zdiv(zv2, za, &zv, &zr);
    zabs(&zv);
    zmul(zv, zd, &zf);
    zadd(zu, zf, &zg);
    zrshift(zg, 1l, ze);
  }
  else goto L2;
  zfree(&zA);
  zfree(&zD);
  zfree(&zQ);
  zfree(&zd);
  zfree(&zf);
  zfree(&zg);
  zfree(&zm);
  zfree(&zn);
  zfree(&zp);
  zfree(&zq);
  zfree(&zr);
  zfree(&zt);
  zfree(&zu);
  zfree(&zv);
  zfree(&zu1);
  zfree(&zu2);
  zfree(&zv1);
  zfree(&zv2);
}

int main(void)
{
  long iterations;
  verylong za = 0, zb = 0, zc = 0, ze = 0;

  zintoz(2l, &za);
  zintoz(8l, &zb);
  zintoz(1l, &zc);
  fundamental(za, zb, zc, &ze, &iterations);
  printf("a = "); zwriteln(za);
  printf("b = "); zwriteln(zb);
  printf("c = "); zwriteln(zc);
  printf("e = "); zwriteln(ze);
  printf("iterations = %ld\n", iterations);
  zfree(&za);
  zfree(&zb);
  zfree(&zc);
  zfree(&ze);
  return 0;
}