/*
 * Copyright (c) 1994-2000 FUJITSU LABORATORIES LIMITED 
 * All rights reserved.
 * 
 * FUJITSU LABORATORIES LIMITED ("FLL") hereby grants you a limited,
 * non-exclusive and royalty-free license to use, copy, modify and
 * redistribute, solely for non-commercial and non-profit purposes, the
 * computer program, "Risa/Asir" ("SOFTWARE"), subject to the terms and
 * conditions of this Agreement. For the avoidance of doubt, you acquire
 * only a limited right to use the SOFTWARE hereunder, and FLL or any
 * third party developer retains all rights, including but not limited to
 * copyrights, in and to the SOFTWARE.
 * 
 * (1) FLL does not grant you a license in any way for commercial
 * purposes. You may use the SOFTWARE only for non-commercial and
 * non-profit purposes only, such as academic, research and internal
 * business use.
 * (2) The SOFTWARE is protected by the Copyright Law of Japan and
 * international copyright treaties. If you make copies of the SOFTWARE,
 * with or without modification, as permitted hereunder, you shall affix
 * to all such copies of the SOFTWARE the above copyright notice.
 * (3) An explicit reference to this SOFTWARE and its copyright owner
 * shall be made on your publication or presentation in any form of the
 * results obtained by use of the SOFTWARE.
 * (4) In the event that you modify the SOFTWARE, you shall notify FLL by
 * e-mail at risa-admin@sec.flab.fujitsu.co.jp of the detailed specification
 * for such modification or the source code of the modified part of the
 * SOFTWARE.
 * 
 * THE SOFTWARE IS PROVIDED AS IS WITHOUT ANY WARRANTY OF ANY KIND. FLL
 * MAKES ABSOLUTELY NO WARRANTIES, EXPRESSED, IMPLIED OR STATUTORY, AND
 * EXPRESSLY DISCLAIMS ANY IMPLIED WARRANTY OF MERCHANTABILITY, FITNESS
 * FOR A PARTICULAR PURPOSE OR NONINFRINGEMENT OF THIRD PARTIES'
 * RIGHTS. NO FLL DEALER, AGENT, EMPLOYEES IS AUTHORIZED TO MAKE ANY
 * MODIFICATIONS, EXTENSIONS, OR ADDITIONS TO THIS WARRANTY.
 * UNDER NO CIRCUMSTANCES AND UNDER NO LEGAL THEORY, TORT, CONTRACT,
 * OR OTHERWISE, SHALL FLL BE LIABLE TO YOU OR ANY OTHER PERSON FOR ANY
 * DIRECT, INDIRECT, SPECIAL, INCIDENTAL, PUNITIVE OR CONSEQUENTIAL
 * DAMAGES OF ANY CHARACTER, INCLUDING, WITHOUT LIMITATION, DAMAGES
 * ARISING OUT OF OR RELATING TO THE SOFTWARE OR THIS AGREEMENT, DAMAGES
 * FOR LOSS OF GOODWILL, WORK STOPPAGE, OR LOSS OF DATA, OR FOR ANY
 * DAMAGES, EVEN IF FLL SHALL HAVE BEEN INFORMED OF THE POSSIBILITY OF
 * SUCH DAMAGES, OR FOR ANY CLAIM BY ANY OTHER PARTY. EVEN IF A PART
 * OF THE SOFTWARE HAS BEEN DEVELOPED BY A THIRD PARTY, THE THIRD PARTY
 * DEVELOPER SHALL HAVE NO LIABILITY IN CONNECTION WITH THE USE,
 * PERFORMANCE OR NON-PERFORMANCE OF THE SOFTWARE.
 *
 * $OpenXM: OpenXM_contrib2/asir2018/builtin/fctr.c,v 1.2 2018/09/28 08:20:27 noro Exp $
*/
#include "ca.h"
#include "parse.h"

void Pfctr(), Pgcd(), Pgcdz(), Plcm(), Psqfr(), Pufctrhint();
void Pptozp(), Pcont(), Psfcont();
void Pafctr(), Pagcd();
void Pmodsqfr(),Pmodfctr(),Pddd(),Pnewddd(),Pddd_tab();
void Psfsqfr(),Psffctr(),Psfbfctr(),Psfufctr(),Psfmintdeg(),Psfgcd();
void Pirred_check(), Pnfctr_mod();
void Pbivariate_hensel_special();

