/*
  Author:  Pate Williams (c) 1997

  14.57 Algorithm Lehmer's gcd algorithm
  See "Handbook of Applied Cryptography" by Alfred
  J. Menezes et al page 607.
*/

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

#define DEBUG

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 b, verylong zx, verylong zy,
                verylong *zv)
{
  long A, B, C, D, q, qp, t;
  long xp, xt, yp, yt, xa[4096], ya[4096];
  verylong zT = 0, za = 0, zb = 0, zc = 0, zd = 0;
  verylong zu = 0;

  zcopy(zx, &za);
  zcopy(zy, &zb);
  while (zscompare(zb, b) > 0) {
    radix_representation(b, za, xa, &xt);
    radix_representation(b, zb, ya, &yt);
    xp = xa[xt - 1];
    yp = ya[yt - 1];
    A = D = 1, B = C = 0;
    if (xt == yt) {
      while (yp + C != 0 && yp + D != 0) {
        q  = (xp + A) / (yp + C);
        qp = (xp + B) / (yp + D);
        #ifdef DEBUG
        printf("%3ld %3ld %2ld %2ld %2ld %2ld %2ld %2ld\n",
               xp, yp, A, B, C, D, q, qp);
        #endif
        if (q != qp) break;
        t = A - q * C, A = C, C = t;
        t = B - q * D, B = D, D = t;
        t = xp - q * yp, xp = yp, yp = t;
      }
    }
    if (B == 0) {
      zmod(za, zb, &zT);
      zcopy(zb, &za);
      zcopy(zT, &zb);
    }
    else {
      zsmul(za, A, &zc);
      zsmul(zb, B, &zd);
      zadd(zc, zd, &zT);
      zsmul(za, C, &zc);
      zsmul(zb, D, &zd);
      zadd(zc, zd, &zu);
      zcopy(zT, &za);
      zcopy(zu, &zb);
    }
    #ifdef DEBUG
    zwrite(za); printf(" "); zwriteln(zb);
    #endif
  }
  while (zcompare(zb, 0l) > 0) {
    zmod(za, zb, &zu);
    zcopy(zb, &za);
    zcopy(zu, &zb);
    #ifdef DEBUG
    zwrite(za); printf(" "); zwriteln(zb);
    #endif
  }
  zcopy(za, zv);
  zfree(&zT);
  zfree(&za);
  zfree(&zb);
  zfree(&zc);
  zfree(&zd);
  zfree(&zu);
}

int main(void)
{
  long b = 1000;
  verylong zv = 0, zx = 0, zy = 0;

  zintoz(768454923l, &zx);
  zintoz(542167814l, &zy);
  Lehmer_gcd(b, zx, zy, &zv);
  printf("x = "); zwriteln(zx);
  printf("y = "); zwriteln(zy);
  printf("gcd(x, y) = "); zwriteln(zv);
  zfree(&zv);
  zfree(&zx);
  zfree(&zy);
  return 0;
}