/*
  Author:  Pate Williams (c) 1997

  "Algorithm 3.2.10 (Primitive Polynomial GCD). Given
  two polynomials A and B with coefficients in a UFD
  R, this algorithm computes a GCD of A and B, using
  only operations in R." -Henri Cohen- See "A Course
  in Computational Algebraic Number Theory" by Henri
  Cohen page 117.
*/

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

#define DEBUG
#define POLY_SIZE 8192l

void cont(long dA, verylong *zA, verylong *zc)
/* content of a polynomial */
{
  long i;
  verylong zd = 0;

  zgcd(zA[0], zA[1], zc);
  for (i = 2; i <= dA; i++) {
    zgcd(*zc, zA[i], &zd);
    zcopy(zd, zc);
  }
  zfree(&zd);
}

void pp(long dA, verylong *zA, verylong *zP)
/* primitive part of a polynomial */
{
  long i;
  verylong zc = 0, zr = 0;

  cont(dA, zA, &zc);
  for (i = 0; i <= dA; i++)
    zdiv(zA[i], zc, &zP[i], &zr);
  zfree(&zc);
  zfree(&zr);
}

void zspow(verylong za, long n, verylong *zb)
{
  int c = zscompare(za, 0l);
  verylong zc = 0;

  if (c > 0)
    zsexp(za, n, zb);
  else if (c < 0) {
    zcopy(za, &zc);
    zabs(&zc);
    zsexp(zc, n, zb);
    if (n & 1) znegate(zb);
  }
  else
    zzero(zb);
  zfree(&zc);
}

void zpoly_div(long m, long n, verylong *zu, verylong *zv,
               verylong *zq, verylong *zr, long *p, long *s)
{
  long j, jk, k, nk;
  verylong za = 0, zb = 0, zvn = 0;

  zcopy(zv[n], &zvn);
  for (j = 0; j <= m; j++)
    zcopy(zu[j], &zr[j]);
  if (m < n) {
    *p = 0, *s = m;
    zzero(&zq[0]);
  }
  else {
    *p = m - n, *s = n - 1;
    for (k = *p; k >= 0; k--) {
      nk = n + k;
      zspow(zvn, k, &za);
      zmul(zr[nk], za, &zq[k]);
      for (j = nk - 1; j >= 0; j--) {
        jk = j - k;
        if (jk >= 0) {
          zmul(zvn, zr[j], &za);
          zmul(zr[nk], zv[jk], &zb);
          zsub(za, zb, &zr[j]);
        }
        else {
          zcopy(zr[j], &za);
          zmul(zvn, za, &zr[j]);
        }
      }
    }
    while (*p > 0 && zscompare(zq[*p], 0l) == 0) *p = *p - 1;
    while (*s > 0 && zscompare(zr[*s], 0l) == 0) *s = *s - 1;
  }
  zfree(&za);
  zfree(&zb);
  zfree(&zvn);
}

void zwrite_poly(long dA, verylong *zA)
{
  long i;

  for (i = dA; i >= 0; i--) {
    zwrite(zA[i]);
    printf(" ");
  }
  printf("\n");
}

void primitive_poly_gcd(long dA, long dB,
                        verylong *zA,
                        verylong *zB,
                        long *dG)
/* returns the primitive polynomial GCD in B */
{
  int zero = 1;
  long dQ, dR, i;
  verylong za = 0, zb = 0, zc = 0, zd = 0, zr = 0;
  verylong zP[POLY_SIZE], zQ[POLY_SIZE], zR[POLY_SIZE];

  for (i = 0; zero && i <= dB; i++)
    zero = zscompare(zB[i], 0l) == 0;
  if (zero) {
    for (i = 0; i <= dA; i++)
      zcopy(zA[i], &zB[i]);
    *dG = dA;
    return;
  }
  cont(dA, zA, &za);
  cont(dB, zB, &zb);
  zgcd(za, zb, &zd);
  /* primitive part of A */
  for (i = 0; i <= dA; i++) {
    zdiv(zA[i], za, &zc, &zr);
    zcopy(zc, &zA[i]);
  }
  /* primitive part of B */
  for (i = 0; i <= dB; i++) {
    zdiv(zB[i], zb, &zc, &zr);
    zcopy(zc, &zB[i]);
  }
  for (i = 0; i < POLY_SIZE; i++)
   zP[i] = zQ[i] = zR[i] = 0;
  L2:
  zspow(zB[dB], dA - dB - 1, &za);
  for (i = 0; i <= dA; i++)
    zmul(zA[i], za, &zP[i]);
  zpoly_div(dA, dB, zP, zB, zQ, zR, &dQ, &dR);
  zero = 1;
  for (i = 0; zero && i <= dR; i++)
    zero = zscompare(zR[i], 0l) == 0;
  if (zero) goto L4;
  if (dR == 0) {
    dB = 0;
    zone(&zB[0]);
    goto L4;
  }
  for (i = 0; i <= dB; i++)
    zcopy(zB[i], &zA[i]);
  dA = dB;
  cont(dR, zR, &zr);
  /* primitive part of R */
  for (i = 0; i <= dR; i++)
    zdiv(zR[i], zr, &zB[i], &za);
  dB = dR;
  #ifdef DEBUG
  printf("A = ");
  zwrite_poly(dA, zA);
  printf("B = ");
  zwrite_poly(dB, zB);
  #endif
  goto L2;
  L4:
  for (i = 0; i <= dB; i++) {
    zmul(zd, zB[i], &zc);
    zcopy(zc, &zB[i]);
  }
  *dG = dB;
  zfree(&za);
  zfree(&zb);
  zfree(&zc);
  zfree(&zd);
  zfree(&zr);
  for (i = 0; i < POLY_SIZE; i++) {
    zfree(&zP[i]);
    zfree(&zQ[i]);
    zfree(&zR[i]);
  }
}

int main(void)
{
  long dA = 8, dB = 6, dG, i;
  long a[9] = {- 5, 2, 8, - 3, - 3, 0, 1, 0, 1};
  long b[7] = {21, - 9, - 4, 0, 5, 0, 3};
  verylong zA[POLY_SIZE], zB[POLY_SIZE];

  for (i = 0; i < POLY_SIZE; i++) zA[i] = zB[i] = 0;
  for (i = 0; i <= 6; i++) {
    zintoz(a[i], &zA[i]);
    zintoz(b[i], &zB[i]);
  }
  zintoz(a[7], &zA[7]);
  zintoz(a[8], &zA[8]);
  printf("A = ");
  zwrite_poly(dA, zA);
  printf("B = ");
  zwrite_poly(dB, zB);
  primitive_poly_gcd(dA, dB, zA, zB, &dG);
  printf("gcd(A, B) = ");
  zwrite_poly(dG, zB);
  for (i = 0; i < POLY_SIZE; i++) {
    zfree(&zA[i]);
    zfree(&zB[i]);
  }
  return 0;
}