void sfmintdeg(VL vl,P fx,int dy,int c,P *fr);

struct ftab fctr_tab[] = {
  {"bivariate_hensel_special",Pbivariate_hensel_special,6},
  {"fctr",Pfctr,-2},
  {"gcd",Pgcd,-3},
  {"gcdz",Pgcdz,2},
  {"lcm",Plcm,2},
  {"sqfr",Psqfr,1},
  {"ufctrhint",Pufctrhint,2},
  {"ptozp",Pptozp,1},
  {"cont",Pcont,-2},
  {"sfcont",Psfcont,-2},
  {"afctr",Pafctr,2},
  {"agcd",Pagcd,3},
  {"modsqfr",Pmodsqfr,2},
  {"modfctr",Pmodfctr,2},
  {"sfsqfr",Psfsqfr,1},
  {"sffctr",Psffctr,1},
  {"sfufctr",Psfufctr,1},
  {"sfbfctr",Psfbfctr,-4},
  {"sfmintdeg",Psfmintdeg,5},
  {"sfgcd",Psfgcd,2},
#if 0
  {"ddd",Pddd,2},
  {"newddd",Pnewddd,2},
#endif
  {"ddd_tab",Pddd_tab,2},
  {"irred_check",Pirred_check,2},
  {"nfctr_mod",Pnfctr_mod,2},
  {0,0,0},
};

/* bivariate_hensel_special(f(x,y):monic in x,g0(x),h0(y),x,y,d) */

void Pbivariate_hensel_special(NODE arg,LIST *rp)
{
  DCP dc;
  struct oVN vn[2];
  P f,g0,h0,ak,bk,gk,hk;
  V vx,vy;
  VL nvl;
  Z qk;
  Q cbd,bb;
  int d;
  NODE n;

  f = (P)ARG0(arg);
  g0 = (P)ARG1(arg);
  h0 = (P)ARG2(arg);
  vx = VR((P)ARG3(arg));
  vy = VR((P)ARG4(arg));
  d = ZTOS((Q)ARG5(arg));
  NEWVL(nvl); nvl->v = vx;
  NEWVL(NEXT(nvl)); NEXT(nvl)->v = vy;
  NEXT(NEXT(nvl)) = 0;
  vn[0].v = vy; vn[0].n = 0;  
  vn[1].v = 0; vn[1].n = 0;
  cbound(nvl,f,&cbd);
  addq(cbd,cbd,&bb);
  henzq1(g0,h0,bb,&bk,&ak,&qk);
  henmv(nvl,vn,f,g0,h0,ak,bk,(P)ONE,(P)ONE,(P)ONE,(P)ONE,qk,d,&gk,&hk);
  n = mknode(2,gk,hk);
  MKLIST(*rp,n);
}

void Pfctr(NODE arg,LIST *rp)
{
  DCP dc;

  asir_assert(ARG0(arg),O_P,"fctr");
  if ( argc(arg) == 1 )
    fctrp(CO,(P)ARG0(arg),&dc);
  else {
    asir_assert(ARG1(arg),O_P,"fctr");
    fctr_wrt_v_p(CO,(P)ARG0(arg),VR((P)ARG1(arg)),&dc);
  }
  dcptolist(dc,rp);
}

void Pgcd(NODE arg,P *rp)
{
  P p1,p2,g1,g2,g;
  Num m;
  int mod;

  p1 = (P)ARG0(arg); p2 = (P)ARG1(arg);
  asir_assert(p1,O_P,"gcd");
  asir_assert(p2,O_P,"gcd");
  if ( !p1 )
    *rp = p2;
  else if ( !p2 )
    *rp = p1;
  else if ( !qpcheck((Obj)p1) || !qpcheck((Obj)p2) )
    gcdprsp(CO,p1,p2,rp);
  else if ( argc(arg) == 2 )
    ezgcdp(CO,p1,p2,rp);
  else {
    m = (Num)ARG2(arg);
    asir_assert(m,O_P,"gcd");
    mod = ZTOS((Q)m);
    ptomp(mod,p1,&g1); ptomp(mod,p2,&g2);
    gcdprsmp(CO,mod,g1,g2,&g);
    mptop(g,rp);
  }
}

