/*
  Author:  Pate Williams (c) 1997

  Algorithm 1.3.7 (Lehmer Extended). Let a and b
  be non-negative integers, such that a >= b.
  This algorithm computes (u, v, d) such that
  a * u + b * v = d = (a, b) = gcd(a, b).
  See "A Course in Computational Algebraic Number
  Theory" by Henri Cohen page 17.
*/

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

void Euclid_extended(verylong za, verylong zb,
                     verylong *zx, verylong *zy,
                     verylong *zd)
{
  verylong zA = 0, zB = 0, ze = 0, zq = 0, zr = 0;
  verylong zx1 = 0, zx2 = 0, zy1 = 0, zy2 = 0;

  if (zscompare(zb, 0l) == 0) {
    zcopy(za, zd);
    zone(zx);
    zzero(zy);
  }
  else {
    zcopy(za, &zA);
    zcopy(zb, &zB);
    zone(&zx2);
    zzero(&zx1);
    zzero(&zy2);
    zone(&zy1);
    while (zscompare(zB, 0l) > 0) {
      zdiv(zA, zB, &zq, &zr);
      zmul(zq, zx1, &ze);
      zsub(zx2, ze, zx);
      zmul(zq, zy1, &ze);
      zsub(zy2, ze, zy);
      zcopy(zB, &zA);
      zcopy(zr, &zB);
      zcopy(zx1, &zx2);
      zcopy(*zx, &zx1);
      zcopy(zy1, &zy2);
      zcopy(*zy, &zy1);
    }
    zcopy(zA, zd);
    zcopy(zx2, zx);
    zcopy(zy2, zy);
  }
  zfree(&zA);
  zfree(&zB);
  zfree(&ze);
  zfree(&zq);
  zfree(&zr);
  zfree(&zx1);
  zfree(&zx2);
  zfree(&zy1);
  zfree(&zy2);
}

void radix_representation(long b, verylong zx,
                          long *a, long *t)
{
  long i = 0;
  verylong za = 0, zq = 0, zr = 0;

  zcopy(zx, &za);
  zabs(&za);
  a[i++] = zsdiv(za, b, &zq);
  while (zscompare(zq, 0l) > 0) {
    zcopy(zq, &za);
    a[i++] = zsdiv(zq, b, &zr);
    zcopy(zr, &zq);
  }
  *t = i;
  zfree(&za);
  zfree(&zq);
  zfree(&zr);
}

void Lehmer_extended(long b, verylong za, verylong zb,
                     verylong *zu, verylong *zv,
                     verylong *zd)
{
  long A, B, C, D, T, ah, bh, q, qp;
  long aa[4096], at, ba[4096], bt;
  verylong ze = 0, zf = 0, zq = 0, zr = 0, zt = 0;
  verylong zv1 = 0;

  zone(zu);
  zzero(&zv1);
  while (zscompare(zb, b) > 0) {
    radix_representation(b, za, aa, &at);
    radix_representation(b, zb, ba, &bt);
    ah = aa[at - 1];
    bh = ba[bt - 1];
    A = D = 1, B = C = 0;
    if (at == bt) {
      while (bh + C != 0 && bh + D != 0) {
        q  = (ah + A) / (bh + C);
        qp = (ah + B) / (bh + D);
        if (q != qp) break;
        T = A - q * C, A = C, C = T;
        T = B - q * D, B = D, D = T;
        T = ah - q * bh, ah = bh, bh = T;
      }
    }
    if (B == 0) {
      zdiv(za, zb, &zq, &zt);
      zcopy(zb, &za);
      zcopy(zt, &zb);
      zmul(zq, zv1, &ze);
      zsub(*zu, ze, &zt);
      zcopy(zv1, zu);
      zcopy(zt, &zv1);
    }
    else {
      zsmul(za, A, &ze);
      zsmul(zb, B, &zf);
      zadd(ze, zf, &zt);
      zsmul(za, C, &ze);
      zsmul(zb, D, &zf);
      zadd(ze, zf, &zr);
      zcopy(zt, &za);
      zcopy(zr, &zb);
      zsmul(*zu, A, &ze);
      zsmul(zv1, B, &zf);
      zadd(ze, zf, &zt);
      zsmul(*zu, C, &ze);
      zsmul(zv1, D, &zf);
      zadd(ze, zf, &zr);
      zcopy(zt, zu);
      zcopy(zr, &zv1);
    }
  }
  Euclid_extended(za, zb, zu, zv, zd);
  zfree(&ze);
  zfree(&zf);
  zfree(&zq);
  zfree(&zr);
  zfree(&zt);
  zfree(&zv1);
}

int main(void)
{
  long a[5] = {4864l, 5000l, 7500, 768454923l, 2000000000l};
  long b[5] = {3458l, 2600l,  250, 542167814l, 1000000000l};
  long i;
  verylong zA = 0, zB = 0;
  verylong za = 0, zb = 0, zd = 0, zu = 0, zv = 0;
  verylong ze = 0, zf = 0, zg = 0;

  for (i = 0; i < 6; i++) {
    if (i < 5) {
      zintoz(a[i], &za);
      zintoz(b[i], &zb);
    }
    else {
      printf("a = "); zread(&za);
      printf("b = "); zread(&zb);
    }
    Euclid_extended(za, zb, &zu, &zv, &zd);
    printf("a = "); zwriteln(za);
    printf("b = "); zwriteln(zb);
    printf("u = "); zwriteln(zu);
    printf("v = "); zwriteln(zv);
    printf("d = "); zwriteln(zd);
    zmul(za, zu, &ze);
    zmul(zb, zv, &zf);
    zadd(ze, zf, &zg);
    if(zcompare(zd, zg) != 0)
      printf("*error*\nin Euclid_extended\n");
    zcopy(za, &zA);
    zcopy(zb, &zB);
    Lehmer_extended(1000l, za, zb, &zu, &zv, &zd);
    printf("a = "); zwriteln(za);
    printf("b = "); zwriteln(zb);
    printf("u = "); zwriteln(zu);
    printf("v = "); zwriteln(zv);
    printf("d = "); zwriteln(zd);
    zmul(za, zu, &ze);
    zmul(zb, zv, &zf);
    zadd(ze, zf, &zg);
    if(zcompare(zd, zg) != 0)
      printf("*error*\nin Lehmer_extended\n");
    zcopy(zA, &za);
    zcopy(zB, &zb);
    zexteucl(za, &zu, zb, &zv, &zd);
    printf("a = "); zwriteln(za);
    printf("b = "); zwriteln(zb);
    printf("u = "); zwriteln(zu);
    printf("v = "); zwriteln(zv);
    printf("d = "); zwriteln(zd);
    zmul(za, zu, &ze);
    zmul(zb, zv, &zf);
    zadd(ze, zf, &zg);
    if(zcompare(zd, zg) != 0)
      printf("*error*\nin exteucl\n");
  }
  zfree(&zA);
  zfree(&zB);
  zfree(&za);
  zfree(&zb);
  zfree(&zd);
  zfree(&ze);
  zfree(&zf);
  zfree(&zg);
  zfree(&zu);
  zfree(&zv);
  return 0;
}