[BACK]Return to int.c CVS log [TXT][DIR] Up to [local] / OpenXM_contrib2 / asir2000 / builtin

Diff for /OpenXM_contrib2/asir2000/builtin/int.c between version 1.13 and 1.14

version 1.13, 2015/08/14 13:51:54 version 1.14, 2018/03/29 01:32:50
Line 45 
Line 45 
  * DEVELOPER SHALL HAVE NO LIABILITY IN CONNECTION WITH THE USE,   * DEVELOPER SHALL HAVE NO LIABILITY IN CONNECTION WITH THE USE,
  * PERFORMANCE OR NON-PERFORMANCE OF THE SOFTWARE.   * PERFORMANCE OR NON-PERFORMANCE OF THE SOFTWARE.
  *   *
  * $OpenXM: OpenXM_contrib2/asir2000/builtin/int.c,v 1.12 2015/08/06 10:01:52 fujimoto Exp $   * $OpenXM: OpenXM_contrib2/asir2000/builtin/int.c,v 1.13 2015/08/14 13:51:54 fujimoto Exp $
 */  */
 #include "ca.h"  #include "ca.h"
 #include "parse.h"  #include "parse.h"
Line 71  void Pimaxrsh(), Pilen();
Line 71  void Pimaxrsh(), Pilen();
 void Ptype_t_NB();  void Ptype_t_NB();
   
 struct ftab int_tab[] = {  struct ftab int_tab[] = {
         {"dp_set_mpi",Pdp_set_mpi,-1},    {"dp_set_mpi",Pdp_set_mpi,-1},
         {"isqrt",Pisqrt,1},    {"isqrt",Pisqrt,1},
         {"idiv",Pidiv,2},    {"idiv",Pidiv,2},
         {"irem",Pirem,2},    {"irem",Pirem,2},
         {"iqr",Piqr,2},    {"iqr",Piqr,2},
         {"igcd",Pigcd,-2},    {"igcd",Pigcd,-2},
         {"ilcm",Pilcm,2},    {"ilcm",Pilcm,2},
         {"up2_inv",Pup2_inv,2},    {"up2_inv",Pup2_inv,2},
         {"up2_init_eg",Pup2_init_eg,0},    {"up2_init_eg",Pup2_init_eg,0},
         {"up2_show_eg",Pup2_show_eg,0},    {"up2_show_eg",Pup2_show_eg,0},
         {"type_t_NB",Ptype_t_NB,2},    {"type_t_NB",Ptype_t_NB,2},
         {"gf2nton",Pgf2nton,1},    {"gf2nton",Pgf2nton,1},
         {"ntogf2n",Pntogf2n,1},    {"ntogf2n",Pntogf2n,1},
         {"set_upkara",Pset_upkara,-1},    {"set_upkara",Pset_upkara,-1},
         {"set_uptkara",Pset_uptkara,-1},    {"set_uptkara",Pset_uptkara,-1},
         {"set_up2kara",Pset_up2kara,-1},    {"set_up2kara",Pset_up2kara,-1},
         {"set_upfft",Pset_upfft,-1},    {"set_upfft",Pset_upfft,-1},
         {"inv",Pinv,2},    {"inv",Pinv,2},
         {"inttorat",Pinttorat,3},    {"inttorat",Pinttorat,3},
         {"fac",Pfac,1},    {"fac",Pfac,1},
         {"prime",Pprime,1},    {"prime",Pprime,1},
         {"lprime",Plprime,1},    {"lprime",Plprime,1},
         {"random",Prandom,-1},    {"random",Prandom,-1},
         {"lrandom",Plrandom,1},    {"lrandom",Plrandom,1},
         {"iand",Piand,2},    {"iand",Piand,2},
         {"ior",Pior,2},    {"ior",Pior,2},
         {"ixor",Pixor,2},    {"ixor",Pixor,2},
         {"ishift",Pishift,2},    {"ishift",Pishift,2},
         {"small_jacobi",Psmall_jacobi,2},    {"small_jacobi",Psmall_jacobi,2},
   
         {"igcdbin",Pigcdbin,2},         /* HM@CCUT extension */    {"igcdbin",Pigcdbin,2},    /* HM@CCUT extension */
         {"igcdbmod",Pigcdbmod,2},       /* HM@CCUT extension */    {"igcdbmod",Pigcdbmod,2},  /* HM@CCUT extension */
         {"igcdeuc",PigcdEuc,2},         /* HM@CCUT extension */    {"igcdeuc",PigcdEuc,2},    /* HM@CCUT extension */
         {"igcdacc",Pigcdacc,2},         /* HM@CCUT extension */    {"igcdacc",Pigcdacc,2},    /* HM@CCUT extension */
         {"igcdcntl",Pigcdcntl,-1},      /* HM@CCUT extension */    {"igcdcntl",Pigcdcntl,-1},  /* HM@CCUT extension */
   
         {"mt_save",Pmt_save,1},    {"mt_save",Pmt_save,1},
         {"mt_load",Pmt_load,1},    {"mt_load",Pmt_load,1},
         {"ntoint32",Pntoint32,1},    {"ntoint32",Pntoint32,1},
         {"int32ton",Pint32ton,1},    {"int32ton",Pint32ton,1},
         {0,0,0},    {0,0,0},
 };  };
   
 static int is_prime_small(unsigned int);  static int is_prime_small(unsigned int);