void Pgcdz(NODE arg,P *rp)
{
  P p1,p2,t;
  Q c1,c2;
  Z n1,n2,n;

  p1 = (P)ARG0(arg); p2 = (P)ARG1(arg);
  asir_assert(p1,O_P,"gcdz");
  asir_assert(p2,O_P,"gcdz");
  if ( !p1 )
    *rp = p2;
  else if ( !p2 )
    *rp = p1;
  else if ( !qpcheck((Obj)p1) || !qpcheck((Obj)p2) )
    error("gcdz : invalid argument");
  else if ( NUM(p1) || NUM(p2) ) {
    if ( NUM(p1) )
      c1 = (Q)p1;
    else
      ptozp(p1,1,&c1,&t);
    if ( NUM(p2) )
      c2 = (Q)p2;
    else
      ptozp(p2,1,&c2,&t);
    /* XXX */
    nmq(c1,&n1); nmq(c2,&n2);
    gcdz(n1,n2,&n); *rp = (P)n;
  } else {
#if 0
    w[0] = p1; w[1] = p2; nezgcdnpz(CO,w,2,rp);
#endif
    ezgcdpz(CO,p1,p2,rp);
  }
}

void Plcm(NODE arg,P *rp)
{
  P t1,t2,p1,p2,g,q;
  Q c;

  p1 = (P)ARG0(arg); p2 = (P)ARG1(arg);
  asir_assert(p1,O_P,"lcm");
  asir_assert(p2,O_P,"lcm");
  if ( !p1 || !p2 )
    *rp = 0;
  else if ( !qpcheck((Obj)p1) || !qpcheck((Obj)p2) )
    error("lcm : invalid argument");
  else {
    ptozp(p1,1,&c,&t1); ptozp(p2,1,&c,&t2);
    ezgcdp(CO,t1,t2,&g); divsp(CO,t1,g,&q); mulp(CO,q,t2,rp);
  }
}

void Psqfr(NODE arg,LIST *rp)
{
  DCP dc;

  asir_assert(ARG0(arg),O_P,"sqfr");
  sqfrp(CO,(P)ARG0(arg),&dc);
  dcptolist(dc,rp);
}

void Pufctrhint(NODE arg,LIST *rp)
{
  DCP dc;

  asir_assert(ARG0(arg),O_P,"ufctrhint");
  asir_assert(ARG1(arg),O_N,"ufctrhint");
  ufctr((P)ARG0(arg),ZTOS((Q)ARG1(arg)),&dc);
  dcptolist(dc,rp);
}

#if 0
Pmgcd(arg,rp)
NODE arg;
Obj *rp;
{
  NODE node,tn;
  int i,m;
  P *l;

  node = BDY((LIST)ARG0(arg));
  for ( i = 0, tn = node; tn; tn = NEXT(tn), i++ );
  m = i; l = (P *)ALLOCA(m*sizeof(P));
  for ( i = 0, tn = node; i < m; tn = NEXT(tn), i++ )
    l[i] = (P)BDY(tn);
  nezgcdnpz(CO,l,m,rp);
}
#endif

void Pcont(NODE arg,P *rp)
{
  DCP dc;
  int m;
  P p,p1;
  P *l;
  V v;

  asir_assert(ARG0(arg),O_P,"cont");
  p = (P)ARG0(arg);
  if ( NUM(p) )
    *rp = p;
  else {
    if ( argc(arg) == 2 ) {
      v = VR((P)ARG1(arg));
      change_mvar(CO,p,v,&p1);
      if ( VR(p1) != v ) {
        *rp = p1; return;
      } else
        p = p1;
    }
    for ( m = 0, dc = DC(p); dc; dc = NEXT(dc), m++ );
    l = (P *)ALLOCA(m*sizeof(P));
    for ( m = 0, dc = DC(p); dc; dc = NEXT(dc), m++ )
      l[m] = COEF(dc);
    nezgcdnpz(CO,l,m,rp);
  }
}

