/*
  Author:  Pate Williams (c) 1997

  Algorithm 3.5.5 (Hensel Lift). See "A Course in
  Computational Algebraic Number Theory" by Henri
  Cohen page 138.
*/

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

#define POLY_SIZE 8192l

long gcd(long a, long b)
{
  long r;

  if (a < 0) a = - a;
  if (b < 0) b = - b;
  while (b > 0)
    r = a % b, a = b, b = r;
  return a;
}

void zpoly_copy(long dA, verylong *zA, verylong *zB, long *dB)
{
  long i;

  *dB = dA;
  for (i = 0; i <= dA; i++) zcopy(zA[i], &zB[i]);
}

void zpoly_add(long dA, long dB, verylong *zA,
               verylong *zB, verylong *zC, long *dC)
{
  int zero;
  long i;

  if (dA >= dB) {
    for (i = 0; i <= dB; i++)
      zadd(zA[i], zB[i], &zC[i]);
    for (i = dB + 1; i <= dA; i++)
      zcopy(zA[i], &zC[i]);
    *dC = dA;
  }
  else {
    for (i = 0; i <= dA; i++)
      zadd(zA[i], zB[i], &zC[i]);
    for (i = dA + 1; i <= dB; i++)
      zcopy(zB[i], &zC[i]);
    *dC = dB;
  }
  zero = zscompare(zC[*dC], 0l) == 0;
  i = *dC;
  while (zero && i >= 0 && *dC > 0)
    *dC = *dC - 1, i--, zero = zscompare(zC[i], 0l) == 0;
}

void zpoly_sub(long dA, long dB, verylong *zA,
               verylong *zB, verylong *zC, long *dC)
{
  int zero;
  long i;

  if (dA >= dB) {
    for (i = 0; i <= dB; i++)
      zsub(zA[i], zB[i], &zC[i]);
    for (i = dB + 1; i <= dA; i++)
      zcopy(zA[i], &zC[i]);
    *dC = dA;
  }
  else {
    for (i = 0; i <= dA; i++)
      zsub(zA[i], zB[i], &zC[i]);
    for (i = dA + 1; i <= dB; i++)
      zcopy(zB[i], &zC[i]);
    for (i = dA + 1; i <= dB; i++)
      znegate(&zC[i]);
    *dC = dB;
  }
  zero = zscompare(zC[*dC], 0l) == 0;
  i = *dC;
  while (zero && i >= 0 && *dC > 0)
    *dC = *dC - 1, i--, zero = zscompare(zC[i], 0l) == 0;
}

void zpoly_mul(long m, long n, verylong *za, verylong *zb,
               verylong *zc, long *p)
{
  int zero;
  long i, j, k;
  verylong zd = 0, zai = 0, zbk = 0, zsum = 0, zterm = 0;

  *p = m + n;
  for (k = 0; k <= *p; k++) {
    zzero(&zsum);
    for (i = 0; i <= k; i++) {
      j = k - i;
      if (i > m) zzero(&zai);
      else zcopy(za[i], &zai);
      if (j > n) zzero(&zbk);
      else zcopy(zb[j], &zbk);
      zmul(zai, zbk, &zterm);
      zcopy(zsum, &zd);
      zadd(zterm, zd, &zsum);
    }
    zcopy(zsum, &zc[k]);
  }
  zero = zscompare(zc[*p], 0l) == 0;
  i = *p;
  while (zero && i >= 0 && *p > 0)
    *p = *p - 1, i--, zero = zscompare(zc[i], 0l) == 0;
  zfree(&zd);
  zfree(&zai);
  zfree(&zbk);
  zfree(&zsum);
  zfree(&zterm);
}

void zpoly_euclid_div(long dA, long dB, verylong *zA,
                      verylong *zB, verylong *zQ,
                      verylong *zR, long *dQ, long *dR)
{
  long dC, dD, dS, i;
  verylong zr = 0;
  verylong zC[POLY_SIZE], zD[POLY_SIZE], zS[POLY_SIZE];

  for (i = 0; i < POLY_SIZE; i++) zC[i] = zD[i] = zS[i] = 0;
  for (i = 0; i <= dA; i++) zR[i] = zA[i];
  *dQ = 0;
  *dR = dA;
  zzero(&zQ[0]);
  while (*dR >= dB) {
    dS = *dR - dB;
    zdiv(zR[*dR], zB[dB], &zS[dS], &zr);
    zpoly_add(*dQ, dS, zQ, zS, zC, &dC);
    zpoly_copy(dC, zC, zQ, dQ);
    zpoly_mul(dS, dB, zS, zB, zC, &dC);
    zpoly_sub(*dR, dC, zR, zC, zD, &dD);
    zpoly_copy(dD, zD, zR, dR);
  }
  zfree(&zr);
  for (i = 0; i < POLY_SIZE; i++) {
    zfree(&zC[i]);
    zfree(&zD[i]);
    zfree(&zS[i]);
  }
}

