/*
  Author:  Pate Williams (c) 1997

  Algorithm 5.4.2 (Reduction of Positive Definite
  Forms). See "A Course in Computational Algebraic
  Number Theory" by Henri Cohen page 243.
*/

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

void reduce(verylong *za, verylong *zb, verylong *zc)
{
  verylong zd = 0, ze = 0, zf = 0, zq = 0, zr = 0;

  zcopy(*za, &zd);
  znegate(&zd);
  if (zcompare(zd, *zb) < 0 && zcompare(*zb, *za) <= 0)
    goto L3;
  L2:
  zlshift(*za, 1l, &zd);
  zdiv(*zb, zd, &zq, &zr);
  if (zcompare(zr, *za) > 0) {
    zsub(zr, zd, &ze);
    zcopy(ze, &zr);
    zsadd(zq, 1l, &ze);
    zcopy(ze, &zq);
  }
  zadd(*zb, zr, &zd);
  zmul(zd, zq, &ze);
  zrshift(ze, 1l, &zf);
  zsub(*zc, zf, &zd);
  zcopy(zd, zc);
  zcopy(zr, zb);
  L3:
  if (zcompare(*za, *zc) > 0) {
    znegate(zb);
    zcopy(*za, &zd);
    zcopy(*zc, za);
    zcopy(zd, zc);
    goto L2;
  }
  if (zcompare(*za, *zc) == 0 && zscompare(*zb, 0l) < 0)
    znegate(zb);
  zfree(&zd);
  zfree(&ze);
  zfree(&zf);
  zfree(&zq);
  zfree(&zr);
}

int main(void)
{
  verylong zD = 0, za = 0, zb = 0, zc = 0, zd = 0;
  verylong ze = 0, zf = 0;

  zintoz(4l, &za);
  zintoz(6l, &zb);
  zintoz(8l, &zc);
  zsq(zb, &zd);
  zmul(za, zc, &ze);
  zsmul(ze, 4l, &zf);
  zsub(zd, zf, &zD);
  printf("a = "); zwriteln(za);
  printf("b = "); zwriteln(zb);
  printf("c = "); zwriteln(zc);
  printf("D = "); zwriteln(zD);
  reduce(&za, &zb, &zc);
  zsq(zb, &zd);
  zmul(za, zc, &ze);
  zsmul(ze, 4l, &zf);
  zsub(zd, zf, &zD);
  printf("a = "); zwriteln(za);
  printf("b = "); zwriteln(zb);
  printf("c = "); zwriteln(zc);
  printf("D = "); zwriteln(zD);
  zfree(&zD);
  zfree(&za);
  zfree(&zb);
  zfree(&zc);
  zfree(&zd);
  zfree(&ze);
  zfree(&zf);
  return 0;
}