void Psfcont(NODE arg,P *rp)
{
  DCP dc;
  MP mp;
  int m;
  Obj obj;
  P p,p1;
  P *l;
  V v;

  obj = (Obj)ARG0(arg);
  if ( !obj || NUM(obj) )
    *rp = (P)obj;
  else if ( OID(obj) == O_P ) {
    p = (P)obj;
    if ( argc(arg) == 2 ) {
      v = VR((P)ARG1(arg));
      change_mvar(CO,p,v,&p1);
      if ( VR(p1) != v ) {
        *rp = p1; return;
      } else
        p = p1;
    }
    for ( m = 0, dc = DC(p); dc; dc = NEXT(dc), m++ );
    l = (P *)ALLOCA(m*sizeof(P));
    for ( m = 0, dc = DC(p); dc; dc = NEXT(dc), m++ )
      l[m] = COEF(dc);
    gcdsf(CO,l,m,rp);
  } else if ( OID(obj) == O_DP ) {
    for ( m = 0, mp = BDY((DP)obj); mp; mp = NEXT(mp), m++ );
    l = (P *)ALLOCA(m*sizeof(P));
    for ( m = 0, mp = BDY((DP)obj); mp; mp = NEXT(mp), m++)
      l[m] = (P)mp->c;
    gcdsf(CO,l,m,rp);
  }
}

void Pptozp(NODE arg,Obj *rp)
{
  Q t;
    NODE tt,p;
    NODE n,n0;
    char *key;
  P pp;
  LIST list;
    int get_factor=0;

  asir_assert(ARG0(arg),O_P,"ptozp");

    /* analyze the option */
    if ( current_option ) {
      for ( tt = current_option; tt; tt = NEXT(tt) ) {
        p = BDY((LIST)BDY(tt));
        key = BDY((STRING)BDY(p));
        /*  value = (Obj)BDY(NEXT(p)); */
        if ( !strcmp(key,"factor") )  get_factor=1;
        else {
          error("ptozp: unknown option.");
        }
      }
    }

  ptozp((P)ARG0(arg),1,&t,&pp);

    /* printexpr(NULL,t); */
  /* if the option factor is given, then it returns the answer
       in the format [zpoly, num] where num*zpoly is equal to the argument.*/
    if (get_factor) {
    n0 = mknode(2,pp,t);
      MKLIST(list,n0);
    *rp = (Obj)list;
    } else
      *rp = (Obj)pp;
}

void Pafctr(NODE arg,LIST *rp)
{
  DCP dc;
  
  asir_assert(ARG0(arg),O_P,"afctr");
  asir_assert(ARG1(arg),O_P,"afctr");
  afctr(CO,(P)ARG0(arg),(P)ARG1(arg),&dc);
  dcptolist(dc,rp);
}

void Pagcd(NODE arg,P *rp)
{
  asir_assert(ARG0(arg),O_P,"agcd");
  asir_assert(ARG1(arg),O_P,"agcd");
  asir_assert(ARG2(arg),O_P,"agcd");
  gcda(CO,(P)ARG0(arg),(P)ARG1(arg),(P)ARG2(arg),rp);
}

#if 1
#define Mulum mulum
#define Divum divum
#define Mulsum mulsum
#define Gcdum gcdum
#endif

void Mulum(), Mulsum(), Gcdum();
int Divum();

#define FCTR 0 /* berlekamp */
#define SQFR 1
#define DDD 2  /* Cantor-Zassenhauss */
#define NEWDDD 3  /* berlekamp + root-finding by Cantor-Zassenhauss */

UM *resberle();

void reduce_sfdc(DCP sfdc, DCP *dc);

void Pmodfctr(NODE arg,LIST *rp)
{
  DCP dc,dcu;
  int mod,i,t;
  P p;
  Obj u;
  VL vl;

  mod = ZTOS((Q)ARG1(arg));
  if ( mod < 0 )
    error("modfctr : invalid modulus");
  p = (P)ARG0(arg);
  clctv(CO,p,&vl);
  if ( !vl ) {
    NEWDC(dc); COEF(dc) = p; DEG(dc) = ONE; NEXT(dc) = 0;
  } else if ( !NEXT(vl) )
    modfctrp(ARG0(arg),mod,NEWDDD,&dc);
  else {
    /* XXX 16384 should be replaced by a macro */
    for ( i = 1, t = mod; t*mod < 16384; t *= mod, i++ );
    current_ff = FF_GFS;
    setmod_sf(mod,i);
    simp_ff((Obj)p,&u);
    mfctrsf(CO,(P)u,&dcu);
    reduce_sfdc(dcu,&dc);
  }
  if ( !dc ) {
    NEWDC(dc); COEF(dc) = 0; DEG(dc) = ONE; NEXT(dc) = 0;
  }
  dcptolist(dc,rp);
}