void zpoly_lift(long dA, long r, verylong *zA)
{
  long i;
  verylong za = 0;

  if (zscompare(zA[dA], 0l) < 0) {
    for (i = 0; i <= dA; i++) {
      zsadd(zA[i], r, &za);
      zcopy(za, &zA[i]);
    }
  }
  else {
    for (i = 0; i <= dA; i++) {
      zsadd(zA[i], - r, &za);
      zcopy(za, &zA[i]);
    }
  }
  zfree(&za);
}

void Hensel_lift(long dA, long dB, long dC,
                 long dU, long dV, long p, long q,
                 verylong *zA, verylong *zB,
                 verylong *zC, verylong *zU,
                 verylong *zV, verylong *zA1,
                 verylong *zB1, long *dA1,
                 long *dB1)
{
  long dD, dE, dQ, dR, df, i, r = gcd(p, q);
  verylong zq = 0;
  verylong zD[POLY_SIZE], zE[POLY_SIZE];
  verylong zQ[POLY_SIZE], zR[POLY_SIZE];
  verylong zf[POLY_SIZE];

  for (i = 0; i < POLY_SIZE; i++)
    zD[i] = zE[i] = zQ[i] = zR[i] = zf[i] = 0;
  zpoly_mul(dA, dB, zA, zB, zD, &dD);
  zpoly_sub(dC, dD, zC, zD, zf, &df);
  for (i = 0; i <= df; i++) {
    zsdiv(zf[i], q, &zq);
    zintoz(zsmod(zq, r), &zf[i]);
  }
  zpoly_mul(dV, df, zV, zf, zD, &dD);
  zpoly_euclid_div(dD, dA, zD, zA, zQ, zR, &dQ, &dR);
  zpoly_lift(dR, r, zR);
  for (i = 0; i <= dR; i++) {
    zsmul(zR[i], q, &zq);
    zcopy(zq, &zR[i]);
  }
  zpoly_add(dA, dR, zA, zR, zA1, dA1);
  zpoly_mul(dU, df, zU, zf, zD, &dD);
  zpoly_copy(dB, zB, zE, &dE);
  for (i = 0; i <= dE; i++) znegate(&zE[i]);
  zpoly_euclid_div(dD, dE, zD, zE, zQ, zR, &dQ, &dR);
  zpoly_lift(dR, r, zR);
  for (i = 0; i <= dR; i++) {
    zsmul(zR[i], q, &zq);
    zcopy(zq, &zR[i]);
  }
  zpoly_add(dB, dR, zB, zR, zB1, dB1);
  zfree(&zq);
  for (i = 0; i < POLY_SIZE; i++) {
    zfree(&zD[i]);
    zfree(&zE[i]);
    zfree(&zQ[i]);
    zfree(&zR[i]);
    zfree(&zf[i]);
  }
}

int main(void)
{
  long dA = 1, dB = 1, dC = 2, dU = 0, dV = 0, dA1, dB1;
  long i, p = 33, q = 9;
  verylong zA[POLY_SIZE], zB[POLY_SIZE];
  verylong zC[POLY_SIZE], zU[POLY_SIZE];
  verylong zV[POLY_SIZE], zA1[POLY_SIZE];
  verylong zB1[POLY_SIZE];

  for (i = 0; i < POLY_SIZE; i++)
    zA[i] = zB[i] = zC[i] = zU[i] = zV[i] = zA1[i] = zB1[i] = 0;
  zone(&zA[1]);
  zintoz(- 3l, &zA[0]);
  zone(&zB[1]);
  zintoz(- 4l, &zB[0]);
  zone(&zC[2]);
  zintoz(2l, &zC[1]);
  zintoz(3l, &zC[0]);
  zone(&zU[0]);
  zone(&zV[0]);
  znegate(&zV[0]);
  Hensel_lift(dA, dB, dC, dU, dV, p, q, zA, zB, zC, zU, zV,
              zA1, zB1, &dA1, &dB1);
  printf("A1 = ");
  for (i = dA1; i >= 0; i--)
    printf("%3d ", ztoint(zA1[i]));
  printf("\n");
  printf("B1 = ");
  for (i = dB1; i >= 0; i--)
    printf("%3d ", ztoint(zB1[i]));
  printf("\n");
  for (i = 0; i < POLY_SIZE; i++) {
    zfree(&zA[i]);
    zfree(&zB[i]);
    zfree(&zC[i]);
    zfree(&zU[i]);
    zfree(&zV[i]);
    zfree(&zA1[i]);
    zfree(&zB1[i]);
  }
  return 0;
}