/* $OpenXM: OpenXM_contrib2/asir2018/engine/gfspn.c,v 1.1 2018/09/19 05:45:07 noro Exp $ */

#include "ca.h"
#include "base.h"

UM current_mod_gfsn;

void setmod_gfsn(UM p)
{
  current_mod_gfsn = p;
}

void getmod_gfsn(UM *up)
{
  *up = current_mod_gfsn;
}

void simpgfsn(GFSN n,GFSN *r)
{
  UM t,q;

  if ( !n )
    *r = 0;
  else if ( NID(n) != N_GFSN )
    *r = n;
  else {
    t = UMALLOC(DEG(BDY(n)));
    q = W_UMALLOC(DEG(BDY(n)));
    cpyum(BDY(n),t);
    DEG(t) = divsfum(t,current_mod_gfsn,q);
    MKGFSN(t,*r);
  }
}

#define NZGFSN(a) ((a)&&(OID(a)==O_N)&&(NID(a)==N_GFSN))

void ntogfsn(Obj a,GFSN *b)
{
  UM t;
  GFS c;

  if ( !a || (OID(a)==O_N && NID(a) == N_GFSN) )
    *b = (GFSN)a;
  else if ( OID(a) == O_N && (NID(a) == N_GFS || NID(a) == N_Q) ) {
    ntogfs((Obj)a,&c);
    if ( !b )
      *b = 0;
    else {
      t = UMALLOC(0); ptosfum((P)c,t);
      MKGFSN(t,*b);
    }
  } else
    error("ntogfsn : invalid argument");
}

void addgfsn(GFSN a,GFSN b,GFSN *c)
{
  UM t,q;
  GFSN z;
  int d;

  ntogfsn((Obj)a,&z); a = z; ntogfsn((Obj)b,&z); b = z;
  if ( !a )
    *c = b;
  else if ( !b )
    *c = a;
  else {
    d = MAX(DEG(BDY(a)),DEG(BDY(b)));
    t = UMALLOC(d);
    q = W_UMALLOC(d);
    addsfum(BDY(a),BDY(b),t);
    DEG(t) = divsfum(t,current_mod_gfsn,q);
    MKGFSN(t,z);
    *c = (GFSN)z;
  }
}

void subgfsn(GFSN a,GFSN b,GFSN *c)
{
  UM t,q;
  GFSN z;
  int d;

  ntogfsn((Obj)a,&z); a = z; ntogfsn((Obj)b,&z); b = z;
  if ( !a )
    chsgngfsn(b,c);
  else if ( !b )
    *c = a;
  else {
    d = MAX(DEG(BDY(a)),DEG(BDY(b)));
    t = UMALLOC(d);
    q = W_UMALLOC(d);
    subsfum(BDY(a),BDY(b),t);
    DEG(t) = divsfum(t,current_mod_gfsn,q);
    MKGFSN(t,*c);
  }
}

extern int up_lazy;

void mulgfsn(GFSN a,GFSN b,GFSN *c)
{
  UM t,q;
  GFSN z;
  int d;

  ntogfsn((Obj)a,&z); a = z; ntogfsn((Obj)b,&z); b = z;
  if ( !a || !b )
    *c = 0;
  else {
    d = DEG(BDY(a))+DEG(BDY(b));
    t = UMALLOC(d);
    q = W_UMALLOC(d);
    mulsfum(BDY(a),BDY(b),t);
    DEG(t) = divsfum(t,current_mod_gfsn,q);
    MKGFSN(t,*c);
  }
}

