/*
  Author:  Pate Williams (c) 1997

  "Algorithm 3.3.7 (Sub-Resultant). Given two
  polynomials with coefficients in a UFD R, this
  algorithm computes the resultant of A and B."
  -Henri Cohen- See "A Course in Computational
  Algebraic Number Theory" by Henri Cohen page
  122.
*/

#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 zpoly_write(char *label, long dA, verylong *zA)
{
  long i;

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

void sub_resultant(long dA, long dB,
                   verylong *zA,
                   verylong *zB,
                   verylong *zu)
/* returns the sub-resultant in u */
{
  int zero_A = 1, zero_B = 1;
  long dQ, dR, delta, i, t;
  verylong za = 0, zb = 0, zg = 0, zh = 0, zr = 0;
  verylong zs = 0, zt = 0, zx = 0, zy = 0, zz = 0;
  verylong zQ[POLY_SIZE] = {0};
  verylong zR[POLY_SIZE] = {0};

  for (i = 0; zero_A && i <= dA; i++)
    zero_A = zscompare(zA[i], 0l) == 0;
  for (i = 0; zero_B && i <= dB; i++)
    zero_B = zscompare(zB[i], 0l) == 0;
  if (zero_A || zero_B) {
    zzero(zu);
    return;
  }
  cont(dA, zA, &za);
  cont(dB, zB, &zb);
  /* primitive part of A */
  for (i = 0; i <= dA; i++) {
    zdiv(zA[i], za, &zx, &zr);
    zcopy(zx, &zA[i]);
  }
  /* primitive part of B */
  for (i = 0; i <= dB; i++) {
    zdiv(zB[i], zb, &zx, &zr);
    zcopy(zx, &zB[i]);
  }
  zone(&zg);
  zone(&zh);
  zone(&zs);
  zsexp(za, dB, &zx);
  zsexp(zb, dA, &zy);
  zmul(zx, zy, &zt);
  if (dA < dB) {
    /* exchange polynomials A and B */
    for (i = 0; i <= dA; i++)
      zcopy(zA[i], &zQ[i]);
    for (i = 0; i <= dB; i++)
      zcopy(zB[i], &zA[i]);
    for (i = 0; i <= dA; i++)
      zcopy(zQ[i], &zB[i]);
    t = dA, dA = dB, dB = t;
    if (dA & 1 && dB & 1) znegate(&zs);
  }
  L2:
  delta = dA - dB;
  zpoly_div(dA, dB, zA, zB, zQ, zR, &dQ, &dR);
  for (i = 0; i <= dB; i++)
    zcopy(zB[i], &zA[i]);
  dA = dB;
  zspow(zh, delta, &zx);
  zmul(zg, zx, &zy);
  for (i = 0; i <= dR; i++)
    zdiv(zR[i], zy, &zB[i], &zr);
  dB = dR;
  zcopy(zA[dA], &zg);
  zspow(zg, delta, &zx);
  if (1 - delta < 0) {
    zspow(zh, delta - 1, &zy);
    zdiv(zx, zy, &zh, &zr);
  }
  else {
    zspow(zh, 1 - delta, &zy);
    zmul(zx, zy, &zh);
  }
  #ifdef DEBUG
  zpoly_write("A = ", dA, zA);
  zpoly_write("B = ", dB, zB);
  zpoly_write("R = ", dR, zR);
  printf("g = "); zwriteln(zg);
  printf("h = "); zwriteln(zh);
  #endif
  if (dB > 0) goto L2;
  zspow(zB[dB], dA, &zx);
  if (1 - dA < 0) {
    zspow(zh, dA - 1, &zy);
    zdiv(zx, zy, &zz, &zr);
  }
  else {
    zspow(zh, 1 - dA, &zy);
    zmul(zx, zy, &zz);
  }
  zmul(zs, zt, &zx);
  zmul(zx, zz, zu);
  zfree(&za);
  zfree(&zb);
  zfree(&zr);
  zfree(&zs);
  zfree(&zt);
  zfree(&zx);
  zfree(&zy);
  for (i = 0; i < POLY_SIZE; i++) {
    zfree(&zQ[i]);
    zfree(&zR[i]);
  }
}

int main(void)
{
  long dA = 3, dB = 2, i;
  long a[4] = {- 1, - 1, - 1, 1}, b[3] = {- 1, - 2, 3};
  verylong zu = 0;
  verylong zA[POLY_SIZE] = {0};
  verylong zB[POLY_SIZE] = {0};

  for (i = 0; i <= dA; i++)
    zintoz(a[i], &zA[i]);
  for (i = 0; i <= dB; i++)
    zintoz(b[i], &zB[i]);
  zpoly_write("A = ", dA, zA);
  zpoly_write("B = ", dB, zB);
  sub_resultant(dA, dB, zA, zB, &zu);
  printf("sub-resultant = ");
  zwriteln(zu);
  zfree(&zu);
  for (i = 0; i < POLY_SIZE; i++) {
    zfree(&zA[i]);
    zfree(&zB[i]);
  }
  return 0;
}