Line 121  int mpi_mag;
Line 121  int mpi_mag;
   
 void Pntoint32(NODE arg,USINT *rp)  void Pntoint32(NODE arg,USINT *rp)
 {  {
         Q q;    Q q;
         unsigned int t;    unsigned int t;
   
         asir_assert(ARG0(arg),O_N,"ntoint32");    asir_assert(ARG0(arg),O_N,"ntoint32");
         q = (Q)ARG0(arg);    q = (Q)ARG0(arg);
         if ( !q ) {    if ( !q ) {
                 MKUSINT(*rp,0);      MKUSINT(*rp,0);
                 return;      return;
         }    }
         if ( NID(q)!=N_Q || !INT(q) || PL(NM(q))>1 )    if ( NID(q)!=N_Q || !INT(q) || PL(NM(q))>1 )
                 error("ntoint32 : invalid argument");      error("ntoint32 : invalid argument");
         t = BD(NM(q))[0];    t = BD(NM(q))[0];
         if ( SGN(q) < 0 )    if ( SGN(q) < 0 )
                 t = -(int)t;      t = -(int)t;
         MKUSINT(*rp,t);    MKUSINT(*rp,t);
 }  }
   
 void Pint32ton(NODE arg,Q *rp)  void Pint32ton(NODE arg,Q *rp)
 {  {
         int t;    int t;
   
         asir_assert(ARG0(arg),O_USINT,"int32ton");    asir_assert(ARG0(arg),O_USINT,"int32ton");
         t = (int)BDY((USINT)ARG0(arg));    t = (int)BDY((USINT)ARG0(arg));
         STOQ(t,*rp);    STOQ(t,*rp);
 }  }
   
 void Pdp_set_mpi(NODE arg,Q *rp)  void Pdp_set_mpi(NODE arg,Q *rp)
 {  {
         if ( arg ) {    if ( arg ) {
                 asir_assert(ARG0(arg),O_N,"dp_set_mpi");      asir_assert(ARG0(arg),O_N,"dp_set_mpi");
                 mpi_mag = QTOS((Q)ARG0(arg));      mpi_mag = QTOS((Q)ARG0(arg));
         }    }
         STOQ(mpi_mag,*rp);    STOQ(mpi_mag,*rp);
 }  }
   
 void Psmall_jacobi(NODE arg,Q *rp)  void Psmall_jacobi(NODE arg,Q *rp)
 {  {
         Q a,m;    Q a,m;
         int a0,m0,s;    int a0,m0,s;
   
         a = (Q)ARG0(arg);    a = (Q)ARG0(arg);
         m = (Q)ARG1(arg);    m = (Q)ARG1(arg);
         asir_assert(a,O_N,"small_jacobi");    asir_assert(a,O_N,"small_jacobi");
         asir_assert(m,O_N,"small_jacobi");    asir_assert(m,O_N,"small_jacobi");
         if ( !a )    if ( !a )
                  *rp = ONE;       *rp = ONE;
         else if ( !m || !INT(m) || !INT(a)    else if ( !m || !INT(m) || !INT(a)
                 || PL(NM(m))>1 || PL(NM(a))>1 || SGN(m) < 0 || EVENN(NM(m)) )      || PL(NM(m))>1 || PL(NM(a))>1 || SGN(m) < 0 || EVENN(NM(m)) )
                 error("small_jacobi : invalid input");      error("small_jacobi : invalid input");
         else {    else {
                 a0 = QTOS(a); m0 = QTOS(m);      a0 = QTOS(a); m0 = QTOS(m);
                 s = small_jacobi(a0,m0);      s = small_jacobi(a0,m0);
                 STOQ(s,*rp);      STOQ(s,*rp);
         }    }
 }  }
   
 int small_jacobi(int a,int m)  int small_jacobi(int a,int m)
 {  {
         int m4,m8,a4,j1,i,s;    int m4,m8,a4,j1,i,s;
   
         a %= m;    a %= m;
         if ( a == 0 || a == 1 )    if ( a == 0 || a == 1 )
                 return 1;      return 1;
         else if ( a < 0 ) {    else if ( a < 0 ) {
                 j1 = small_jacobi(-a,m);      j1 = small_jacobi(-a,m);
                 m4 = m%4;      m4 = m%4;
                 return m4==1?j1:-j1;      return m4==1?j1:-j1;
         } else {    } else {
                 for ( i = 0; a && !(a&1); i++, a >>= 1 );      for ( i = 0; a && !(a&1); i++, a >>= 1 );
                 if ( i&1 ) {      if ( i&1 ) {
                         m8 = m%8;        m8 = m%8;
                         s = (m8==1||m8==7) ? 1 : -1;        s = (m8==1||m8==7) ? 1 : -1;
                 } else      } else
                         s = 1;        s = 1;
                 /* a, m are odd */      /* a, m are odd */
                 j1 = small_jacobi(m%a,a);      j1 = small_jacobi(m%a,a);
                 m4 = m%4; a4 = a%4;      m4 = m%4; a4 = a%4;
                 s *= (m4==1||a4==1) ? 1 : -1;      s *= (m4==1||a4==1) ? 1 : -1;
                 return j1*s;      return j1*s;
         }    }
 }  }
   
 void Ptype_t_NB(NODE arg,Q *rp)  void Ptype_t_NB(NODE arg,Q *rp)
 {  {
         if ( TypeT_NB_check(QTOS((Q)ARG0(arg)),QTOS((Q)ARG1(arg))))    if ( TypeT_NB_check(QTOS((Q)ARG0(arg)),QTOS((Q)ARG1(arg))))
                 *rp = ONE;      *rp = ONE;
         else    else
                 *rp = 0;      *rp = 0;
 }  }
   
 int TypeT_NB_check(unsigned int m, unsigned int t)  int TypeT_NB_check(unsigned int m, unsigned int t)
 {  {
         unsigned int p,k,u,h,d;    unsigned int p,k,u,h,d;
   
         if ( !(m%8) )    if ( !(m%8) )
                 return 0;      return 0;
         p = t*m+1;    p = t*m+1;
         if ( !is_prime_small(p) )    if ( !is_prime_small(p) )
                 return 0;      return 0;
         for ( k = 1, u = 2%p; ; k++ )    for ( k = 1, u = 2%p; ; k++ )
                 if ( u == 1 )      if ( u == 1 )
                         break;        break;
                 else      else
                         u = (2*u)%p;        u = (2*u)%p;
         h = t*m/k;    h = t*m/k;
         d = gcd_small(h,m);    d = gcd_small(h,m);
         return d == 1 ? 1 : 0;    return d == 1 ? 1 : 0;
 }  }
   
 /*  /*
Line 238  int TypeT_NB_check(unsigned int m, unsigned int t)
Line 238  int TypeT_NB_check(unsigned int m, unsigned int t)
   
 static int is_prime_small(unsigned int a)  static int is_prime_small(unsigned int a)
 {  {
         unsigned int b,t,i;    unsigned int b,t,i;
   
         if ( !(a%2) ) return 0;    if ( !(a%2) ) return 0;
         for ( t = a, i = 0; t; i++, t >>= 1 );    for ( t = a, i = 0; t; i++, t >>= 1 );
         /* b >= sqrt(a) */    /* b >= sqrt(a) */
         b = 1<<((i+1)/2);    b = 1<<((i+1)/2);
   
         /* divisibility test by all odd numbers <= b */    /* divisibility test by all odd numbers <= b */
         for ( i = 3; i <= b; i += 2 )    for ( i = 3; i <= b; i += 2 )
                 if ( !(a%i) )      if ( !(a%i) )
                         return 0;        return 0;
         return 1;    return 1;
 }  }
   
 /*  /*
Line 261  static int is_prime_small(unsigned int a)
Line 261  static int is_prime_small(unsigned int a)
   
 static unsigned int gcd_small(unsigned int a,unsigned int b)  static unsigned int gcd_small(unsigned int a,unsigned int b)
 {  {
         unsigned int t;    unsigned int t;
   
         if ( b > a ) {    if ( b > a ) {
                 t = a; a = b; b = t;      t = a; a = b; b = t;
         }    }
         /* Euclid's algorithm */    /* Euclid's algorithm */
         while ( 1 )    while ( 1 )
                 if ( !(t = a%b) ) return b;      if ( !(t = a%b) ) return b;
                 else {      else {
                         a = b; b = t;        a = b; b = t;
                 }      }
 }  }
   
 void Pmt_save(NODE arg,Q *rp)  void Pmt_save(NODE arg,Q *rp)
 {  {
         int ret;    int ret;
   
         ret = mt_save(BDY((STRING)ARG0(arg)));    ret = mt_save(BDY((STRING)ARG0(arg)));
         STOQ(ret,*rp);    STOQ(ret,*rp);
 }  }
   
 void Pmt_load(NODE arg,Q *rp)  void Pmt_load(NODE arg,Q *rp)
 {  {
         int ret;    int ret;
   
         ret = mt_load(BDY((STRING)ARG0(arg)));    ret = mt_load(BDY((STRING)ARG0(arg)));
         STOQ(ret,*rp);    STOQ(ret,*rp);
 }  }
   
 void Pisqrt(NODE arg,Q *rp)  void Pisqrt(NODE arg,Q *rp)
 {  {
         Q a;    Q a;
         N r;    N r;
   
         a = (Q)ARG0(arg);    a = (Q)ARG0(arg);
         asir_assert(a,O_N,"isqrt");    asir_assert(a,O_N,"isqrt");
         if ( !a )    if ( !a )
                 *rp = 0;      *rp = 0;
         else if ( SGN(a) < 0 )    else if ( SGN(a) < 0 )
                 error("isqrt : negative argument");      error("isqrt : negative argument");
         else {    else {
                 isqrt(NM(a),&r);      isqrt(NM(a),&r);
                 NTOQ(r,1,*rp);      NTOQ(r,1,*rp);
         }    }
 }  }
   
 void Pidiv(NODE arg,Obj *rp)  void Pidiv(NODE arg,Obj *rp)
 {  {
         N q,r;    N q,r;
         Q a;    Q a;
         Q dnd,dvr;    Q dnd,dvr;
   
         dnd = (Q)ARG0(arg); dvr = (Q)ARG1(arg);    dnd = (Q)ARG0(arg); dvr = (Q)ARG1(arg);
         asir_assert(dnd,O_N,"idiv");    asir_assert(dnd,O_N,"idiv");
         asir_assert(dvr,O_N,"idiv");    asir_assert(dvr,O_N,"idiv");
         if ( !dvr )    if ( !dvr )
                 error("idiv: division by 0");      error("idiv: division by 0");
         else if ( !dnd )    else if ( !dnd )
                 *rp = 0;      *rp = 0;
         else {    else {
                 divn(NM(dnd),NM(dvr),&q,&r);      divn(NM(dnd),NM(dvr),&q,&r);
                 NTOQ(q,SGN(dnd)*SGN(dvr),a); *rp = (Obj)a;      NTOQ(q,SGN(dnd)*SGN(dvr),a); *rp = (Obj)a;
         }    }
 }  }
   
 void Pirem(NODE arg,Obj *rp)  void Pirem(NODE arg,Obj *rp)
 {  {
         N q,r;    N q,r;
         Q a;    Q a;
         Q dnd,dvr;    Q dnd,dvr;
   
         dnd = (Q)ARG0(arg); dvr = (Q)ARG1(arg);    dnd = (Q)ARG0(arg); dvr = (Q)ARG1(arg);
         asir_assert(dnd,O_N,"irem");    asir_assert(dnd,O_N,"irem");
         asir_assert(dvr,O_N,"irem");    asir_assert(dvr,O_N,"irem");
         if ( !dvr )    if ( !dvr )
                 error("irem: division by 0");      error("irem: division by 0");
         else if ( !dnd )    else if ( !dnd )
                 *rp = 0;      *rp = 0;
         else {    else {
                 divn(NM(dnd),NM(dvr),&q,&r);      divn(NM(dnd),NM(dvr),&q,&r);
                 NTOQ(r,SGN(dnd),a); *rp = (Obj)a;      NTOQ(r,SGN(dnd),a); *rp = (Obj)a;
         }    }
 }  }
   
 void Piqr(NODE arg,LIST *rp)  void Piqr(NODE arg,LIST *rp)
 {  {
         N q,r;    N q,r;
         Q a,b;    Q a,b;
         Q dnd,dvr;    Q dnd,dvr;
         NODE node1,node2;    NODE node1,node2;
   
         dnd = (Q)ARG0(arg); dvr = (Q)ARG1(arg);    dnd = (Q)ARG0(arg); dvr = (Q)ARG1(arg);
         if ( !dvr )    if ( !dvr )
                 error("iqr: division by 0");      error("iqr: division by 0");
         else if ( !dnd )    else if ( !dnd )
                 a = b = 0;      a = b = 0;
         else if ( OID(dnd) == O_VECT ) {    else if ( OID(dnd) == O_VECT ) {
                 iqrv((VECT)dnd,dvr,rp); return;      iqrv((VECT)dnd,dvr,rp); return;
         } else {    } else {
                 asir_assert(dnd,O_N,"iqr");      asir_assert(dnd,O_N,"iqr");
                 asir_assert(dvr,O_N,"iqr");      asir_assert(dvr,O_N,"iqr");
                 divn(NM(dnd),NM(dvr),&q,&r);      divn(NM(dnd),NM(dvr),&q,&r);
                 NTOQ(q,SGN(dnd)*SGN(dvr),a);      NTOQ(q,SGN(dnd)*SGN(dvr),a);
                 NTOQ(r,SGN(dnd),b);      NTOQ(r,SGN(dnd),b);
         }    }
         MKNODE(node2,b,0); MKNODE(node1,a,node2); MKLIST(*rp,node1);    MKNODE(node2,b,0); MKNODE(node1,a,node2); MKLIST(*rp,node1);
 }  }
   
 void Pinttorat(NODE arg,LIST *rp)  void Pinttorat(NODE arg,LIST *rp)
 {  {
         Q cq,qq,t,u1,v1,r1,nm;    Q cq,qq,t,u1,v1,r1,nm;
         N m,b,q,r,c,u2,v2,r2;    N m,b,q,r,c,u2,v2,r2;
         NODE node1,node2;    NODE node1,node2;
         P p;    P p;
   
         asir_assert(ARG0(arg),O_N,"inttorat");    asir_assert(ARG0(arg),O_N,"inttorat");
         asir_assert(ARG1(arg),O_N,"inttorat");    asir_assert(ARG1(arg),O_N,"inttorat");
         asir_assert(ARG2(arg),O_N,"inttorat");    asir_assert(ARG2(arg),O_N,"inttorat");
         cq = (Q)ARG0(arg); m = NM((Q)ARG1(arg)); b = NM((Q)ARG2(arg));    cq = (Q)ARG0(arg); m = NM((Q)ARG1(arg)); b = NM((Q)ARG2(arg));
         if ( !cq ) {    if ( !cq ) {
                 MKNODE(node2,(pointer)ONE,0); MKNODE(node1,0,node2); MKLIST(*rp,node1);      MKNODE(node2,(pointer)ONE,0); MKNODE(node1,0,node2); MKLIST(*rp,node1);
                 return;      return;
         }    }
         divn(NM(cq),m,&q,&r);    divn(NM(cq),m,&q,&r);
         if ( !r ) {    if ( !r ) {
                 MKNODE(node2,(pointer)ONE,0); MKNODE(node1,0,node2); MKLIST(*rp,node1);      MKNODE(node2,(pointer)ONE,0); MKNODE(node1,0,node2); MKLIST(*rp,node1);
                 return;      return;
         } else if ( SGN(cq) < 0 ) {    } else if ( SGN(cq) < 0 ) {
                 subn(m,r,&c);      subn(m,r,&c);
         } else    } else
                 c = r;      c = r;
         u1 = 0; v1 = ONE; u2 = m; v2 = c;    u1 = 0; v1 = ONE; u2 = m; v2 = c;
         while ( cmpn(v2,b) >= 0 ) {    while ( cmpn(v2,b) >= 0 ) {
                 divn(u2,v2,&q,&r2); u2 = v2; v2 = r2;      divn(u2,v2,&q,&r2); u2 = v2; v2 = r2;
                 NTOQ(q,1,qq); mulq(qq,v1,&t); subq(u1,t,&r1); u1 = v1; v1 = r1;      NTOQ(q,1,qq); mulq(qq,v1,&t); subq(u1,t,&r1); u1 = v1; v1 = r1;
         }    }
         if ( cmpn(NM(v1),b) >= 0 )    if ( cmpn(NM(v1),b) >= 0 )
                 *rp = 0;      *rp = 0;
         else {    else {
                 if ( SGN(v1) < 0 ) {      if ( SGN(v1) < 0 ) {
                         chsgnp((P)v1,&p); v1 = (Q)p; NTOQ(v2,-1,nm);        chsgnp((P)v1,&p); v1 = (Q)p; NTOQ(v2,-1,nm);
                 } else      } else
                         NTOQ(v2,1,nm);        NTOQ(v2,1,nm);
                 MKNODE(node2,v1,0); MKNODE(node1,nm,node2); MKLIST(*rp,node1);      MKNODE(node2,v1,0); MKNODE(node1,nm,node2); MKLIST(*rp,node1);
         }    }
 }  }
   
 void Pigcd(NODE arg,Q *rp)  void Pigcd(NODE arg,Q *rp)
 {  {
         N g;    N g;
         Q n1,n2,a;    Q n1,n2,a;
   
         if ( argc(arg) == 1 ) {    if ( argc(arg) == 1 ) {
                 igcdv((VECT)ARG0(arg),rp);      igcdv((VECT)ARG0(arg),rp);
                 return;      return;
         }    }
         n1 = (Q)ARG0(arg); n2 = (Q)ARG1(arg);    n1 = (Q)ARG0(arg); n2 = (Q)ARG1(arg);
         asir_assert(n1,O_N,"igcd");    asir_assert(n1,O_N,"igcd");
         asir_assert(n2,O_N,"igcd");    asir_assert(n2,O_N,"igcd");
         if ( !n1 )    if ( !n1 )
                 *rp = n2;      *rp = n2;
         else if ( !n2 )    else if ( !n2 )
                 *rp = n1;      *rp = n1;
         else {    else {
                 gcdn(NM(n1),NM(n2),&g);      gcdn(NM(n1),NM(n2),&g);
                 NTOQ(g,1,a); *rp = a;      NTOQ(g,1,a); *rp = a;
         }    }
 }  }
   
 int comp_n(N *a,N *b)  int comp_n(N *a,N *b)
 {  {
         return cmpn(*a,*b);    return cmpn(*a,*b);
 }  }
   
 void iqrv(VECT a,Q dvr,LIST *rp)  void iqrv(VECT a,Q dvr,LIST *rp)
 {  {
         int i,n;    int i,n;
         VECT q,r;    VECT q,r;
         Q dnd,qi,ri;    Q dnd,qi,ri;
         Q *b;    Q *b;
         N qn,rn;    N qn,rn;
         NODE n0,n1;    NODE n0,n1;
   
         if ( !dvr )    if ( !dvr )
                 error("iqrv: division by 0");      error("iqrv: division by 0");
         n = a->len; b = (Q *)BDY(a);    n = a->len; b = (Q *)BDY(a);
         MKVECT(q,n); MKVECT(r,n);    MKVECT(q,n); MKVECT(r,n);
         for ( i = 0; i < n; i++ ) {    for ( i = 0; i < n; i++ ) {
                 dnd = b[i];      dnd = b[i];
                 if ( !dnd ) {      if ( !dnd ) {
                         qi = ri = 0;        qi = ri = 0;
                 } else {      } else {
                         divn(NM(dnd),NM(dvr),&qn,&rn);        divn(NM(dnd),NM(dvr),&qn,&rn);
                         NTOQ(qn,SGN(dnd)*SGN(dvr),qi);        NTOQ(qn,SGN(dnd)*SGN(dvr),qi);
                         NTOQ(rn,SGN(dnd),ri);        NTOQ(rn,SGN(dnd),ri);
                 }      }
                 BDY(q)[i] = (pointer)qi; BDY(r)[i] = (pointer)ri;      BDY(q)[i] = (pointer)qi; BDY(r)[i] = (pointer)ri;
         }    }
         MKNODE(n1,r,0); MKNODE(n0,q,n1); MKLIST(*rp,n0);    MKNODE(n1,r,0); MKNODE(n0,q,n1); MKLIST(*rp,n0);
 }  }
   
 /*  /*
Line 468  void iqrv(VECT a,Q dvr,LIST *rp)
Line 468  void iqrv(VECT a,Q dvr,LIST *rp)
   
 void igcd_cofactor(Q a,Q b,Q *gcd,Q *ca,Q *cb)  void igcd_cofactor(Q a,Q b,Q *gcd,Q *ca,Q *cb)
 {  {
         N gn,tn;    N gn,tn;
   
         if ( !a ) {    if ( !a ) {
                 if ( !b )      if ( !b )
                         error("igcd_cofactor : invalid input");        error("igcd_cofactor : invalid input");
                 else {      else {
                         *ca = 0;        *ca = 0;
                         *cb = ONE;        *cb = ONE;
                         *gcd = b;        *gcd = b;
                 }      }
         } else if ( !b ) {    } else if ( !b ) {
                 *ca = ONE;      *ca = ONE;
                 *cb = 0;      *cb = 0;
                 *gcd = a;      *gcd = a;
         } else {    } else {
                 gcdn(NM(a),NM(b),&gn); NTOQ(gn,1,*gcd);      gcdn(NM(a),NM(b),&gn); NTOQ(gn,1,*gcd);
                 divsn(NM(a),gn,&tn); NTOQ(tn,SGN(a),*ca);      divsn(NM(a),gn,&tn); NTOQ(tn,SGN(a),*ca);
                 divsn(NM(b),gn,&tn); NTOQ(tn,SGN(b),*cb);      divsn(NM(b),gn,&tn); NTOQ(tn,SGN(b),*cb);
         }    }
 }  }
   
 void igcdv(VECT a,Q *rp)  void igcdv(VECT a,Q *rp)
 {  {
         int i,j,n,nz;    int i,j,n,nz;
         N g,gt,q,r;    N g,gt,q,r;
         N *c;    N *c;
   
         n = a->len;    n = a->len;
         c = (N *)ALLOCA(n*sizeof(N));    c = (N *)ALLOCA(n*sizeof(N));
         for ( i = 0; i < n; i++ )    for ( i = 0; i < n; i++ )
                 c[i] = BDY(a)[i]?NM((Q)(BDY(a)[i])):0;      c[i] = BDY(a)[i]?NM((Q)(BDY(a)[i])):0;
         qsort(c,n,sizeof(N),(int (*) (const void *,const void *))comp_n);    qsort(c,n,sizeof(N),(int (*) (const void *,const void *))comp_n);
         for ( ; n && ! *c; n--, c++ );    for ( ; n && ! *c; n--, c++ );
   
         if ( !n ) {    if ( !n ) {
                 *rp = 0; return;      *rp = 0; return;
         } else if ( n == 1 ) {    } else if ( n == 1 ) {
                 NTOQ(*c,1,*rp); return;      NTOQ(*c,1,*rp); return;
         }    }
         gcdn(c[0],c[1],&g);    gcdn(c[0],c[1],&g);
 #if 0  #if 0
         for ( i = 2; i < n; i++ ) {    for ( i = 2; i < n; i++ ) {
                 divn(c[i],g,&q,&r);      divn(c[i],g,&q,&r);
                 gcdn(g,r,&gt);      gcdn(g,r,&gt);
                 if ( !cmpn(g,gt) ) {      if ( !cmpn(g,gt) ) {
                         for ( j = i+1, nz = 0; j < n; j++ ) {        for ( j = i+1, nz = 0; j < n; j++ ) {
                                 divn(c[j],g,&q,&r); c[j] = r;          divn(c[j],g,&q,&r); c[j] = r;
                                 if ( r )          if ( r )
                                         nz = 1;            nz = 1;
                         }        }
                 } else      } else
                         g = gt;        g = gt;
         }    }
 #else  #else
         for ( i = 2; i < n; i++ ) {    for ( i = 2; i < n; i++ ) {
                 gcdn(g,c[i],&gt); g = gt;      gcdn(g,c[i],&gt); g = gt;
         }    }
 #endif  #endif
         NTOQ(g,1,*rp);    NTOQ(g,1,*rp);
 }  }
   
 void igcdv_estimate(VECT a,Q *rp)  void igcdv_estimate(VECT a,Q *rp)
 {  {
         int n,i,m;    int n,i,m;
         N s,t,u,g;    N s,t,u,g;
         Q *q;    Q *q;
   
         n = a->len; q = (Q *)a->body;    n = a->len; q = (Q *)a->body;
         if ( n == 1 ) {    if ( n == 1 ) {
                 NTOQ(NM(q[0]),1,*rp); return;      NTOQ(NM(q[0]),1,*rp); return;
         }    }
   
         m = n/2;    m = n/2;
         for ( i = 0 , s = 0; i < m; i++ ) {    for ( i = 0 , s = 0; i < m; i++ ) {
                 addn(s,NM(q[i]),&u); s = u;      addn(s,NM(q[i]),&u); s = u;
         }    }
         for ( t = 0; i < n; i++ ) {    for ( t = 0; i < n; i++ ) {
                 addn(t,NM(q[i]),&u); t = u;      addn(t,NM(q[i]),&u); t = u;
         }    }
         gcdn(s,t,&g); NTOQ(g,1,*rp);    gcdn(s,t,&g); NTOQ(g,1,*rp);
 }  }
   
 void Pilcm(NODE arg,Obj *rp)  void Pilcm(NODE arg,Obj *rp)
 {  {
         N g,q,r,l;    N g,q,r,l;
         Q n1,n2,a;    Q n1,n2,a;
   
         n1 = (Q)ARG0(arg); n2 = (Q)ARG1(arg);    n1 = (Q)ARG0(arg); n2 = (Q)ARG1(arg);
         asir_assert(n1,O_N,"ilcm");    asir_assert(n1,O_N,"ilcm");
         asir_assert(n2,O_N,"ilcm");    asir_assert(n2,O_N,"ilcm");
         if ( !n1 || !n2 )    if ( !n1 || !n2 )
                 *rp = 0;      *rp = 0;
         else {    else {
                 gcdn(NM(n1),NM(n2),&g); divn(NM(n1),g,&q,&r);      gcdn(NM(n1),NM(n2),&g); divn(NM(n1),g,&q,&r);
                 muln(q,NM(n2),&l); NTOQ(l,1,a); *rp = (Obj)a;      muln(q,NM(n2),&l); NTOQ(l,1,a); *rp = (Obj)a;
         }    }
 }  }
   
 void Piand(NODE arg,Q *rp)  void Piand(NODE arg,Q *rp)
 {  {
         N g;    N g;
         Q n1,n2,a;    Q n1,n2,a;
   
         n1 = (Q)ARG0(arg); n2 = (Q)ARG1(arg);    n1 = (Q)ARG0(arg); n2 = (Q)ARG1(arg);
         asir_assert(n1,O_N,"iand");    asir_assert(n1,O_N,"iand");
         asir_assert(n2,O_N,"iand");    asir_assert(n2,O_N,"iand");
         if ( !n1 || !n2 )    if ( !n1 || !n2 )
                 *rp = 0;      *rp = 0;
         else {    else {
                 iand(NM(n1),NM(n2),&g);      iand(NM(n1),NM(n2),&g);
                 NTOQ(g,1,a); *rp = a;      NTOQ(g,1,a); *rp = a;
         }    }
 }  }
   
 void Pior(NODE arg,Q *rp)  void Pior(NODE arg,Q *rp)
 {  {
         N g;    N g;
         Q n1,n2,a;    Q n1,n2,a;
   
         n1 = (Q)ARG0(arg); n2 = (Q)ARG1(arg);    n1 = (Q)ARG0(arg); n2 = (Q)ARG1(arg);
         asir_assert(n1,O_N,"ior");    asir_assert(n1,O_N,"ior");
         asir_assert(n2,O_N,"ior");    asir_assert(n2,O_N,"ior");
         if ( !n1 )    if ( !n1 )
                 *rp = n2;      *rp = n2;
         else if ( !n2 )    else if ( !n2 )
                 *rp = n1;      *rp = n1;
         else {    else {
                 ior(NM(n1),NM(n2),&g);      ior(NM(n1),NM(n2),&g);
                 NTOQ(g,1,a); *rp = a;      NTOQ(g,1,a); *rp = a;
         }    }
 }  }
   
 void Pixor(NODE arg,Q *rp)  void Pixor(NODE arg,Q *rp)
 {  {
         N g;    N g;
         Q n1,n2,a;    Q n1,n2,a;
   
         n1 = (Q)ARG0(arg); n2 = (Q)ARG1(arg);    n1 = (Q)ARG0(arg); n2 = (Q)ARG1(arg);
         asir_assert(n1,O_N,"ixor");    asir_assert(n1,O_N,"ixor");
         asir_assert(n2,O_N,"ixor");    asir_assert(n2,O_N,"ixor");
         if ( !n1 )    if ( !n1 )
                 *rp = n2;      *rp = n2;
         else if ( !n2 )    else if ( !n2 )
                 *rp = n1;      *rp = n1;
         else {    else {
                 ixor(NM(n1),NM(n2),&g);      ixor(NM(n1),NM(n2),&g);
                 NTOQ(g,1,a); *rp = a;      NTOQ(g,1,a); *rp = a;
         }    }
 }  }
   
 void Pishift(NODE arg,Q *rp)  void Pishift(NODE arg,Q *rp)
 {  {
         N g;    N g;
         Q n1,s,a;    Q n1,s,a;
   
         n1 = (Q)ARG0(arg); s = (Q)ARG1(arg);    n1 = (Q)ARG0(arg); s = (Q)ARG1(arg);
         asir_assert(n1,O_N,"ixor");    asir_assert(n1,O_N,"ixor");
         asir_assert(s,O_N,"ixor");    asir_assert(s,O_N,"ixor");
         if ( !n1 )    if ( !n1 )
                 *rp = 0;      *rp = 0;
         else if ( !s )    else if ( !s )
                 *rp = n1;      *rp = n1;
         else {    else {
                 bshiftn(NM(n1),QTOS(s),&g);      bshiftn(NM(n1),QTOS(s),&g);
                 NTOQ(g,1,a); *rp = a;      NTOQ(g,1,a); *rp = a;
         }    }
 }  }
   
 void isqrt(N a,N *r)  void isqrt(N a,N *r)
 {  {
         int k;    int k;
         N x,t,x2,xh,quo,rem;    N x,t,x2,xh,quo,rem;
   
         if ( !a )    if ( !a )
                 *r = 0;      *r = 0;
         else if ( UNIN(a) )    else if ( UNIN(a) )
                 *r = ONEN;      *r = ONEN;
         else {    else {
                 k = n_bits(a); /* a <= 2^k-1 */      k = n_bits(a); /* a <= 2^k-1 */
                 bshiftn(ONEN,-((k>>1)+(k&1)),&x); /* a <= x^2 */      bshiftn(ONEN,-((k>>1)+(k&1)),&x); /* a <= x^2 */
                 while ( 1 ) {      while ( 1 ) {
                         pwrn(x,2,&t);        pwrn(x,2,&t);
                         if ( cmpn(t,a) <= 0 ) {        if ( cmpn(t,a) <= 0 ) {
                                 *r = x; return;          *r = x; return;
                         } else {        } else {
                                 if ( BD(x)[0] & 1 )          if ( BD(x)[0] & 1 )
                                         addn(x,a,&t);            addn(x,a,&t);
                                 else          else
                                         t = a;            t = a;
                                 bshiftn(x,-1,&x2); divn(t,x2,&quo,&rem);          bshiftn(x,-1,&x2); divn(t,x2,&quo,&rem);
                                 bshiftn(x,1,&xh); addn(quo,xh,&x);          bshiftn(x,1,&xh); addn(quo,xh,&x);
                         }        }
                 }      }
         }    }
 }  }
   
 void iand(N n1,N n2,N *r)  void iand(N n1,N n2,N *r)
 {  {
         int d1,d2,d,i;    int d1,d2,d,i;
         N nr;    N nr;
         int *p1,*p2,*pr;    int *p1,*p2,*pr;
   
         d1 = PL(n1); d2 = PL(n2);    d1 = PL(n1); d2 = PL(n2);
         d = MIN(d1,d2);    d = MIN(d1,d2);
         nr = NALLOC(d);    nr = NALLOC(d);
         for ( i = 0, p1 = BD(n1), p2 = BD(n2), pr = BD(nr); i < d; i++ )    for ( i = 0, p1 = BD(n1), p2 = BD(n2), pr = BD(nr); i < d; i++ )
                 pr[i] = p1[i] & p2[i];      pr[i] = p1[i] & p2[i];
         for ( i = d-1; i >= 0 && !pr[i]; i-- );    for ( i = d-1; i >= 0 && !pr[i]; i-- );
         if ( i < 0 )    if ( i < 0 )
                 *r = 0;      *r = 0;
         else {    else {
                 PL(nr) = i+1; *r = nr;      PL(nr) = i+1; *r = nr;
         }    }
 }  }
   
 void ior(N n1,N n2,N *r)  void ior(N n1,N n2,N *r)
 {  {
         int d1,d2,i;    int d1,d2,i;
         N nr;    N nr;
         int *p1,*p2,*pr;    int *p1,*p2,*pr;
   
         if ( PL(n1) < PL(n2) ) {    if ( PL(n1) < PL(n2) ) {
                 nr = n1; n1 = n2; n2 = nr;      nr = n1; n1 = n2; n2 = nr;
         }    }
         d1 = PL(n1); d2 = PL(n2);    d1 = PL(n1); d2 = PL(n2);
         *r = nr = NALLOC(d1);    *r = nr = NALLOC(d1);
         for ( i = 0, p1 = BD(n1), p2 = BD(n2), pr = BD(nr); i < d2; i++ )    for ( i = 0, p1 = BD(n1), p2 = BD(n2), pr = BD(nr); i < d2; i++ )
                 pr[i] = p1[i] | p2[i];      pr[i] = p1[i] | p2[i];
         for ( ; i < d1; i++ )    for ( ; i < d1; i++ )
                 pr[i] = p1[i];      pr[i] = p1[i];
         for ( i = d1-1; i >= 0 && !pr[i]; i-- );    for ( i = d1-1; i >= 0 && !pr[i]; i-- );
         if ( i < 0 )    if ( i < 0 )
                 *r = 0;      *r = 0;
         else {    else {
                 PL(nr) = i+1; *r = nr;      PL(nr) = i+1; *r = nr;
         }    }
 }  }
   
 void ixor(N n1,N n2,N *r)  void ixor(N n1,N n2,N *r)
 {  {
         int d1,d2,i;    int d1,d2,i;
         N nr;    N nr;
         int *p1,*p2,*pr;    int *p1,*p2,*pr;
   
         if ( PL(n1) < PL(n2) ) {    if ( PL(n1) < PL(n2) ) {
                 nr = n1; n1 = n2; n2 = nr;      nr = n1; n1 = n2; n2 = nr;
         }    }
         d1 = PL(n1); d2 = PL(n2);    d1 = PL(n1); d2 = PL(n2);
         *r = nr = NALLOC(d1);    *r = nr = NALLOC(d1);
         for ( i = 0, p1 = BD(n1), p2 = BD(n2), pr = BD(nr); i < d2; i++ )    for ( i = 0, p1 = BD(n1), p2 = BD(n2), pr = BD(nr); i < d2; i++ )
                 pr[i] = p1[i] ^ p2[i];      pr[i] = p1[i] ^ p2[i];
         for ( ; i < d1; i++ )    for ( ; i < d1; i++ )
                 pr[i] = p1[i];      pr[i] = p1[i];
         for ( i = d1-1; i >= 0 && !pr[i]; i-- );    for ( i = d1-1; i >= 0 && !pr[i]; i-- );
         if ( i < 0 )    if ( i < 0 )
                 *r = 0;      *r = 0;
         else {    else {
                 PL(nr) = i+1; *r = nr;      PL(nr) = i+1; *r = nr;
         }    }
 }  }
   
 void Pup2_init_eg(Obj *rp)  void Pup2_init_eg(Obj *rp)
 {  {
         up2_init_eg();    up2_init_eg();
         *rp = 0;    *rp = 0;
 }  }
   
 void Pup2_show_eg(Obj *rp)  void Pup2_show_eg(Obj *rp)
 {  {
         up2_show_eg();    up2_show_eg();
         *rp = 0;    *rp = 0;
 }  }
   
 void Pgf2nton(NODE arg,Q *rp)  void Pgf2nton(NODE arg,Q *rp)
 {  {
         if ( !ARG0(arg) )    if ( !ARG0(arg) )
                 *rp = 0;      *rp = 0;
         else    else
                 up2ton(((GF2N)ARG0(arg))->body,rp);      up2ton(((GF2N)ARG0(arg))->body,rp);
 }  }
   
 void Pntogf2n(NODE arg,GF2N *rp)  void Pntogf2n(NODE arg,GF2N *rp)
 {  {
         UP2 t;    UP2 t;
   
         ntoup2(ARG0(arg),&t);    ntoup2(ARG0(arg),&t);
         MKGF2N(t,*rp);    MKGF2N(t,*rp);
 }  }
   
 void Pup2_inv(NODE arg,P *rp)  void Pup2_inv(NODE arg,P *rp)
 {  {
         UP2 a,b,t;    UP2 a,b,t;
   
         ptoup2(ARG0(arg),&a);    ptoup2(ARG0(arg),&a);
         ptoup2(ARG1(arg),&b);    ptoup2(ARG1(arg),&b);
         invup2(a,b,&t);    invup2(a,b,&t);
         up2top(t,rp);    up2top(t,rp);
 }  }
   
 void Pinv(NODE arg,Num *rp)  void Pinv(NODE arg,Num *rp)
 {  {
         Num n;    Num n;
         Q mod;    Q mod;
         MQ r;    MQ r;
         int inv;    int inv;
   
         n = (Num)ARG0(arg); mod = (Q)ARG1(arg);    n = (Num)ARG0(arg); mod = (Q)ARG1(arg);
         asir_assert(n,O_N,"inv");    asir_assert(n,O_N,"inv");
         asir_assert(mod,O_N,"inv");    asir_assert(mod,O_N,"inv");
         if ( !n || !mod )    if ( !n || !mod )
                 error("inv: invalid input");      error("inv: invalid input");
         else    else
                 switch ( NID(n) ) {      switch ( NID(n) ) {
                         case N_Q:        case N_Q:
                                 invl((Q)n,mod,(Q *)rp);          invl((Q)n,mod,(Q *)rp);
                                 break;          break;
                         case N_M:        case N_M:
                                 inv = invm(CONT((MQ)n),QTOS(mod));          inv = invm(CONT((MQ)n),QTOS(mod));
                                 STOMQ(inv,r);          STOMQ(inv,r);
                                 *rp = (Num)r;          *rp = (Num)r;
                                 break;          break;
                         default:        default:
                                 error("inv: invalid input");          error("inv: invalid input");
                 }      }
 }  }
   
 void Pfac(NODE arg,Q *rp)  void Pfac(NODE arg,Q *rp)
 {  {
         asir_assert(ARG0(arg),O_N,"fac");    asir_assert(ARG0(arg),O_N,"fac");
         factorial(QTOS((Q)ARG0(arg)),rp);    factorial(QTOS((Q)ARG0(arg)),rp);
 }  }
   
 void Plrandom(NODE arg,Q *rp)  void Plrandom(NODE arg,Q *rp)
 {  {
         N r;    N r;
   
         asir_assert(ARG0(arg),O_N,"lrandom");    asir_assert(ARG0(arg),O_N,"lrandom");
         randomn(QTOS((Q)ARG0(arg)),&r);    randomn(QTOS((Q)ARG0(arg)),&r);
         NTOQ(r,1,*rp);    NTOQ(r,1,*rp);
 }  }
   
 void Prandom(NODE arg,Q *rp)  void Prandom(NODE arg,Q *rp)
 {  {
         unsigned int r;    unsigned int r;
   
 #if 0  #if 0
 #if defined(_PA_RISC1_1)  #if defined(_PA_RISC1_1)
         r = mrand48()&BMASK;    r = mrand48()&BMASK;
 #else  #else
         if ( arg )    if ( arg )
                 srandom(QTOS((Q)ARG0(arg)));      srandom(QTOS((Q)ARG0(arg)));
         r = random()&BMASK;    r = random()&BMASK;
 #endif  #endif
 #endif  #endif
         if ( arg )    if ( arg )
                 mt_sgenrand(QTOS((Q)ARG0(arg)));      mt_sgenrand(QTOS((Q)ARG0(arg)));
         r = mt_genrand();    r = mt_genrand();
         UTOQ(r,*rp);    UTOQ(r,*rp);
 }  }
   
 #if defined(VISUAL) || defined(__MINGW32__)  #if defined(VISUAL) || defined(__MINGW32__)