void Psfgcd(NODE arg,P *rp)
{
  P ps[2];

  ps[0] = (P)ARG0(arg);
  ps[1] = (P)ARG1(arg);
  gcdsf(CO,ps,2,rp);
}

void Psffctr(NODE arg,LIST *rp)
{
  DCP dc;

  mfctrsf(CO,ARG0(arg),&dc);
  dcptolist(dc,rp);
}

void Psfsqfr(NODE arg,LIST *rp)
{
  DCP dc;

  sqfrsf(CO,ARG0(arg),&dc);
  dcptolist(dc,rp);
}

void Psfufctr(NODE arg,LIST *rp)
{
  DCP dc;

  ufctrsf(ARG0(arg),&dc);
  dcptolist(dc,rp);
}

void Psfbfctr(NODE arg,LIST *rp)
{
  V x,y;
  DCP dc,dct;
  P t;
  struct oVL vl1,vl2;
  VL vl;
  int degbound;

  x = VR((P)ARG1(arg));
  y = VR((P)ARG2(arg));
  vl1.v = x; vl1.next = &vl2;
  vl2.v = y; vl2.next = 0;
  vl = &vl1;
  if ( argc(arg) == 4 )
    degbound = ZTOS((Q)ARG3(arg));
  else
    degbound = -1;

  sfbfctr((P)ARG0(arg),x,y,degbound,&dc);
  for ( dct = dc; dct; dct = NEXT(dct) ) {
    reorderp(CO,vl,COEF(dct),&t); COEF(dct) = t;
  }
  dcptolist(dc,rp);
}

void Psfmintdeg(NODE arg,P *rp)
{
  V x,y;
  P r;
  struct oVL vl1,vl2;
  VL vl;
  int dy,c;

  x = VR((P)ARG1(arg));
  y = VR((P)ARG2(arg));
  vl1.v = x; vl1.next = &vl2;
  vl2.v = y; vl2.next = 0;
  vl = &vl1;
  dy = ZTOS((Q)ARG3(arg));
  c = ZTOS((Q)ARG4(arg));
  sfmintdeg(vl,(P)ARG0(arg),dy,c,&r);
  reorderp(CO,vl,r,rp);
}

void Pmodsqfr(NODE arg,LIST *rp)
{
  DCP dc;

  if ( !ARG0(arg) ) {
    NEWDC(dc); COEF(dc) = 0; DEG(dc) = ONE; NEXT(dc) = 0;
  } else
    modfctrp(ARG0(arg),ZTOS((Q)ARG1(arg)),SQFR,&dc);
  dcptolist(dc,rp);
}

void Pddd(NODE arg,LIST *rp)
{
  DCP dc;

  if ( !ARG0(arg) ) {
    NEWDC(dc); COEF(dc) = 0; DEG(dc) = ONE; NEXT(dc) = 0;
  } else
    modfctrp(ARG0(arg),ZTOS((Q)ARG1(arg)),DDD,&dc);
  dcptolist(dc,rp);
}

void Pnewddd(NODE arg,LIST *rp)
{
  DCP dc=0;

  if ( !ARG0(arg) ) {
    NEWDC(dc); COEF(dc) = 0; DEG(dc) = ONE; NEXT(dc) = 0;
  } else 
    modfctrp(ARG0(arg),ZTOS((Q)ARG1(arg)),NEWDDD,&dc);
  dcptolist(dc,rp);
}

void Pirred_check(NODE arg,Z *rp)
{
  P p;
  UM mp;
  int r,mod;

  p = (P)ARG0(arg);
  if ( !p ) {
    *rp = 0; return;
  }
  mp = W_UMALLOC(UDEG(p));
  mod = ZTOS((Q)ARG1(arg));
  ptoum(mod,p,mp);
  r = irred_check(mp,mod);
  if ( r )
    *rp = ONE;
  else
    *rp = 0;
}

