/*
  Author:  Pate Williams (c) 1997

  "Algorithm 1.3.3 (Lehmer). Let a and b be
  non-negative multi-precision integers, and
  assume a >= b. This algorithm computes
  (a, b) = gcd(a, b)." -Henri Cohen- See "A
  Course in Computational Algebraic Number
  Theory" by Henri Cohen pages 13-14.
*/

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

#define DEBUG
#define MAX_DIGITS 4096

void gcd(verylong zA, verylong zB, verylong *za)
{
  verylong zb = 0, zr = 0;

  zcopy(zA, za);
  zcopy(zB, &zb);
  while (zscompare(zb, 0l) != 0) {
    zmod(*za, zb, &zr);
    zcopy(zb, za);
    zcopy(zr, &zb);
  }
  zfree(&zb);
  zfree(&zr);
}

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

  zcopy(zx, &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_gcd(long M, verylong zA, verylong zB,
                verylong *zg)
{
  long A, B, C, D, T, ah, bh, q;
  long xt, yt, xa[MAX_DIGITS], ya[MAX_DIGITS];
  verylong za = 0, zb = 0, zr = 0, zs = 0, zt = 0;
  verylong zu = 0;

  zcopy(zA, &za);
  zcopy(zB, &zb);
  L1:
  if (zscompare(zb, M) < 0) {
    #ifdef DEBUG
    printf("a = "); zwriteln(za);
    printf("b = "); zwriteln(zb);
    #endif
    gcd(za, zb, zg);
    zfree(&za);
    zfree(&zb);
    zfree(&zr);
    zfree(&zs);
    zfree(&zt);
    zfree(&zu);
    return;
  }
  radix_representation(M, za, xa, &xt);
  radix_representation(M, zb, ya, &yt);
  ah = xa[xt - 1];
  bh = ya[yt - 1];
  A = 1, B = 0, C = 0, D = 1;
  L2:
  if (xt == yt) {
    if (bh + C == 0 || bh + D == 0) goto L4;
    q = (ah + A) / (bh + C);
    if (q != (ah + B) / (bh + D)) goto L4;
    T = A - q * C, A = C, C = T;
    T = B - q * D, B = D, D = T;
    T = ah - q * bh, ah = bh, bh = T;
    #ifdef DEBUG
    printf("%3ld %3ld %2ld %2ld %2ld", ah, bh, A, B, C);
    printf(" %2ld %2ld\n", D, q);
    #endif
    goto L2;
  }
  L4:
  if (B == 0) {
    zmod(za, zb, &zt);
    zcopy(zb, &za);
    zcopy(zt, &zb);
  }
  else {
    #ifdef DEBUG
    printf("a = "); zwriteln(za);
    printf("b = "); zwriteln(zb);
    #endif
    zsmul(za, A, &zr);
    zsmul(zb, B, &zs);
    zadd(zr, zs, &zt);
    zsmul(za, C, &zr);
    zsmul(zb, D, &zs);
    zadd(zr, zs, &zu);
    zcopy(zt, &za);
    zcopy(zu, &zb);
    #ifdef DEBUG
    printf("a = "); zwriteln(za);
    printf("b = "); zwriteln(zb);
    #endif
  }
  goto L1;
}

int main(void)
{
  long M = 1000;
  verylong za = 0, zb = 0, zg = 0;

  zintoz(768454923l, &za);
  zintoz(542167814l, &zb);
  Lehmer_gcd(M, za, zb, &zg);
  printf("a = "); zwriteln(za);
  printf("b = "); zwriteln(zb);
  printf("gcd(a, b) = "); zwriteln(zg);
  zfree(&za);
  zfree(&zb);
  zfree(&zg);
  return 0;
}