Line 841  unsigned int random() {
Line 841  unsigned int random() {
   
 void srandom(unsigned int s)  void srandom(unsigned int s)
 {  {
                 if ( s )      if ( s )
                         R_Next = s;        R_Next = s;
         else if ( !R_Next )          else if ( !R_Next )
             R_Next = 1;              R_Next = 1;
 }  }
Line 850  void srandom(unsigned int s)
Line 850  void srandom(unsigned int s)
   
 void Pprime(NODE arg,Q *rp)  void Pprime(NODE arg,Q *rp)
 {  {
         int i;    int i;
   
         asir_assert(ARG0(arg),O_N,"prime");    asir_assert(ARG0(arg),O_N,"prime");
         i = QTOS((Q)ARG0(arg));    i = QTOS((Q)ARG0(arg));
         if ( i < 0 || i >= 1900 )    if ( i < 0 || i >= 1900 )
                 *rp = 0;      *rp = 0;
         else    else
                 STOQ(sprime[i],*rp);      STOQ(sprime[i],*rp);
 }  }
   
 void Plprime(NODE arg,Q *rp)  void Plprime(NODE arg,Q *rp)
 {  {
         int i,p;    int i,p;
   
         asir_assert(ARG0(arg),O_N,"lprime");    asir_assert(ARG0(arg),O_N,"lprime");
         i = QTOS((Q)ARG0(arg));    i = QTOS((Q)ARG0(arg));
         if ( i < 0 )    if ( i < 0 )
                 *rp = 0;      *rp = 0;
         else    else
                 p = get_lprime(i);      p = get_lprime(i);
         STOQ(p,*rp);    STOQ(p,*rp);
 }  }
   
 extern int up_kara_mag, up_tkara_mag, up_fft_mag;  extern int up_kara_mag, up_tkara_mag, up_fft_mag;
   
 void Pset_upfft(NODE arg,Q *rp)  void Pset_upfft(NODE arg,Q *rp)
 {  {
         if ( arg ) {    if ( arg ) {
                 asir_assert(ARG0(arg),O_N,"set_upfft");      asir_assert(ARG0(arg),O_N,"set_upfft");
                 up_fft_mag = QTOS((Q)ARG0(arg));      up_fft_mag = QTOS((Q)ARG0(arg));
         }    }
         STOQ(up_fft_mag,*rp);    STOQ(up_fft_mag,*rp);
 }  }
   
 void Pset_upkara(NODE arg,Q *rp)  void Pset_upkara(NODE arg,Q *rp)
 {  {
         if ( arg ) {    if ( arg ) {
                 asir_assert(ARG0(arg),O_N,"set_upkara");      asir_assert(ARG0(arg),O_N,"set_upkara");
                 up_kara_mag = QTOS((Q)ARG0(arg));      up_kara_mag = QTOS((Q)ARG0(arg));
         }    }
         STOQ(up_kara_mag,*rp);    STOQ(up_kara_mag,*rp);
 }  }
   
 void Pset_uptkara(NODE arg,Q *rp)  void Pset_uptkara(NODE arg,Q *rp)
 {  {
         if ( arg ) {    if ( arg ) {
                 asir_assert(ARG0(arg),O_N,"set_uptkara");      asir_assert(ARG0(arg),O_N,"set_uptkara");
                 up_tkara_mag = QTOS((Q)ARG0(arg));      up_tkara_mag = QTOS((Q)ARG0(arg));
         }    }
         STOQ(up_tkara_mag,*rp);    STOQ(up_tkara_mag,*rp);
 }  }
   
 extern int up2_kara_mag;  extern int up2_kara_mag;
   
 void Pset_up2kara(NODE arg,Q *rp)  void Pset_up2kara(NODE arg,Q *rp)
 {  {
         if ( arg ) {    if ( arg ) {
                 asir_assert(ARG0(arg),O_N,"set_up2kara");      asir_assert(ARG0(arg),O_N,"set_up2kara");
                 up2_kara_mag = QTOS((Q)ARG0(arg));      up2_kara_mag = QTOS((Q)ARG0(arg));
         }    }
         STOQ(up2_kara_mag,*rp);    STOQ(up2_kara_mag,*rp);
 }  }
   
 void Pigcdbin(NODE arg,Obj *rp)  void Pigcdbin(NODE arg,Obj *rp)
 {  {
         N g;    N g;
         Q n1,n2,a;    Q n1,n2,a;
   
         n1 = (Q)ARG0(arg); n2 = (Q)ARG1(arg);    n1 = (Q)ARG0(arg); n2 = (Q)ARG1(arg);
         asir_assert(n1,O_N,"igcd");    asir_assert(n1,O_N,"igcd");
         asir_assert(n2,O_N,"igcd");    asir_assert(n2,O_N,"igcd");
         if ( !n1 )    if ( !n1 )
                 *rp = (Obj)n2;      *rp = (Obj)n2;
         else if ( !n2 )    else if ( !n2 )
                 *rp = (Obj)n1;      *rp = (Obj)n1;
         else {    else {
                 gcdbinn(NM(n1),NM(n2),&g);      gcdbinn(NM(n1),NM(n2),&g);
                 NTOQ(g,1,a); *rp = (Obj)a;      NTOQ(g,1,a); *rp = (Obj)a;
         }    }
 }  }
   
 void Pigcdbmod(NODE arg,Obj *rp)  void Pigcdbmod(NODE arg,Obj *rp)
 {  {
         N g;    N g;
         Q n1,n2,a;    Q n1,n2,a;
   
         n1 = (Q)ARG0(arg); n2 = (Q)ARG1(arg);    n1 = (Q)ARG0(arg); n2 = (Q)ARG1(arg);
         asir_assert(n1,O_N,"igcdbmod");    asir_assert(n1,O_N,"igcdbmod");
         asir_assert(n2,O_N,"igcdbmod");    asir_assert(n2,O_N,"igcdbmod");
         if ( !n1 )    if ( !n1 )
                 *rp = (Obj)n2;      *rp = (Obj)n2;
         else if ( !n2 )    else if ( !n2 )
                 *rp = (Obj)n1;      *rp = (Obj)n1;
         else {    else {
                 gcdbmodn(NM(n1),NM(n2),&g);      gcdbmodn(NM(n1),NM(n2),&g);
                 NTOQ(g,1,a); *rp = (Obj)a;      NTOQ(g,1,a); *rp = (Obj)a;
         }    }
 }  }
   
 void Pigcdacc(NODE arg,Obj *rp)  void Pigcdacc(NODE arg,Obj *rp)
 {  {
         N g;    N g;
         Q n1,n2,a;    Q n1,n2,a;
   
         n1 = (Q)ARG0(arg); n2 = (Q)ARG1(arg);    n1 = (Q)ARG0(arg); n2 = (Q)ARG1(arg);
         asir_assert(n1,O_N,"igcdacc");    asir_assert(n1,O_N,"igcdacc");
         asir_assert(n2,O_N,"igcdacc");    asir_assert(n2,O_N,"igcdacc");
         if ( !n1 )    if ( !n1 )
                 *rp = (Obj)n2;      *rp = (Obj)n2;
         else if ( !n2 )    else if ( !n2 )
                 *rp = (Obj)n1;      *rp = (Obj)n1;
         else {    else {
                 gcdaccn(NM(n1),NM(n2),&g);      gcdaccn(NM(n1),NM(n2),&g);
                 NTOQ(g,1,a); *rp = (Obj)a;      NTOQ(g,1,a); *rp = (Obj)a;
         }    }
 }  }
   
 void PigcdEuc(NODE arg,Obj *rp)  void PigcdEuc(NODE arg,Obj *rp)
 {  {
         N g;    N g;
         Q n1,n2,a;    Q n1,n2,a;
   
         n1 = (Q)ARG0(arg); n2 = (Q)ARG1(arg);    n1 = (Q)ARG0(arg); n2 = (Q)ARG1(arg);
         asir_assert(n1,O_N,"igcdbmod");    asir_assert(n1,O_N,"igcdbmod");
         asir_assert(n2,O_N,"igcdbmod");    asir_assert(n2,O_N,"igcdbmod");
         if ( !n1 )    if ( !n1 )
                 *rp = (Obj)n2;      *rp = (Obj)n2;
         else if ( !n2 )    else if ( !n2 )
                 *rp = (Obj)n1;      *rp = (Obj)n1;
         else {    else {
                 gcdEuclidn(NM(n1),NM(n2),&g);      gcdEuclidn(NM(n1),NM(n2),&g);
                 NTOQ(g,1,a); *rp = (Obj)a;      NTOQ(g,1,a); *rp = (Obj)a;
         }    }
 }  }
   
 extern int igcd_algorithm;  extern int igcd_algorithm;
     /*  == 0 : Euclid,      /*  == 0 : Euclid,
      *  == 1 : binary,       *  == 1 : binary,
      *  == 2 : bmod,       *  == 2 : bmod,
      *  >= 3 : (Weber's accelerated)/(Jebelean's generalized binary) algorithm,       *  >= 3 : (Weber's accelerated)/(Jebelean's generalized binary) algorithm,
      */       */
 extern int igcd_thre_inidiv;  extern int igcd_thre_inidiv;
     /*      /*
Line 1006  extern int igcdacc_thre;
Line 1006  extern int igcdacc_thre;
   
 void Pigcdcntl(NODE arg,Obj *rp)  void Pigcdcntl(NODE arg,Obj *rp)
 {  {
         Obj p;    Obj p;
         Q a;    Q a;
         int k, m;    int k, m;
   
         if ( arg ) {    if ( arg ) {
                 p = (Obj)ARG0(arg);      p = (Obj)ARG0(arg);
                 if ( !p ) {      if ( !p ) {
                         igcd_algorithm = 0;        igcd_algorithm = 0;
                         *rp = p;        *rp = p;
                         return;        return;
                 } else if ( OID(p) == O_N ) {      } else if ( OID(p) == O_N ) {
                         k = QTOS((Q)p);        k = QTOS((Q)p);
                         a = (Q)p;        a = (Q)p;
                         if ( k >= 0 ) igcd_algorithm = k;        if ( k >= 0 ) igcd_algorithm = k;
                         else if ( k == -1 ) {        else if ( k == -1 ) {
                         ret_thre:        ret_thre:
                                 k = - igcd_thre_inidiv - igcdacc_thre*10000;          k = - igcd_thre_inidiv - igcdacc_thre*10000;
                                 STOQ(k,a);          STOQ(k,a);
                                 *rp = (Obj)a;          *rp = (Obj)a;
                                 return;          return;
                         } else {        } else {
                                 if ( (m = (-k)%10000) != 0 ) igcd_thre_inidiv = m;          if ( (m = (-k)%10000) != 0 ) igcd_thre_inidiv = m;
                                 if ( (m = -k/10000) != 0 ) igcdacc_thre = m;          if ( (m = -k/10000) != 0 ) igcdacc_thre = m;
                                 goto ret_thre;          goto ret_thre;
                         }        }
                 } else if ( OID(p) == O_STR ) {      } else if ( OID(p) == O_STR ) {
                         char *n = BDY((STRING) p);        char *n = BDY((STRING) p);
   
                         if ( !strcmp( n, "binary" ) || !strcmp( n, "Binary" )        if ( !strcmp( n, "binary" ) || !strcmp( n, "Binary" )
                              || !strcmp( n, "bin" ) || !strcmp( n, "Bin" ) )             || !strcmp( n, "bin" ) || !strcmp( n, "Bin" ) )
                                 k = igcd_algorithm = 1;          k = igcd_algorithm = 1;
                         else if ( !strcmp( n, "bmod" ) || !strcmp( n, "Bmod" ) )        else if ( !strcmp( n, "bmod" ) || !strcmp( n, "Bmod" ) )
                                 igcd_algorithm = 2;          igcd_algorithm = 2;
                         else if ( !strcmp( n, "euc" ) || !strcmp( n, "Euc" )        else if ( !strcmp( n, "euc" ) || !strcmp( n, "Euc" )
                                   || !strcmp( n, "euclid" ) || !strcmp( n, "Euclid" ) )            || !strcmp( n, "euclid" ) || !strcmp( n, "Euclid" ) )
                                 igcd_algorithm = 0;          igcd_algorithm = 0;
                         else if ( !strcmp( n, "acc" ) || !strcmp( n, "Acc" )        else if ( !strcmp( n, "acc" ) || !strcmp( n, "Acc" )
                              || !strcmp( n, "gen" ) || !strcmp( n, "Gen" )             || !strcmp( n, "gen" ) || !strcmp( n, "Gen" )
                              || !strcmp( n, "genbin" ) || !strcmp( n, "GenBin" ) )             || !strcmp( n, "genbin" ) || !strcmp( n, "GenBin" ) )
                                 igcd_algorithm = 3;          igcd_algorithm = 3;
                         else error( "igcdhow : invalid algorithm specification" );        else error( "igcdhow : invalid algorithm specification" );
                 }      }
         }    }
         STOQ(igcd_algorithm,a);    STOQ(igcd_algorithm,a);
         *rp = (Obj)a;    *rp = (Obj)a;
         return;    return;
 }  }
Legend:
Removed from v.1.13  
changed lines
  Added in v.1.14
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>