void Pnfctr_mod(NODE arg,Z *rp)
{
  P p;
  UM mp;
  int r,mod;

  p = (P)ARG0(arg);
  if ( !p ) {
    *rp = 0; return;
  }
  mp = W_UMALLOC(UDEG(p));
  mod = ZTOS((Q)ARG1(arg));
  ptoum(mod,p,mp);
  r = nfctr_mod(mp,mod);
  STOZ(r,*rp);
}

void Pddd_tab(NODE arg,VECT *rp)
{
  P p;
  UM mp,t,q,r1,w,w1;
  UM *r,*s;
  int dr,mod,n,i;
  VECT result;
  V v;

  p = (P)ARG0(arg); mod = ZTOS((Q)ARG1(arg));
  v = VR(p);
  n = UDEG(p); mp = W_UMALLOC(n);
  ptoum(mod,p,mp);
  r = (UM *)W_ALLOC(n); s = (UM *)W_ALLOC(n);
  r[0] = UMALLOC(0); DEG(r[0]) = 0; COEF(r[0])[0] = 1;
  t = W_UMALLOC(mod); bzero(COEF(t),sizeof(int)*(mod+1));
  DEG(t) = mod; COEF(t)[mod] = 1;
  q = W_UMALLOC(mod);
  dr = divum(mod,t,mp,q);
  DEG(t) = dr; r[1] = r1 = UMALLOC(dr); cpyum(t,r1);
  s[0] = W_UMALLOC(dr); cpyum(t,s[0]);
  w = W_UMALLOC(n); bzero(COEF(w),sizeof(int)*(n+1));
  w1 = W_UMALLOC(2*n); bzero(COEF(w1),sizeof(int)*(2*n+1));
  for ( i = 1; i < n; i++ ) {
    DEG(w) = i; COEF(w)[i-1] = 0; COEF(w)[i] = 1;
    mulum(mod,r1,w,w1);
    dr = divum(mod,w1,mp,q); DEG(w1) = dr;
    s[i] = W_UMALLOC(dr); cpyum(w1,s[i]);
  }
  for ( i = 2; i < n; i++ ) {
    mult_mod_tab(r[i-1],mod,s,w,n);
    r[i] = UMALLOC(DEG(w)); cpyum(w,r[i]);
  }
  MKVECT(result,n);
  for ( i = 0; i < n; i++ )
    umtop(v,r[i],(P *)&BDY(result)[i]);
  *rp = result;
}

void reduce_sfdc(DCP sfdc,DCP *dcr)
{
  P c,t,s,u,f;
  DCP dc0,dc,tdc;
  DCP *a;
  int i,j,n;

  if ( !current_gfs_ext ) {
    /* we simply apply sfptop() */
    for ( dc0 = 0; sfdc; sfdc = NEXT(sfdc) ) {
      NEXTDC(dc0,dc);
      DEG(dc) = DEG(sfdc);
      sfptop(COEF(sfdc),&COEF(dc));
    }
    NEXT(dc) = 0;
    *dcr = dc0;
    return;
  }

  if ( NUM(COEF(sfdc)) ) {
    sfptop(COEF(sfdc),&c);
    sfdc = NEXT(sfdc);
  } else
    c = (P)ONE;

  for ( n = 0, tdc = sfdc; tdc; tdc = NEXT(tdc), n++ );
  a = (DCP *)ALLOCA(n*sizeof(DCP));
  for ( i = 0, tdc = sfdc; i < n; tdc = NEXT(tdc), i++ )
    a[i] = tdc;

  dc0 = 0; NEXTDC(dc0,dc); DEG(dc) = ONE; COEF(dc) = c;
  for ( i = 0; i < n; i++ ) {
    if ( !a[i] )
      continue;
    t = COEF(a[i]);
    f = t;
    while ( 1 ) {
      sf_galois_action(t,(Q)ONE,&s);
      for ( j = i; j < n; j++ )
        if ( a[j] && !compp(CO,s,COEF(a[j])) )
          break;
      if ( j == n )
        error("reduce_sfdc : cannot happen");
      if ( j == i ) {
        NEXTDC(dc0,dc); DEG(dc) = DEG(a[i]);
        sfptop(f,&COEF(dc));
        break;
      } else {
        mulp(CO,f,s,&u); f = u;
        t = s;
        a[j] = 0;
      }
    }
  }
  *dcr = dc0;
}