void divgfsn(GFSN a,GFSN b,GFSN *c)
{
  GFSN z;
  int d;
  UM wb,wc,wd,we,t,q;

  ntogfsn((Obj)a,&z); a = z; ntogfsn((Obj)b,&z); b = z;
  if ( !b )
    error("divgfsn: division by 0");
  else if ( !a )
    *c = 0;
  else {
    wb = W_UMALLOC(DEG(BDY(b))); cpyum(BDY(b),wb); 
    d = DEG(current_mod_gfsn);
    wc = W_UMALLOC(d); cpyum(current_mod_gfsn,wc);
    wd = W_UMALLOC(2*d); we = W_UMALLOC(2*d);
    /* wd*wb+we*wc=1 */
    eucsfum(wb,wc,wd,we);
    d = DEG(BDY(a))+DEG(wd);
    t = UMALLOC(d);
    q = W_UMALLOC(d);
    mulsfum(BDY(a),wd,t);
    DEG(t) = divsfum(t,current_mod_gfsn,q);
    MKGFSN(t,*c);
  }
}

void invgfsn(GFSN b,GFSN *c)
{
  GFSN z;
  int d;
  UM wb,wc,wd,we,t;

  ntogfsn((Obj)b,&z); b = z;
  if ( !b )
    error("divgfsn: division by 0");
  else {
    wb = W_UMALLOC(DEG(BDY(b))); cpyum(BDY(b),wb); 
    d = DEG(current_mod_gfsn);
    wc = W_UMALLOC(d); cpyum(current_mod_gfsn,wc);
    wd = W_UMALLOC(2*d); we = W_UMALLOC(2*d);
    /* wd*wb+we*wc=1 */
    eucsfum(wb,wc,wd,we);
    d = DEG(wd);
    t = UMALLOC(d);
    cpyum(wd,t);
    MKGFSN(t,*c);
  }
}

void chsgngfsn(GFSN a,GFSN *c)
{
  GFSN z;
  int d;
  struct oUM zero;
  UM t;

  ntogfsn((Obj)a,&z); a = z;
  if ( !a )
    *c = 0;
  else {
    d = DEG(BDY(a));
    t = UMALLOC(d);
    DEG(&zero) = -1;
    subsfum(&zero,BDY(a),t);
    MKGFSN(t,*c);
  }
}

void pwrgfsn(GFSN a,Z b,GFSN *c)
{
  GFSN z;
  UM t,x,y,q;
  int d,k,blen;

  ntogfsn((Obj)a,&z); a = z;
  if ( !b ) {
    t = UMALLOC(0); DEG(t) = 0; COEF(t)[0] = _onesf();
    MKGFSN(t,*c);
  } else if ( !a )
    *c = 0;
  else {
    d = DEG(current_mod_gfsn);

    /* y = 1 */
    y = UMALLOC(d); DEG(y) = 0; COEF(y)[0] = _onesf();

    t = W_UMALLOC(2*d); q = W_UMALLOC(2*d);

    /* x = simplify(a) */
    x = W_UMALLOC(DEG(BDY(a))); cpyum(BDY(a),x);
    DEG(x) = divsfum(x,current_mod_gfsn,q);
    if ( DEG(x) < 0 ) {
      *c = 0;
    } else {
      blen = z_bits((Q)b);
      for ( k = blen-1; k >= 0; k-- ) {
        mulsfum(y,y,t);
        DEG(t) = divsfum(t,current_mod_gfsn,q);
        cpyum(t,y);
        if ( tstbitz(b,k) ) {
          mulsfum(y,x,t);
          DEG(t) = divsfum(t,current_mod_gfsn,q);
          cpyum(t,y);
        }
      }
      MKGFSN(y,*c);
    }
  }
}

int cmpgfsn(GFSN a,GFSN b)
{
  GFSN z;

  ntogfsn((Obj)a,&z); a = z; ntogfsn((Obj)b,&z); b = z;
  if ( !a )
    if ( !b )
      return 0;
    else
      return -1;
  else if ( !b )
      return 1;
  else
    return compsfum(BDY(a),BDY(b));
}

void randomgfsn(GFSN *r)
{
  int d;
  UM t;

  if ( !current_mod_gfsn )
    error("randomgfsn : current_mod_gfsn is not set");
  d = DEG(current_mod_gfsn);
  t = UMALLOC(d-1);
  randsfum(d,t);
  degum(t,d-1);
  if ( DEG(t) < 0 )
    *r = 0;
  else {
    MKGFSN(t,*r);
  }
}