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

Diff for /OpenXM_contrib2/asir2000/engine/gmpq.c between version 1.9 and 1.10

version 1.9, 2017/09/15 01:52:51 version 1.10, 2018/03/29 01:32:52
Line 14  void bshiftgz(GZ a,int n,GZ *r);
Line 14  void bshiftgz(GZ a,int n,GZ *r);
   
 void *gc_realloc(void *p,size_t osize,size_t nsize)  void *gc_realloc(void *p,size_t osize,size_t nsize)
 {  {
         return (void *)Risa_GC_realloc(p,nsize);    return (void *)Risa_GC_realloc(p,nsize);
 }  }
   
 void gc_free(void *p,size_t size)  void gc_free(void *p,size_t size)
 {  {
         Risa_GC_free(p);    Risa_GC_free(p);
 }  }
   
 void init_gmpq()  void init_gmpq()
 {  {
         mp_set_memory_functions(Risa_GC_malloc,gc_realloc,gc_free);    mp_set_memory_functions(Risa_GC_malloc,gc_realloc,gc_free);
   
         mpz_init(ONEMPZ); mpz_set_ui(ONEMPZ,1); MPZTOGZ(ONEMPZ,ONEGZ);    mpz_init(ONEMPZ); mpz_set_ui(ONEMPZ,1); MPZTOGZ(ONEMPZ,ONEGZ);
 }  }
   
 GZ utogz(unsigned int u)  GZ utogz(unsigned int u)
 {  {
         mpz_t z;    mpz_t z;
         GZ r;    GZ r;
   
         if ( !u ) return 0;    if ( !u ) return 0;
         mpz_init(z); mpz_set_ui(z,u); MPZTOGZ(z,r); return r;    mpz_init(z); mpz_set_ui(z,u); MPZTOGZ(z,r); return r;
 }  }
   
 GZ stogz(int s)  GZ stogz(int s)
 {  {
         mpz_t z;    mpz_t z;
         GZ r;    GZ r;
   
         if ( !s ) return 0;    if ( !s ) return 0;
         mpz_init(z); mpz_set_si(z,s); MPZTOGZ(z,r); return r;    mpz_init(z); mpz_set_si(z,s); MPZTOGZ(z,r); return r;
 }  }
   
 GQ mpqtogzq(mpq_t a)  GQ mpqtogzq(mpq_t a)
 {  {
         GZ z;    GZ z;
         GQ q;    GQ q;
   
         if ( INTMPQ(a) ) {    if ( INTMPQ(a) ) {
                 MPZTOGZ(mpq_numref(a),z); return (GQ)z;      MPZTOGZ(mpq_numref(a),z); return (GQ)z;
         } else {    } else {
                 MPQTOGQ(a,q); return q;      MPQTOGQ(a,q); return q;
         }    }
 }  }
   
 GZ ztogz(Q a)  GZ ztogz(Q a)
 {  {
         mpz_t z;    mpz_t z;
         mpq_t b;    mpq_t b;
         N nm;    N nm;
         GZ s;    GZ s;
   
         if ( !a ) return 0;    if ( !a ) return 0;
         nm = NM(a);    nm = NM(a);
         mpz_init(z);    mpz_init(z);
         mpz_import(z,PL(nm),-1,sizeof(BD(nm)[0]),0,0,BD(nm));    mpz_import(z,PL(nm),-1,sizeof(BD(nm)[0]),0,0,BD(nm));
         if ( SGN(a)<0 ) mpz_neg(z,z);    if ( SGN(a)<0 ) mpz_neg(z,z);
         MPZTOGZ(z,s);    MPZTOGZ(z,s);
         return s;    return s;
 }  }
   
 Q gztoz(GZ a)  Q gztoz(GZ a)
 {  {
         N nm;    N nm;
         Q q;    Q q;
         int sgn;    int sgn;
         size_t len;    size_t len;
   
         if ( !a ) return 0;    if ( !a ) return 0;
         len = WORDSIZE_IN_N(BDY((GZ)a)); nm = NALLOC(len);    len = WORDSIZE_IN_N(BDY((GZ)a)); nm = NALLOC(len);
         mpz_export(BD(nm),&len,-1,sizeof(int),0,0,BDY((GZ)a));    mpz_export(BD(nm),&len,-1,sizeof(int),0,0,BDY((GZ)a));
         PL(nm) = len;    PL(nm) = len;
         sgn = mpz_sgn(BDY((GZ)a)); NTOQ(nm,sgn,q);    sgn = mpz_sgn(BDY((GZ)a)); NTOQ(nm,sgn,q);
         return q;    return q;
 }  }
   
 void dupgz(GZ a,GZ *b)  void dupgz(GZ a,GZ *b)
Line 102  void dupgz(GZ a,GZ *b)
Line 102  void dupgz(GZ a,GZ *b)
   
 int n_bits_gz(GZ a)  int n_bits_gz(GZ a)
 {  {
         return a ? mpz_sizeinbase(BDY(a),2) : 0;    return a ? mpz_sizeinbase(BDY(a),2) : 0;
 }  }
   
 GQ qtogq(Q a)  GQ qtogq(Q a)
 {  {
         mpz_t z;    mpz_t z;
         mpq_t b;    mpq_t b;
         N nm,dn;    N nm,dn;
         GZ s;    GZ s;
         GQ r;    GQ r;
   
         if ( !a ) return 0;    if ( !a ) return 0;
         if ( INT(a) ) {    if ( INT(a) ) {
                 nm = NM(a);      nm = NM(a);
                 mpz_init(z);      mpz_init(z);
                 mpz_import(z,PL(nm),-1,sizeof(BD(nm)[0]),0,0,BD(nm));      mpz_import(z,PL(nm),-1,sizeof(BD(nm)[0]),0,0,BD(nm));
                 if ( SGN(a)<0 ) mpz_neg(z,z);      if ( SGN(a)<0 ) mpz_neg(z,z);
                 MPZTOGZ(z,s);      MPZTOGZ(z,s);
                 return (GQ)s;      return (GQ)s;
         } else {    } else {
                 nm = NM(a); dn = DN(a);      nm = NM(a); dn = DN(a);
                 mpq_init(b);      mpq_init(b);
                 mpz_import(mpq_numref(b),PL(nm),-1,sizeof(BD(nm)[0]),0,0,BD(nm));      mpz_import(mpq_numref(b),PL(nm),-1,sizeof(BD(nm)[0]),0,0,BD(nm));
                 mpz_import(mpq_denref(b),PL(dn),-1,sizeof(BD(dn)[0]),0,0,BD(dn));      mpz_import(mpq_denref(b),PL(dn),-1,sizeof(BD(dn)[0]),0,0,BD(dn));
                 if ( SGN(a)<0 ) mpq_neg(b,b);      if ( SGN(a)<0 ) mpq_neg(b,b);
                 MPQTOGQ(b,r);      MPQTOGQ(b,r);
                 return r;      return r;
         }    }
 }  }
   
 Q gqtoq(GQ a)  Q gqtoq(GQ a)
 {  {
         N nm,dn;    N nm,dn;
         Q q;    Q q;
         int sgn;    int sgn;
         size_t len;    size_t len;
   
         if ( !a ) return 0;    if ( !a ) return 0;
         if ( NID(a) == N_GZ ) {    if ( NID(a) == N_GZ ) {
                 len = WORDSIZE_IN_N(BDY((GZ)a)); nm = NALLOC(len);      len = WORDSIZE_IN_N(BDY((GZ)a)); nm = NALLOC(len);
                 mpz_export(BD(nm),&len,-1,sizeof(int),0,0,BDY((GZ)a));      mpz_export(BD(nm),&len,-1,sizeof(int),0,0,BDY((GZ)a));
                 PL(nm) = len;      PL(nm) = len;
                 sgn = mpz_sgn(BDY((GZ)a)); NTOQ(nm,sgn,q);      sgn = mpz_sgn(BDY((GZ)a)); NTOQ(nm,sgn,q);
         } else {    } else {
                 len = WORDSIZE_IN_N(mpq_numref(BDY(a))); nm = NALLOC(len);      len = WORDSIZE_IN_N(mpq_numref(BDY(a))); nm = NALLOC(len);
                 mpz_export(BD(nm),&len,-1,sizeof(int),0,0,mpq_numref(BDY(a)));      mpz_export(BD(nm),&len,-1,sizeof(int),0,0,mpq_numref(BDY(a)));
                 PL(nm) = len;      PL(nm) = len;
                 len = WORDSIZE_IN_N(mpq_denref(BDY(a))); dn = NALLOC(len);      len = WORDSIZE_IN_N(mpq_denref(BDY(a))); dn = NALLOC(len);
                 mpz_export(BD(dn),&len,-1,sizeof(int),0,0,mpq_denref(BDY(a)));      mpz_export(BD(dn),&len,-1,sizeof(int),0,0,mpq_denref(BDY(a)));
                 PL(dn) = len;      PL(dn) = len;
                 sgn = mpz_sgn(mpq_numref(BDY(a))); NDTOQ(nm,dn,sgn,q);      sgn = mpz_sgn(mpq_numref(BDY(a))); NDTOQ(nm,dn,sgn,q);
         }    }
         return q;    return q;
 }  }
   
 P ptogp(P a)  P ptogp(P a)
 {  {
         DCP dc,dcr,dcr0;    DCP dc,dcr,dcr0;
         P b;    P b;
   
         if ( !a ) return 0;    if ( !a ) return 0;
         if ( NUM(a) ) return (P)qtogq((Q)a);    if ( NUM(a) ) return (P)qtogq((Q)a);
         for ( dc = DC(a), dcr0 = 0; dc; dc = NEXT(dc) ) {    for ( dc = DC(a), dcr0 = 0; dc; dc = NEXT(dc) ) {
                 NEXTDC(dcr0,dcr); DEG(dcr) = DEG(dc); COEF(dcr) = (P)ptogp(COEF(dc));      NEXTDC(dcr0,dcr); DEG(dcr) = DEG(dc); COEF(dcr) = (P)ptogp(COEF(dc));
         }    }
         NEXT(dcr) = 0; MKP(VR(a),dcr0,b);    NEXT(dcr) = 0; MKP(VR(a),dcr0,b);
         return b;    return b;
 }  }
   
 P gptop(P a)  P gptop(P a)
 {  {
         DCP dc,dcr,dcr0;    DCP dc,dcr,dcr0;
         P b;    P b;
   
         if ( !a ) return 0;    if ( !a ) return 0;
         if ( NUM(a) ) b = (P)gqtoq((GQ)a);    if ( NUM(a) ) b = (P)gqtoq((GQ)a);
         else {    else {
                 for ( dc = DC(a), dcr0 = 0; dc; dc = NEXT(dc) ) {      for ( dc = DC(a), dcr0 = 0; dc; dc = NEXT(dc) ) {
                         NEXTDC(dcr0,dcr); DEG(dcr) = DEG(dc);        NEXTDC(dcr0,dcr); DEG(dcr) = DEG(dc);
                         COEF(dcr) = (P)gptop(COEF(dc));        COEF(dcr) = (P)gptop(COEF(dc));
                 }      }
                 NEXT(dcr) = 0; MKP(VR(a),dcr0,b);      NEXT(dcr) = 0; MKP(VR(a),dcr0,b);
         }    }
         return b;    return b;
 }  }
   
 void addgz(GZ n1,GZ n2,GZ *nr)  void addgz(GZ n1,GZ n2,GZ *nr)
 {  {
         mpz_t t;    mpz_t t;
         int s1,s2;    int s1,s2;
   
         if ( !n1 ) *nr = n2;    if ( !n1 ) *nr = n2;
         else if ( !n2 ) *nr = n1;    else if ( !n2 ) *nr = n1;
         else {    else {
                 mpz_init(t); mpz_add(t,BDY(n1),BDY(n2)); MPZTOGZ(t,*nr);      mpz_init(t); mpz_add(t,BDY(n1),BDY(n2)); MPZTOGZ(t,*nr);
         }    }
 }  }
   
 void subgz(GZ n1,GZ n2,GZ *nr)  void subgz(GZ n1,GZ n2,GZ *nr)
 {  {
         mpz_t t;    mpz_t t;
   
         if ( !n1 )    if ( !n1 )
                 if ( !n2 )      if ( !n2 )
                         *nr = 0;        *nr = 0;
                 else {      else {
                         t[0] = BDY(n2)[0]; mpz_neg(t,t); MPZTOGZ(t,*nr);        t[0] = BDY(n2)[0]; mpz_neg(t,t); MPZTOGZ(t,*nr);
                 }      }
         else if ( !n2 )    else if ( !n2 )
                 *nr = n1;      *nr = n1;
         else if ( n1 == n2 )    else if ( n1 == n2 )
                 *nr = 0;      *nr = 0;
         else {    else {
                 mpz_init(t); mpz_sub(t,BDY(n1),BDY(n2)); MPZTOGZ(t,*nr);      mpz_init(t); mpz_sub(t,BDY(n1),BDY(n2)); MPZTOGZ(t,*nr);
         }    }
 }  }
   
 void mulgz(GZ n1,GZ n2,GZ *nr)  void mulgz(GZ n1,GZ n2,GZ *nr)
 {  {
         mpz_t t;    mpz_t t;
   
         if ( !n1 || !n2 ) *nr = 0;    if ( !n1 || !n2 ) *nr = 0;
 #if 1  #if 1
         else if ( UNIGZ(n1) ) *nr = n2;    else if ( UNIGZ(n1) ) *nr = n2;
         else if ( UNIGZ(n2) ) *nr = n1;    else if ( UNIGZ(n2) ) *nr = n1;
         else if ( MUNIGZ(n1) ) chsgngz(n2,nr);    else if ( MUNIGZ(n1) ) chsgngz(n2,nr);
         else if ( MUNIGZ(n2) ) chsgngz(n1,nr);    else if ( MUNIGZ(n2) ) chsgngz(n1,nr);
 #endif  #endif
         else {    else {
                 mpz_init(t); mpz_mul(t,BDY(n1),BDY(n2)); MPZTOGZ(t,*nr);      mpz_init(t); mpz_mul(t,BDY(n1),BDY(n2)); MPZTOGZ(t,*nr);
         }    }
 }  }
   
 /* nr += n1*n2 */  /* nr += n1*n2 */
Line 241  void muladdtogz(GZ n1,GZ n2,GZ *nr)
Line 241  void muladdtogz(GZ n1,GZ n2,GZ *nr)
 {  {
     GZ t;      GZ t;
   
         if ( n1 && n2 ) {    if ( n1 && n2 ) {
         if ( !(*nr) ) {          if ( !(*nr) ) {
           NEWGZ(t); mpz_init(BDY(t)); *nr = t;            NEWGZ(t); mpz_init(BDY(t)); *nr = t;
         }          }
Line 251  void muladdtogz(GZ n1,GZ n2,GZ *nr)
Line 251  void muladdtogz(GZ n1,GZ n2,GZ *nr)
   
 void mul1gz(GZ n1,int n2,GZ *nr)  void mul1gz(GZ n1,int n2,GZ *nr)
 {  {
         mpz_t t;    mpz_t t;
   
         if ( !n1 || !n2 ) *nr = 0;    if ( !n1 || !n2 ) *nr = 0;
         else {    else {
                 mpz_init(t); mpz_mul_ui(t,BDY(n1),(long)n2); MPZTOGZ(t,*nr);      mpz_init(t); mpz_mul_ui(t,BDY(n1),(long)n2); MPZTOGZ(t,*nr);
         }    }
 }  }
   
 void divgz(GZ n1,GZ n2,GZ *nq)  void divgz(GZ n1,GZ n2,GZ *nq)
 {  {
         mpz_t t;    mpz_t t;
         mpq_t a, b, q;    mpq_t a, b, q;
   
         if ( !n2 ) {    if ( !n2 ) {
                 error("division by 0");      error("division by 0");
                 *nq = 0;      *nq = 0;
         } else if ( !n1 )    } else if ( !n1 )
                 *nq = 0;      *nq = 0;
         else if ( n1 == n2 ) {    else if ( n1 == n2 ) {
                 mpz_init(t); mpz_set_ui(t,1); MPZTOGZ(t,*nq);      mpz_init(t); mpz_set_ui(t,1); MPZTOGZ(t,*nq);
         } else {    } else {
                 MPZTOMPQ(BDY(n1),a); MPZTOMPQ(BDY(n2),b);      MPZTOMPQ(BDY(n1),a); MPZTOMPQ(BDY(n2),b);
                 mpq_init(q); mpq_div(q,a,b); *nq = (GZ)mpqtogzq(q);      mpq_init(q); mpq_div(q,a,b); *nq = (GZ)mpqtogzq(q);
         }    }
 }  }
   
 void remgz(GZ n1,GZ n2,GZ *nr)  void remgz(GZ n1,GZ n2,GZ *nr)
 {  {
         mpz_t r;    mpz_t r;
   
         if ( !n2 ) {    if ( !n2 ) {
                 error("division by 0");      error("division by 0");
                 *nr = 0;      *nr = 0;
         } else if ( !n1 || n1 == n2 )    } else if ( !n1 || n1 == n2 )
                 *nr = 0;      *nr = 0;
         else {    else {
                 mpz_init(r);      mpz_init(r);
                 mpz_mod(r,BDY(n1),BDY(n2));      mpz_mod(r,BDY(n1),BDY(n2));
                 if ( !mpz_sgn(r) ) *nr = 0;      if ( !mpz_sgn(r) ) *nr = 0;
                 else MPZTOGZ(r,*nr);      else MPZTOGZ(r,*nr);
         }    }
 }  }
   
 void divqrgz(GZ n1,GZ n2,GZ *nq,GZ *nr)  void divqrgz(GZ n1,GZ n2,GZ *nq,GZ *nr)
 {  {
         mpz_t t, a, b, q, r;    mpz_t t, a, b, q, r;
   
         if ( !n2 ) {    if ( !n2 ) {
                 error("division by 0");      error("division by 0");
                 *nq = 0; *nr = 0;      *nq = 0; *nr = 0;
         } else if ( !n1 ) {    } else if ( !n1 ) {
                 *nq = 0; *nr = 0;      *nq = 0; *nr = 0;
         } else if ( n1 == n2 ) {    } else if ( n1 == n2 ) {
                 mpz_init(t); mpz_set_ui(t,1); MPZTOGZ(t,*nq); *nr = 0;      mpz_init(t); mpz_set_ui(t,1); MPZTOGZ(t,*nq); *nr = 0;
         } else {    } else {
                 mpz_init(q); mpz_init(r);      mpz_init(q); mpz_init(r);
                 mpz_fdiv_qr(q,r,BDY(n1),BDY(n2));      mpz_fdiv_qr(q,r,BDY(n1),BDY(n2));
                 if ( !mpz_sgn(q) ) *nq = 0;      if ( !mpz_sgn(q) ) *nq = 0;
                 else MPZTOGZ(q,*nq);      else MPZTOGZ(q,*nq);
                 if ( !mpz_sgn(r) ) *nr = 0;      if ( !mpz_sgn(r) ) *nr = 0;
                 else MPZTOGZ(r,*nr);      else MPZTOGZ(r,*nr);
         }    }
 }  }
   
 void divsgz(GZ n1,GZ n2,GZ *nq)  void divsgz(GZ n1,GZ n2,GZ *nq)
 {  {
         mpz_t t;    mpz_t t;
         mpq_t a, b, q;    mpq_t a, b, q;
   
         if ( !n2 ) {    if ( !n2 ) {
                 error("division by 0");      error("division by 0");
                 *nq = 0;      *nq = 0;
         } else if ( !n1 )    } else if ( !n1 )
                 *nq = 0;      *nq = 0;
         else if ( n1 == n2 ) {    else if ( n1 == n2 ) {
                 mpz_init(t); mpz_set_ui(t,1); MPZTOGZ(t,*nq);      mpz_init(t); mpz_set_ui(t,1); MPZTOGZ(t,*nq);
         } else {    } else {
                 mpz_init(t); mpz_divexact(t,BDY(n1),BDY(n2)); MPZTOGZ(t,*nq);      mpz_init(t); mpz_divexact(t,BDY(n1),BDY(n2)); MPZTOGZ(t,*nq);
         }    }
 }  }
   
 void chsgngz(GZ n,GZ *nr)  void chsgngz(GZ n,GZ *nr)
 {  {
         mpz_t t;    mpz_t t;
   
         if ( !n )    if ( !n )
                 *nr = 0;      *nr = 0;
         else {    else {
                 t[0] = BDY(n)[0]; mpz_neg(t,t); MPZTOGZ(t,*nr);      t[0] = BDY(n)[0]; mpz_neg(t,t); MPZTOGZ(t,*nr);
         }    }
 }  }
   
 void pwrgz(GZ n1,Q n,GZ *nr)  void pwrgz(GZ n1,Q n,GZ *nr)
 {  {
         mpq_t t,q;    mpq_t t,q;
         mpz_t z;    mpz_t z;
         GQ p,r;    GQ p,r;
   
         if ( !n || UNIGZ(n1) ) *nr = ONEGZ;    if ( !n || UNIGZ(n1) ) *nr = ONEGZ;
         else if ( !n1 ) *nr = 0;    else if ( !n1 ) *nr = 0;
         else if ( !INT(n) ) {    else if ( !INT(n) ) {
                 error("can't calculate fractional power."); *nr = 0;      error("can't calculate fractional power."); *nr = 0;
         } else if ( MUNIGZ(n1) ) *nr = BD(NM(n))[0]%2 ? n1 : ONEGZ;    } else if ( MUNIGZ(n1) ) *nr = BD(NM(n))[0]%2 ? n1 : ONEGZ;
         else if ( PL(NM(n)) > 1 ) {    else if ( PL(NM(n)) > 1 ) {
                 error("exponent too big."); *nr = 0;      error("exponent too big."); *nr = 0;
         } else if ( NID(n1)==N_GZ && SGN(n)>0 ) {    } else if ( NID(n1)==N_GZ && SGN(n)>0 ) {
                 mpz_init(z); mpz_pow_ui(z,BDY(n1),QTOS(n)); MPZTOGZ(z,*nr);      mpz_init(z); mpz_pow_ui(z,BDY(n1),QTOS(n)); MPZTOGZ(z,*nr);
         } else {    } else {
                 MPZTOMPQ(BDY(n1),q); MPQTOGQ(q,r);      MPZTOMPQ(BDY(n1),q); MPQTOGQ(q,r);
                 pwrgq(r,n,&p); *nr = (GZ)p;      pwrgq(r,n,&p); *nr = (GZ)p;
         }    }
 }  }
   
 int cmpgz(GZ q1,GZ q2)  int cmpgz(GZ q1,GZ q2)
 {  {
         int sgn;    int sgn;
   
         if ( !q1 )    if ( !q1 )
                 if ( !q2 )      if ( !q2 )
                         return 0;        return 0;
                 else      else
                         return -mpz_sgn(BDY(q2));        return -mpz_sgn(BDY(q2));
         else if ( !q2 )    else if ( !q2 )
                 return mpz_sgn(BDY(q1));      return mpz_sgn(BDY(q1));
         else if ( (sgn = mpz_sgn(BDY(q1))) != mpz_sgn(BDY(q2)) )    else if ( (sgn = mpz_sgn(BDY(q1))) != mpz_sgn(BDY(q2)) )
                         return sgn;        return sgn;
         else {    else {
                 sgn = mpz_cmp(BDY(q1),BDY(q2));      sgn = mpz_cmp(BDY(q1),BDY(q2));
                 if ( sgn > 0 ) return 1;      if ( sgn > 0 ) return 1;
                 else if ( sgn < 0 ) return -1;      else if ( sgn < 0 ) return -1;
                 else return 0;      else return 0;
         }    }
 }  }
   
 void gcdgz(GZ n1,GZ n2,GZ *nq)  void gcdgz(GZ n1,GZ n2,GZ *nq)
 {  {
         mpz_t t;    mpz_t t;
   
         if ( !n1 ) *nq = n2;    if ( !n1 ) *nq = n2;
         else if ( !n2 ) *nq = n1;    else if ( !n2 ) *nq = n1;
         else {    else {
                 mpz_init(t); mpz_gcd(t,BDY(n1),BDY(n2));      mpz_init(t); mpz_gcd(t,BDY(n1),BDY(n2));
                 MPZTOGZ(t,*nq);      MPZTOGZ(t,*nq);
         }    }
 }  }
   
 void invgz(GZ n1,GZ *nq)  void invgz(GZ n1,GZ *nq)
 {  {
         mpz_t t;    mpz_t t;
   
         mpz_init(t); mpz_invert(t,BDY(n1),BDY(current_mod_lf));    mpz_init(t); mpz_invert(t,BDY(n1),BDY(current_mod_lf));
         MPZTOGZ(t,*nq);    MPZTOGZ(t,*nq);
 }  }
   
 void lcmgz(GZ n1,GZ n2,GZ *nq)  void lcmgz(GZ n1,GZ n2,GZ *nq)
 {  {
         GZ g,t;    GZ g,t;
   
         if ( !n1 || !n2 ) *nq = 0;    if ( !n1 || !n2 ) *nq = 0;
         else {    else {
                 gcdgz(n1,n2,&g); divsgz(n1,g,&t);      gcdgz(n1,n2,&g); divsgz(n1,g,&t);
                 mulgz(n2,t,nq);      mulgz(n2,t,nq);
         }    }
 }  }
   
 void gcdvgz(VECT v,GZ *q)  void gcdvgz(VECT v,GZ *q)
 {  {
         int n,i;    int n,i;
         GZ *b;    GZ *b;
         GZ g,g1;    GZ g,g1;
   
         n = v->len;    n = v->len;
         b = (GZ *)v->body;    b = (GZ *)v->body;
         g = b[0];    g = b[0];
         for ( i = 1; i < n; i++ ) {    for ( i = 1; i < n; i++ ) {
                 gcdgz(g,b[i],&g1); g = g1;      gcdgz(g,b[i],&g1); g = g1;
         }    }
         *q = g;    *q = g;
 }  }
   
 void gcdvgz_estimate(VECT v,GZ *q)  void gcdvgz_estimate(VECT v,GZ *q)
 {  {
         int n,m,i;    int n,m,i;
         GZ s,t,u;    GZ s,t,u;
         GZ *b;    GZ *b;
   
         n = v->len;    n = v->len;
         b = (GZ *)v->body;    b = (GZ *)v->body;
         if ( n == 1 ) {    if ( n == 1 ) {
                 if ( mpz_sgn(BDY(b[0]))<0 ) chsgngz(b[0],q);      if ( mpz_sgn(BDY(b[0]))<0 ) chsgngz(b[0],q);
                 else *q = b[0];      else *q = b[0];
         }    }
         m = n/2;    m = n/2;
         for ( i = 0, s = 0; i < m; i++ ) {    for ( i = 0, s = 0; i < m; i++ ) {
                 if ( b[i] && mpz_sgn(BDY(b[i]))<0 ) subgz(s,b[i],&u);      if ( b[i] && mpz_sgn(BDY(b[i]))<0 ) subgz(s,b[i],&u);
                 else addgz(s,b[i],&u);      else addgz(s,b[i],&u);
                 s = u;      s = u;
         }    }
         for ( i = 0, t = 0; i < n; i++ ) {    for ( i = 0, t = 0; i < n; i++ ) {
                 if ( b[i] && mpz_sgn(BDY(b[i]))<0 ) subgz(t,b[i],&u);      if ( b[i] && mpz_sgn(BDY(b[i]))<0 ) subgz(t,b[i],&u);
                 else addgz(t,b[i],&u);      else addgz(t,b[i],&u);
                 t = u;      t = u;
         }    }
         gcdgz(s,t,q);    gcdgz(s,t,q);
 }  }
   
 void addgq(GQ n1,GQ n2,GQ *nr)  void addgq(GQ n1,GQ n2,GQ *nr)
 {  {
         mpq_t q1,q2,t;    mpq_t q1,q2,t;
   
         if ( !n1 ) *nr = n2;    if ( !n1 ) *nr = n2;
         else if ( !n2 ) *nr = n1;    else if ( !n2 ) *nr = n1;
         else {    else {
                 if ( NID(n1) == N_GZ ) MPZTOMPQ(BDY((GZ)n1),q1);      if ( NID(n1) == N_GZ ) MPZTOMPQ(BDY((GZ)n1),q1);
                 else q1[0] = BDY(n1)[0];      else q1[0] = BDY(n1)[0];
                 if ( NID(n2) == N_GZ ) MPZTOMPQ(BDY((GZ)n2),q2);      if ( NID(n2) == N_GZ ) MPZTOMPQ(BDY((GZ)n2),q2);
                 else q2[0] = BDY(n2)[0];      else q2[0] = BDY(n2)[0];
                 mpq_init(t); mpq_add(t,q1,q2); *nr = mpqtogzq(t);      mpq_init(t); mpq_add(t,q1,q2); *nr = mpqtogzq(t);
         }    }
 }  }
   
 void subgq(GQ n1,GQ n2,GQ *nr)  void subgq(GQ n1,GQ n2,GQ *nr)
 {  {
         mpq_t q1,q2,t;    mpq_t q1,q2,t;
   
         if ( !n1 )    if ( !n1 )
                 if ( !n2 ) *nr = 0;      if ( !n2 ) *nr = 0;
                 else {      else {
                         if ( NID(n1) == N_GZ ) chsgngz((GZ)n1,(GZ *)nr);        if ( NID(n1) == N_GZ ) chsgngz((GZ)n1,(GZ *)nr);
                         else {        else {
                                 mpq_init(t); mpq_neg(t,BDY(n2)); MPQTOGQ(t,*nr);          mpq_init(t); mpq_neg(t,BDY(n2)); MPQTOGQ(t,*nr);
                         }        }
                 }      }
         else if ( !n2 ) *nr = n1;    else if ( !n2 ) *nr = n1;
         else if ( n1 == n2 ) *nr = 0;    else if ( n1 == n2 ) *nr = 0;
         else {    else {
                 if ( NID(n1) == N_GZ ) MPZTOMPQ(BDY((GZ)n1),q1);      if ( NID(n1) == N_GZ ) MPZTOMPQ(BDY((GZ)n1),q1);
                 else q1[0] = BDY(n1)[0];      else q1[0] = BDY(n1)[0];
                 if ( NID(n2) == N_GZ ) MPZTOMPQ(BDY((GZ)n2),q2);      if ( NID(n2) == N_GZ ) MPZTOMPQ(BDY((GZ)n2),q2);
                 else q2[0] = BDY(n2)[0];      else q2[0] = BDY(n2)[0];
                 mpq_init(t); mpq_sub(t,q1,q2); *nr = mpqtogzq(t);      mpq_init(t); mpq_sub(t,q1,q2); *nr = mpqtogzq(t);
         }    }
 }  }
   
 void mulgq(GQ n1,GQ n2,GQ *nr)  void mulgq(GQ n1,GQ n2,GQ *nr)
 {  {
         mpq_t t,q1,q2;    mpq_t t,q1,q2;
   
         if ( !n1 || !n2 ) *nr = 0;    if ( !n1 || !n2 ) *nr = 0;
         else {    else {
                 if ( NID(n1) == N_GZ ) MPZTOMPQ(BDY((GZ)n1),q1);      if ( NID(n1) == N_GZ ) MPZTOMPQ(BDY((GZ)n1),q1);
                 else q1[0] = BDY(n1)[0];      else q1[0] = BDY(n1)[0];
                 if ( NID(n2) == N_GZ ) MPZTOMPQ(BDY((GZ)n2),q2);      if ( NID(n2) == N_GZ ) MPZTOMPQ(BDY((GZ)n2),q2);
                 else q2[0] = BDY(n2)[0];      else q2[0] = BDY(n2)[0];
                 mpq_init(t); mpq_mul(t,q1,q2); *nr = mpqtogzq(t);      mpq_init(t); mpq_mul(t,q1,q2); *nr = mpqtogzq(t);
         }    }
 }  }
   
 void divgq(GQ n1,GQ n2,GQ *nq)  void divgq(GQ n1,GQ n2,GQ *nq)
 {  {
         mpq_t t,q1,q2;    mpq_t t,q1,q2;
   
         if ( !n2 ) {    if ( !n2 ) {
                 error("division by 0");      error("division by 0");
                 *nq = 0;      *nq = 0;
                 return;      return;
         } else if ( !n1 ) *nq = 0;    } else if ( !n1 ) *nq = 0;
         else if ( n1 == n2 ) *nq = (GQ)ONEGZ;    else if ( n1 == n2 ) *nq = (GQ)ONEGZ;
         else {    else {
                 if ( NID(n1) == N_GZ ) MPZTOMPQ(BDY((GZ)n1),q1);      if ( NID(n1) == N_GZ ) MPZTOMPQ(BDY((GZ)n1),q1);
                 else q1[0] = BDY(n1)[0];      else q1[0] = BDY(n1)[0];
                 if ( NID(n2) == N_GZ ) MPZTOMPQ(BDY((GZ)n2),q2);      if ( NID(n2) == N_GZ ) MPZTOMPQ(BDY((GZ)n2),q2);
                 else q2[0] = BDY(n2)[0];      else q2[0] = BDY(n2)[0];
                 mpq_init(t); mpq_div(t,q1,q2); *nq = mpqtogzq(t);      mpq_init(t); mpq_div(t,q1,q2); *nq = mpqtogzq(t);
         }    }
 }  }
   
 void chsgngq(GQ n,GQ *nr)  void chsgngq(GQ n,GQ *nr)
 {  {
         mpq_t t;    mpq_t t;
   
         if ( !n ) *nr = 0;    if ( !n ) *nr = 0;
         else if ( NID(n) == N_GZ ) chsgngz((GZ)n,(GZ *)nr);    else if ( NID(n) == N_GZ ) chsgngz((GZ)n,(GZ *)nr);
         else {    else {
                 mpq_init(t); mpq_neg(t,BDY(n)); MPQTOGQ(t,*nr);      mpq_init(t); mpq_neg(t,BDY(n)); MPQTOGQ(t,*nr);
         }    }
 }  }
   
 void pwrgq(GQ n1,Q n,GQ *nr)  void pwrgq(GQ n1,Q n,GQ *nr)
 {  {
         int e;    int e;
         mpz_t nm,dn;    mpz_t nm,dn;
         mpq_t t;    mpq_t t;
   
         if ( !n || UNIGZ((GZ)n1) || UNIGQ(n1) ) *nr = (GQ)ONEGZ;    if ( !n || UNIGZ((GZ)n1) || UNIGQ(n1) ) *nr = (GQ)ONEGZ;
         else if ( !n1 ) *nr = 0;    else if ( !n1 ) *nr = 0;
         else if ( !INT(n) ) {    else if ( !INT(n) ) {
                 error("can't calculate fractional power."); *nr = 0;      error("can't calculate fractional power."); *nr = 0;
         } else if ( PL(NM(n)) > 1 ) {    } else if ( PL(NM(n)) > 1 ) {
                 error("exponent too big."); *nr = 0;      error("exponent too big."); *nr = 0;
         } else {    } else {
                 e = QTOS(n);      e = QTOS(n);
                 if ( e < 0 ) {      if ( e < 0 ) {
                         e = -e;        e = -e;
                         if ( NID(n1)==N_GZ ) {        if ( NID(n1)==N_GZ ) {
                                 nm[0] = ONEMPZ[0];          nm[0] = ONEMPZ[0];
                                 dn[0] = BDY((GZ)n1)[0];          dn[0] = BDY((GZ)n1)[0];
                         } else {        } else {
                                 nm[0] = mpq_denref(BDY(n1))[0]; dn[0] = mpq_numref(BDY(n1))[0];          nm[0] = mpq_denref(BDY(n1))[0]; dn[0] = mpq_numref(BDY(n1))[0];
                         }        }
                 } else {      } else {
                         if ( NID(n1)==N_GZ ) {        if ( NID(n1)==N_GZ ) {
                                 nm[0] = BDY((GZ)n1)[0]; dn[0] = ONEMPZ[0];          nm[0] = BDY((GZ)n1)[0]; dn[0] = ONEMPZ[0];
                         } else {        } else {
                                 nm[0] = mpq_numref(BDY(n1))[0]; dn[0] = mpq_denref(BDY(n1))[0];          nm[0] = mpq_numref(BDY(n1))[0]; dn[0] = mpq_denref(BDY(n1))[0];
                         }        }
                 }      }
                 mpq_init(t);      mpq_init(t);
                 mpz_pow_ui(mpq_numref(t),nm,e); mpz_pow_ui(mpq_denref(t),dn,e);      mpz_pow_ui(mpq_numref(t),nm,e); mpz_pow_ui(mpq_denref(t),dn,e);
                 *nr = mpqtogzq(t);      *nr = mpqtogzq(t);
         }    }
 }  }
   
 int cmpgq(GQ n1,GQ n2)  int cmpgq(GQ n1,GQ n2)
 {  {
         mpq_t q1,q2;    mpq_t q1,q2;
         int sgn;    int sgn;
   
         if ( !n1 )    if ( !n1 )
                 if ( !n2 ) return 0;      if ( !n2 ) return 0;
                 else return (NID(n2)==N_GZ) ? -mpz_sgn(BDY((GZ)n2)) : -mpq_sgn(BDY(n2));      else return (NID(n2)==N_GZ) ? -mpz_sgn(BDY((GZ)n2)) : -mpq_sgn(BDY(n2));
         if ( !n2 ) return (NID(n1)==N_GZ) ? mpz_sgn(BDY((GZ)n1)) : mpq_sgn(BDY(n1));    if ( !n2 ) return (NID(n1)==N_GZ) ? mpz_sgn(BDY((GZ)n1)) : mpq_sgn(BDY(n1));
         else if ( (sgn = mpq_sgn(BDY(n1))) != mpq_sgn(BDY(n2)) ) return sgn;    else if ( (sgn = mpq_sgn(BDY(n1))) != mpq_sgn(BDY(n2)) ) return sgn;
         else {    else {
                 if ( NID(n1) == N_GZ ) MPZTOMPQ(BDY((GZ)n1),q1);      if ( NID(n1) == N_GZ ) MPZTOMPQ(BDY((GZ)n1),q1);
                 else q1[0] = BDY(n1)[0];      else q1[0] = BDY(n1)[0];
                 if ( NID(n2) == N_GZ ) MPZTOMPQ(BDY((GZ)n2),q2);      if ( NID(n2) == N_GZ ) MPZTOMPQ(BDY((GZ)n2),q2);
                 else q2[0] = BDY(n2)[0];      else q2[0] = BDY(n2)[0];
                 sgn = mpq_cmp(q1,q2);      sgn = mpq_cmp(q1,q2);
                 if ( sgn > 0 ) return 1;      if ( sgn > 0 ) return 1;
                 else if ( sgn < 0 ) return -1;      else if ( sgn < 0 ) return -1;
                 else return 0;      else return 0;
         }    }
 }  }
   
 void mkgwc(int k,int l,GZ *t)  void mkgwc(int k,int l,GZ *t)
 {  {
         mpz_t a,b,q,nm,z,u;    mpz_t a,b,q,nm,z,u;
         int i,n;    int i,n;
   
         n = MIN(k,l);    n = MIN(k,l);
         mpz_init_set_ui(z,1);    mpz_init_set_ui(z,1);
         mpz_init(u); mpz_set(u,z); MPZTOGZ(u,t[0]);    mpz_init(u); mpz_set(u,z); MPZTOGZ(u,t[0]);
         mpz_init(a); mpz_init(b); mpz_init(nm);    mpz_init(a); mpz_init(b); mpz_init(nm);
         for ( i = 1; i <= n; i++ ) {    for ( i = 1; i <= n; i++ ) {
                 mpz_set_ui(a,k-i+1); mpz_set_ui(b,l-i+1); mpz_mul(nm,a,b);      mpz_set_ui(a,k-i+1); mpz_set_ui(b,l-i+1); mpz_mul(nm,a,b);
                 mpz_mul(z,BDY(t[i-1]),nm); mpz_fdiv_q_ui(z,z,i);      mpz_mul(z,BDY(t[i-1]),nm); mpz_fdiv_q_ui(z,z,i);
                 mpz_init(u); mpz_set(u,z); MPZTOGZ(u,t[i]);      mpz_init(u); mpz_set(u,z); MPZTOGZ(u,t[i]);
         }    }
 }  }
   
 void gz_lgp(P p,GZ *g,GZ *l);  void gz_lgp(P p,GZ *g,GZ *l);
   
 void gz_ptozp(P p,int sgn,GQ *c,P *pr)  void gz_ptozp(P p,int sgn,GQ *c,P *pr)
 {  {
         GZ nm,dn;    GZ nm,dn;
   
         if ( !p ) {    if ( !p ) {
                 *c = 0; *pr = 0;      *c = 0; *pr = 0;
         } else {    } else {
                 gz_lgp(p,&nm,&dn);      gz_lgp(p,&nm,&dn);
                 divgz(nm,dn,(GZ *)c);      divgz(nm,dn,(GZ *)c);
                 divsp(CO,p,(P)*c,pr);      divsp(CO,p,(P)*c,pr);
         }    }
 }  }
   
 void gz_lgp(P p,GZ *g,GZ *l)  void gz_lgp(P p,GZ *g,GZ *l)
 {  {
         DCP dc;    DCP dc;
         GZ g1,g2,l1,l2,l3,l4;    GZ g1,g2,l1,l2,l3,l4;
   
         if ( NUM(p) ) {    if ( NUM(p) ) {
                 if ( NID((GZ)p)==N_GZ ) {      if ( NID((GZ)p)==N_GZ ) {
                         MPZTOGZ(BDY((GZ)p),*g);        MPZTOGZ(BDY((GZ)p),*g);
                         *l = ONEGZ;        *l = ONEGZ;
                 } else {      } else {
                         MPZTOGZ(mpq_numref(BDY((GQ)p)),*g);        MPZTOGZ(mpq_numref(BDY((GQ)p)),*g);
                         MPZTOGZ(mpq_denref(BDY((GQ)p)),*l);        MPZTOGZ(mpq_denref(BDY((GQ)p)),*l);
                 }      }
         } else {    } else {
                 dc = DC(p); gz_lgp(COEF(dc),g,l);      dc = DC(p); gz_lgp(COEF(dc),g,l);
                 for ( dc = NEXT(dc); dc; dc = NEXT(dc) ) {      for ( dc = NEXT(dc); dc; dc = NEXT(dc) ) {
                         gz_lgp(COEF(dc),&g1,&l1); gcdgz(*g,g1,&g2); *g = g2;        gz_lgp(COEF(dc),&g1,&l1); gcdgz(*g,g1,&g2); *g = g2;
                         gcdgz(*l,l1,&l2); mulgz(*l,l1,&l3); divgz(l3,l2,l);        gcdgz(*l,l1,&l2); mulgz(*l,l1,&l3); divgz(l3,l2,l);
                 }      }
         }    }
 }  }
   
 void gz_qltozl(GQ *w,int n,GZ *dvr)  void gz_qltozl(GQ *w,int n,GZ *dvr)
 {  {
         GZ nm,dn;    GZ nm,dn;
         GZ g,g1,l1,l2,l3;    GZ g,g1,l1,l2,l3;
         GQ c;    GQ c;
         int i;    int i;
         struct oVECT v;    struct oVECT v;
   
         for ( i = 0; i < n; i++ )    for ( i = 0; i < n; i++ )
                 if ( w[i] && NID(w[i])==N_GQ )      if ( w[i] && NID(w[i])==N_GQ )
                         break;        break;
         if ( i == n ) {    if ( i == n ) {
                 v.id = O_VECT; v.len = n; v.body = (pointer *)w;      v.id = O_VECT; v.len = n; v.body = (pointer *)w;
                 gcdvgz(&v,dvr); return;      gcdvgz(&v,dvr); return;
         }    }
         for ( i = 0; !w[i]; i++ );    for ( i = 0; !w[i]; i++ );
         c = w[i];    c = w[i];
         if ( NID(c)==N_GQ ) {    if ( NID(c)==N_GQ ) {
                 MPZTOGZ(mpq_numref(BDY(c)),nm); MPZTOGZ(mpq_denref(BDY(c)),dn);      MPZTOGZ(mpq_numref(BDY(c)),nm); MPZTOGZ(mpq_denref(BDY(c)),dn);
         } else {    } else {
                 MPZTOGZ(BDY((GZ)c),nm); dn = ONEGZ;      MPZTOGZ(BDY((GZ)c),nm); dn = ONEGZ;
         }    }
         for ( i++; i < n; i++ ) {    for ( i++; i < n; i++ ) {
                 c = w[i];      c = w[i];
                 if ( !c ) continue;      if ( !c ) continue;
                 if ( NID(c)==N_GQ ) {      if ( NID(c)==N_GQ ) {
                         MPZTOGZ(mpq_numref(BDY(c)),g1); MPZTOGZ(mpq_denref(BDY(c)),l1);        MPZTOGZ(mpq_numref(BDY(c)),g1); MPZTOGZ(mpq_denref(BDY(c)),l1);
                 } else {      } else {
                         MPZTOGZ(BDY((GZ)c),g1); l1 = ONEGZ;        MPZTOGZ(BDY((GZ)c),g1); l1 = ONEGZ;
                 }      }
                 gcdgz(nm,g1,&g); nm = g;      gcdgz(nm,g1,&g); nm = g;
                 gcdgz(dn,l1,&l2); mulgz(dn,l1,&l3); divgz(l3,l2,&dn);      gcdgz(dn,l1,&l2); mulgz(dn,l1,&l3); divgz(l3,l2,&dn);
         }    }
         divgz(nm,dn,dvr);    divgz(nm,dn,dvr);
 }  }
   
 int gz_bits(GQ q)  int gz_bits(GQ q)
 {  {
         if ( !q ) return 0;    if ( !q ) return 0;
         else if ( NID(q)==N_Q )    else if ( NID(q)==N_Q )
                 return n_bits(NM((Q)q))+(INT((Q)q)?0:n_bits(DN((Q)q)));      return n_bits(NM((Q)q))+(INT((Q)q)?0:n_bits(DN((Q)q)));
         else if ( NID(q)==N_GZ ) return mpz_sizeinbase(BDY((GZ)q),2);    else if ( NID(q)==N_GZ ) return mpz_sizeinbase(BDY((GZ)q),2);
         else    else
                 return mpz_sizeinbase(mpq_numref(BDY(q)),2)      return mpz_sizeinbase(mpq_numref(BDY(q)),2)
                         + mpz_sizeinbase(mpq_denref(BDY(q)),2);        + mpz_sizeinbase(mpq_denref(BDY(q)),2);
 }  }
   
 int gzp_mag(P p)  int gzp_mag(P p)
 {  {
         int s;    int s;
         DCP dc;    DCP dc;
   
         if ( !p ) return 0;    if ( !p ) return 0;
         else if ( OID(p) == O_N ) return gz_bits((GQ)p);    else if ( OID(p) == O_N ) return gz_bits((GQ)p);
         else {    else {
                 for ( dc = DC(p), s = 0; dc; dc = NEXT(dc) ) s += gzp_mag(COEF(dc));      for ( dc = DC(p), s = 0; dc; dc = NEXT(dc) ) s += gzp_mag(COEF(dc));
                 return s;      return s;
         }    }
 }  }
   
 void makesubstgz(VL v,NODE *s)  void makesubstgz(VL v,NODE *s)
 {  {
         NODE r,r0;    NODE r,r0;
         GZ q;    GZ q;
         unsigned int n;    unsigned int n;
   
         for ( r0 = 0; v; v = NEXT(v) ) {    for ( r0 = 0; v; v = NEXT(v) ) {
                 NEXTNODE(r0,r); BDY(r) = (pointer)v->v;      NEXTNODE(r0,r); BDY(r) = (pointer)v->v;
 #if defined(_PA_RISC1_1)  #if defined(_PA_RISC1_1)
                 n = mrand48()&BMASK; q = utogz(n);      n = mrand48()&BMASK; q = utogz(n);
 #else  #else
                 n = random(); q = utogz(n);      n = random(); q = utogz(n);
 #endif  #endif
                 NEXTNODE(r0,r); BDY(r) = (pointer)q;      NEXTNODE(r0,r); BDY(r) = (pointer)q;
         }    }
         if ( r0 ) NEXT(r) = 0;    if ( r0 ) NEXT(r) = 0;
         *s = r0;    *s = r0;
 }  }
   
 unsigned int remgq(GQ a,unsigned int mod)  unsigned int remgq(GQ a,unsigned int mod)
 {  {
         unsigned int c,nm,dn;    unsigned int c,nm,dn;
         mpz_t r;    mpz_t r;
   
         if ( !a ) return 0;    if ( !a ) return 0;
         else if ( NID(a)==N_GZ ) {    else if ( NID(a)==N_GZ ) {
                 mpz_init(r);      mpz_init(r);
                 c = mpz_fdiv_r_ui(r,BDY((GZ)a),mod);      c = mpz_fdiv_r_ui(r,BDY((GZ)a),mod);
         } else {    } else {
                 mpz_init(r);      mpz_init(r);
                 nm = mpz_fdiv_r_ui(r,mpq_numref(BDY(a)),mod);      nm = mpz_fdiv_r_ui(r,mpq_numref(BDY(a)),mod);
                 dn = mpz_fdiv_r_ui(r,mpq_denref(BDY(a)),mod);      dn = mpz_fdiv_r_ui(r,mpq_denref(BDY(a)),mod);
                 dn = invm(dn,mod);      dn = invm(dn,mod);
                 DMAR(nm,dn,0,mod,c);      DMAR(nm,dn,0,mod,c);
         }    }
         return c;    return c;
 }  }
   
 extern int DP_Print;  extern int DP_Print;
Line 753  extern int DP_Print;
Line 753  extern int DP_Print;
   
 int gz_generic_gauss_elim(MAT mat,MAT *nm,GZ *dn,int **rindp,int **cindp)  int gz_generic_gauss_elim(MAT mat,MAT *nm,GZ *dn,int **rindp,int **cindp)
 {  {
         int **wmat;    int **wmat;
         GZ **bmat,**tmat,*bmi,*tmi;    GZ **bmat,**tmat,*bmi,*tmi;
         GZ q,m1,m2,m3,s,u;    GZ q,m1,m2,m3,s,u;
         int *wmi,*colstat,*wcolstat,*rind,*cind;    int *wmi,*colstat,*wcolstat,*rind,*cind;
         int row,col,ind,md,i,j,k,l,t,t1,rank,rank0,inv;    int row,col,ind,md,i,j,k,l,t,t1,rank,rank0,inv;
         MAT r,crmat;    MAT r,crmat;
         int ret;    int ret;
   
         bmat = (GZ **)mat->body;    bmat = (GZ **)mat->body;
         row = mat->row; col = mat->col;    row = mat->row; col = mat->col;
         wmat = (int **)almat(row,col);    wmat = (int **)almat(row,col);
         colstat = (int *)MALLOC_ATOMIC(col*sizeof(int));    colstat = (int *)MALLOC_ATOMIC(col*sizeof(int));
         wcolstat = (int *)MALLOC_ATOMIC(col*sizeof(int));    wcolstat = (int *)MALLOC_ATOMIC(col*sizeof(int));
         for ( ind = 0; ; ind++ ) {    for ( ind = 0; ; ind++ ) {
                 if ( DP_Print ) {      if ( DP_Print ) {
                         fprintf(asir_out,"."); fflush(asir_out);        fprintf(asir_out,"."); fflush(asir_out);
                 }      }
                 md = get_lprime(ind);      md = get_lprime(ind);
                 for ( i = 0; i < row; i++ )      for ( i = 0; i < row; i++ )
                         for ( j = 0, bmi = bmat[i], wmi = wmat[i]; j < col; j++ )        for ( j = 0, bmi = bmat[i], wmi = wmat[i]; j < col; j++ )
                                 wmi[j] = remgq((GQ)bmi[j],md);          wmi[j] = remgq((GQ)bmi[j],md);
                 rank = generic_gauss_elim_mod(wmat,row,col,md,wcolstat);      rank = generic_gauss_elim_mod(wmat,row,col,md,wcolstat);
                 if ( !ind ) {      if ( !ind ) {
 RESET:  RESET:
                         m1 = utogz(md);        m1 = utogz(md);
                         rank0 = rank;        rank0 = rank;
                         bcopy(wcolstat,colstat,col*sizeof(int));        bcopy(wcolstat,colstat,col*sizeof(int));
                         MKMAT(crmat,rank,col-rank);        MKMAT(crmat,rank,col-rank);
                         MKMAT(r,rank,col-rank); *nm = r;        MKMAT(r,rank,col-rank); *nm = r;
                         tmat = (GZ **)crmat->body;        tmat = (GZ **)crmat->body;
                         for ( i = 0; i < rank; i++ )        for ( i = 0; i < rank; i++ )
                                 for ( j = k = 0, tmi = tmat[i], wmi = wmat[i]; j < col; j++ )          for ( j = k = 0, tmi = tmat[i], wmi = wmat[i]; j < col; j++ )
                                         if ( !colstat[j] ) tmi[k++] = utogz(wmi[j]);            if ( !colstat[j] ) tmi[k++] = utogz(wmi[j]);
                 } else {      } else {
                         if ( rank < rank0 ) {        if ( rank < rank0 ) {
                                 if ( DP_Print ) {          if ( DP_Print ) {
                                         fprintf(asir_out,"lower rank matrix; continuing...\n");            fprintf(asir_out,"lower rank matrix; continuing...\n");
                                         fflush(asir_out);            fflush(asir_out);
                                 }          }
                                 continue;          continue;
                         } else if ( rank > rank0 ) {        } else if ( rank > rank0 ) {
                                 if ( DP_Print ) {          if ( DP_Print ) {
                                         fprintf(asir_out,"higher rank matrix; resetting...\n");            fprintf(asir_out,"higher rank matrix; resetting...\n");
                                         fflush(asir_out);            fflush(asir_out);
                                 }          }
                                 goto RESET;          goto RESET;
                         } else {        } else {
                                 for ( j = 0; (j<col) && (colstat[j]==wcolstat[j]); j++ );          for ( j = 0; (j<col) && (colstat[j]==wcolstat[j]); j++ );
                                 if ( j < col ) {          if ( j < col ) {
                                         if ( DP_Print ) {            if ( DP_Print ) {
                                                 fprintf(asir_out,"inconsitent colstat; resetting...\n");              fprintf(asir_out,"inconsitent colstat; resetting...\n");
                                                 fflush(asir_out);              fflush(asir_out);
                                         }            }
                                         goto RESET;            goto RESET;
                                 }          }
                         }        }
   
                         inv = invm(remgq((GQ)m1,md),md);        inv = invm(remgq((GQ)m1,md),md);
                         m2 = utogz(md); mulgz(m1,m2,&m3);        m2 = utogz(md); mulgz(m1,m2,&m3);
                         for ( i = 0; i < rank; i++ )        for ( i = 0; i < rank; i++ )
                                 for ( j = k = 0, tmi = tmat[i], wmi = wmat[i]; j < col; j++ )          for ( j = k = 0, tmi = tmat[i], wmi = wmat[i]; j < col; j++ )
                                         if ( !colstat[j] ) {            if ( !colstat[j] ) {
                                                 if ( tmi[k] ) {              if ( tmi[k] ) {
                                                 /* f3 = f1+m1*(m1 mod m2)^(-1)*(f2 - f1 mod m2) */              /* f3 = f1+m1*(m1 mod m2)^(-1)*(f2 - f1 mod m2) */
                                                         t = remgq((GQ)tmi[k],md);                t = remgq((GQ)tmi[k],md);
                                                         if ( wmi[j] >= t ) t = wmi[j]-t;                if ( wmi[j] >= t ) t = wmi[j]-t;
                                                         else t = md-(t-wmi[j]);                else t = md-(t-wmi[j]);
                                                         DMAR(t,inv,0,md,t1)                DMAR(t,inv,0,md,t1)
                                                         u = utogz(t1); mulgz(m1,u,&s);                u = utogz(t1); mulgz(m1,u,&s);
                                                         addgz(tmi[k],s,&u); tmi[k] = u;                addgz(tmi[k],s,&u); tmi[k] = u;
                                                 } else if ( wmi[j] ) {              } else if ( wmi[j] ) {
                                                 /* f3 = m1*(m1 mod m2)^(-1)*f2 */              /* f3 = m1*(m1 mod m2)^(-1)*f2 */
                                                         DMAR(wmi[j],inv,0,md,t)                DMAR(wmi[j],inv,0,md,t)
                                                         u = utogz(t); mulgz(m1,u,&s); tmi[k] = s;                u = utogz(t); mulgz(m1,u,&s); tmi[k] = s;
                                                 }              }
                                                 k++;              k++;
                                         }            }
                         m1 = m3;        m1 = m3;
                         if ( ind % GZ_F4_INTRAT_PERIOD )        if ( ind % GZ_F4_INTRAT_PERIOD )
                                 ret = 0;          ret = 0;
                         else        else
                                 ret = gz_intmtoratm(crmat,m1,*nm,dn);          ret = gz_intmtoratm(crmat,m1,*nm,dn);
                         if ( ret ) {        if ( ret ) {
                                 *rindp = rind = (int *)MALLOC_ATOMIC(rank*sizeof(int));          *rindp = rind = (int *)MALLOC_ATOMIC(rank*sizeof(int));
                                 *cindp = cind = (int *)MALLOC_ATOMIC((col-rank)*sizeof(int));          *cindp = cind = (int *)MALLOC_ATOMIC((col-rank)*sizeof(int));
                                 for ( j = k = l = 0; j < col; j++ )          for ( j = k = l = 0; j < col; j++ )
                                         if ( colstat[j] ) rind[k++] = j;            if ( colstat[j] ) rind[k++] = j;
                                         else cind[l++] = j;            else cind[l++] = j;
                                 if ( gz_gensolve_check(mat,*nm,*dn,rind,cind) )          if ( gz_gensolve_check(mat,*nm,*dn,rind,cind) )
                                         return rank;            return rank;
                         }        }
                 }      }
         }    }
 }  }
   
 int gz_generic_gauss_elim2(MAT mat,MAT *nm,GZ *dn,int **rindp,int **cindp)  int gz_generic_gauss_elim2(MAT mat,MAT *nm,GZ *dn,int **rindp,int **cindp)
 {  {
   
         MAT full;    MAT full;
         GZ **bmat,**b;    GZ **bmat,**b;
         GZ *bmi;    GZ *bmi;
         GZ dn0;    GZ dn0;
         int row,col,md,i,j,rank,ret;    int row,col,md,i,j,rank,ret;
         int **wmat;    int **wmat;
         int *wmi;    int *wmi;
         int *colstat,*rowstat;    int *colstat,*rowstat;
   
         bmat = (GZ **)mat->body;    bmat = (GZ **)mat->body;
         row = mat->row; col = mat->col;    row = mat->row; col = mat->col;
         wmat = (int **)almat(row,col);    wmat = (int **)almat(row,col);
         colstat = (int *)MALLOC_ATOMIC(col*sizeof(int));    colstat = (int *)MALLOC_ATOMIC(col*sizeof(int));
         rowstat = (int *)MALLOC_ATOMIC(row*sizeof(int));    rowstat = (int *)MALLOC_ATOMIC(row*sizeof(int));
         /* XXX */    /* XXX */
         md = get_lprime(0);    md = get_lprime(0);
         for ( i = 0; i < row; i++ )    for ( i = 0; i < row; i++ )
                 for ( j = 0, bmi = bmat[i], wmi = wmat[i]; j < col; j++ )      for ( j = 0, bmi = bmat[i], wmi = wmat[i]; j < col; j++ )
                         wmi[j] = remgq((GQ)bmi[j],md);        wmi[j] = remgq((GQ)bmi[j],md);
         rank = generic_gauss_elim_mod2(wmat,row,col,md,colstat,rowstat);    rank = generic_gauss_elim_mod2(wmat,row,col,md,colstat,rowstat);
         b = (GZ **)MALLOC(rank*sizeof(GZ));    b = (GZ **)MALLOC(rank*sizeof(GZ));
         for ( i = 0; i < rank; i++ ) b[i] = bmat[rowstat[i]];    for ( i = 0; i < rank; i++ ) b[i] = bmat[rowstat[i]];
         NEWMAT(full); full->row = rank; full->col = col; full->body = (pointer **)b;    NEWMAT(full); full->row = rank; full->col = col; full->body = (pointer **)b;
         ret = gz_generic_gauss_elim_full(full,nm,dn,rindp,cindp);    ret = gz_generic_gauss_elim_full(full,nm,dn,rindp,cindp);
         if ( !ret ) {    if ( !ret ) {
                 rank = gz_generic_gauss_elim(mat,nm,&dn0,rindp,cindp);      rank = gz_generic_gauss_elim(mat,nm,&dn0,rindp,cindp);
                 for ( i = 0; i < rank; i++ ) dn[i] = dn0;      for ( i = 0; i < rank; i++ ) dn[i] = dn0;
         }    }
         return rank;    return rank;
 }  }
   
 int gz_generic_gauss_elim_full(MAT mat,MAT *nm,GZ *dn,int **rindp,int **cindp)  int gz_generic_gauss_elim_full(MAT mat,MAT *nm,GZ *dn,int **rindp,int **cindp)
 {  {
         int **wmat;    int **wmat;
         int *stat;    int *stat;
         GZ **bmat,**tmat,*bmi,*tmi,*ri;    GZ **bmat,**tmat,*bmi,*tmi,*ri;
         GZ q,m1,m2,m3,s,u;    GZ q,m1,m2,m3,s,u;
         int *wmi,*colstat,*wcolstat,*rind,*cind;    int *wmi,*colstat,*wcolstat,*rind,*cind;
         int row,col,ind,md,i,j,k,l,t,t1,rank,rank0,inv,h;    int row,col,ind,md,i,j,k,l,t,t1,rank,rank0,inv,h;
         MAT r,crmat;    MAT r,crmat;
         int ret,initialized,done;    int ret,initialized,done;
   
         initialized = 0;    initialized = 0;
         bmat = (GZ **)mat->body;    bmat = (GZ **)mat->body;
         row = mat->row; col = mat->col;    row = mat->row; col = mat->col;
         wmat = (int **)almat(row,col);    wmat = (int **)almat(row,col);
         stat = (int *)MALLOC_ATOMIC(row*sizeof(int));    stat = (int *)MALLOC_ATOMIC(row*sizeof(int));
         for ( i = 0; i < row; i++ ) stat[i] = 0;    for ( i = 0; i < row; i++ ) stat[i] = 0;
         colstat = (int *)MALLOC_ATOMIC(col*sizeof(int));    colstat = (int *)MALLOC_ATOMIC(col*sizeof(int));
         wcolstat = (int *)MALLOC_ATOMIC(col*sizeof(int));    wcolstat = (int *)MALLOC_ATOMIC(col*sizeof(int));
         for ( ind = 0; ; ind++ ) {    for ( ind = 0; ; ind++ ) {
                 if ( DP_Print ) {      if ( DP_Print ) {
                         fprintf(asir_out,"."); fflush(asir_out);        fprintf(asir_out,"."); fflush(asir_out);
                 }      }
                 md = get_lprime(ind);      md = get_lprime(ind);
                 for ( i = 0; i < row; i++ )      for ( i = 0; i < row; i++ )
                         for ( j = 0, bmi = bmat[i], wmi = wmat[i]; j < col; j++ )        for ( j = 0, bmi = bmat[i], wmi = wmat[i]; j < col; j++ )
                                 wmi[j] = remgq((GQ)bmi[j],md);          wmi[j] = remgq((GQ)bmi[j],md);
                 rank = generic_gauss_elim_mod(wmat,row,col,md,wcolstat);      rank = generic_gauss_elim_mod(wmat,row,col,md,wcolstat);
                 if ( rank < row ) continue;      if ( rank < row ) continue;
                 if ( !initialized ) {      if ( !initialized ) {
                         m1 = utogz(md);        m1 = utogz(md);
                         bcopy(wcolstat,colstat,col*sizeof(int));        bcopy(wcolstat,colstat,col*sizeof(int));
                         MKMAT(crmat,row,col-row);        MKMAT(crmat,row,col-row);
                         MKMAT(r,row,col-row); *nm = r;        MKMAT(r,row,col-row); *nm = r;
                         tmat = (GZ **)crmat->body;        tmat = (GZ **)crmat->body;
                         for ( i = 0; i < row; i++ )        for ( i = 0; i < row; i++ )
                                 for ( j = k = 0, tmi = tmat[i], wmi = wmat[i]; j < col; j++ )          for ( j = k = 0, tmi = tmat[i], wmi = wmat[i]; j < col; j++ )
                                         if ( !colstat[j] ) tmi[k++] = utogz(wmi[j]);            if ( !colstat[j] ) tmi[k++] = utogz(wmi[j]);
                         initialized = 1;        initialized = 1;
                 } else {      } else {
                         for ( j = 0; (j<col) && (colstat[j]==wcolstat[j]); j++ );        for ( j = 0; (j<col) && (colstat[j]==wcolstat[j]); j++ );
                         if ( j < col ) continue;        if ( j < col ) continue;
   
                         inv = invm(remgq((GQ)m1,md),md);        inv = invm(remgq((GQ)m1,md),md);
                         m2 = utogz(md); mulgz(m1,m2,&m3);        m2 = utogz(md); mulgz(m1,m2,&m3);
                         for ( i = 0; i < row; i++ )        for ( i = 0; i < row; i++ )
                                 switch ( stat[i] ) {          switch ( stat[i] ) {
                                         case 1:            case 1:
                                                 /* consistency check */              /* consistency check */
                                                 ri = (GZ *)BDY(r)[i]; wmi = wmat[i];              ri = (GZ *)BDY(r)[i]; wmi = wmat[i];
                                                 for ( j = 0; j < col; j++ ) if ( colstat[j] ) break;              for ( j = 0; j < col; j++ ) if ( colstat[j] ) break;
                                                 h = md-remgq((GQ)dn[i],md);              h = md-remgq((GQ)dn[i],md);
                                                 for ( j++, k = 0; j < col; j++ )              for ( j++, k = 0; j < col; j++ )
                                                         if ( !colstat[j] ) {                if ( !colstat[j] ) {
                                                                 t = remgq((GQ)ri[k],md);                  t = remgq((GQ)ri[k],md);
                                                                 DMAR(wmi[i],h,t,md,t1);                  DMAR(wmi[i],h,t,md,t1);
                                                                 if ( t1 ) break;                  if ( t1 ) break;
                                                         }                }
                                                 if ( j == col ) { stat[i]++; break; }              if ( j == col ) { stat[i]++; break; }
                                                 else {              else {
                                                         /* fall to the case 0 */                /* fall to the case 0 */
                                                         stat[i] = 0;                stat[i] = 0;
                                                 }              }
                                         case 0:            case 0:
                                                 tmi = tmat[i]; wmi = wmat[i];              tmi = tmat[i]; wmi = wmat[i];
                                                 for ( j = k = 0; j < col; j++ )              for ( j = k = 0; j < col; j++ )
                                                         if ( !colstat[j] ) {                if ( !colstat[j] ) {
                                                                 if ( tmi[k] ) {                  if ( tmi[k] ) {
                                                                 /* f3 = f1+m1*(m1 mod m2)^(-1)*(f2 - f1 mod m2) */                  /* f3 = f1+m1*(m1 mod m2)^(-1)*(f2 - f1 mod m2) */
                                                                         t = remgq((GQ)tmi[k],md);                    t = remgq((GQ)tmi[k],md);
                                                                         if ( wmi[j] >= t ) t = wmi[j]-t;                    if ( wmi[j] >= t ) t = wmi[j]-t;
                                                                         else t = md-(t-wmi[j]);                    else t = md-(t-wmi[j]);
                                                                         DMAR(t,inv,0,md,t1)                    DMAR(t,inv,0,md,t1)
                                                                         u = utogz(t1); mulgz(m1,u,&s);                    u = utogz(t1); mulgz(m1,u,&s);
                                                                         addgz(tmi[k],s,&u); tmi[k] = u;                    addgz(tmi[k],s,&u); tmi[k] = u;
                                                                 } else if ( wmi[j] ) {                  } else if ( wmi[j] ) {
                                                                 /* f3 = m1*(m1 mod m2)^(-1)*f2 */                  /* f3 = m1*(m1 mod m2)^(-1)*f2 */
                                                                         DMAR(wmi[j],inv,0,md,t)                    DMAR(wmi[j],inv,0,md,t)
                                                                         u = utogz(t); mulgz(m1,u,&s); tmi[k] = s;                    u = utogz(t); mulgz(m1,u,&s); tmi[k] = s;
                                                                 }                  }
                                                                 k++;                  k++;
                                                         }                }
                                                 break;              break;
                                         case 2: default:            case 2: default:
                                                 break;              break;
                                 }          }
                         m1 = m3;        m1 = m3;
                         if ( ind % 4 )        if ( ind % 4 )
                                 ret = 0;          ret = 0;
                         else        else
                                 ret = gz_intmtoratm2(crmat,m1,*nm,dn,stat);          ret = gz_intmtoratm2(crmat,m1,*nm,dn,stat);
                         if ( ret ) {        if ( ret ) {
                                 *rindp = rind = (int *)MALLOC_ATOMIC(row*sizeof(int));          *rindp = rind = (int *)MALLOC_ATOMIC(row*sizeof(int));
                                 *cindp = cind = (int *)MALLOC_ATOMIC((col-row)*sizeof(int));          *cindp = cind = (int *)MALLOC_ATOMIC((col-row)*sizeof(int));
                                 for ( j = k = l = 0; j < col; j++ )          for ( j = k = l = 0; j < col; j++ )
                                         if ( colstat[j] ) rind[k++] = j;            if ( colstat[j] ) rind[k++] = j;
                                         else cind[l++] = j;            else cind[l++] = j;
                                 return gz_gensolve_check2(mat,*nm,dn,rind,cind);          return gz_gensolve_check2(mat,*nm,dn,rind,cind);
                         }        }
                 }      }
         }    }
 }  }
   
 int gz_generic_gauss_elim_direct(MAT mat,MAT *nm,GZ *dn,int **rindp,int **cindp){  int gz_generic_gauss_elim_direct(MAT mat,MAT *nm,GZ *dn,int **rindp,int **cindp){
         GZ **bmat,*s;    GZ **bmat,*s;
         GZ u,v,w,x,d,t,y;    GZ u,v,w,x,d,t,y;
         int row,col,i,j,k,l,m,rank;    int row,col,i,j,k,l,m,rank;
         int *colstat,*colpos,*cind;    int *colstat,*colpos,*cind;
         MAT r,in;    MAT r,in;
   
         row = mat->row; col = mat->col;    row = mat->row; col = mat->col;
         MKMAT(in,row,col);    MKMAT(in,row,col);
         for ( i = 0; i < row; i++ )    for ( i = 0; i < row; i++ )
                 for ( j = 0; j < col; j++ ) in->body[i][j] = mat->body[i][j];      for ( j = 0; j < col; j++ ) in->body[i][j] = mat->body[i][j];
         bmat = (GZ **)in->body;    bmat = (GZ **)in->body;
         colstat = (int *)MALLOC_ATOMIC(col*sizeof(int));    colstat = (int *)MALLOC_ATOMIC(col*sizeof(int));
         *rindp = colpos = (int *)MALLOC_ATOMIC(row*sizeof(int));    *rindp = colpos = (int *)MALLOC_ATOMIC(row*sizeof(int));
         for ( j = 0, rank = 0, d = ONEGZ; j < col; j++ ) {    for ( j = 0, rank = 0, d = ONEGZ; j < col; j++ ) {
                 for ( i = rank; i < row && !bmat[i][j]; i++  );      for ( i = rank; i < row && !bmat[i][j]; i++  );
                 if ( i == row ) { colstat[j] = 0; continue; }      if ( i == row ) { colstat[j] = 0; continue; }
                 else { colstat[j] = 1; colpos[rank] = j; }      else { colstat[j] = 1; colpos[rank] = j; }
                 if ( i != rank )      if ( i != rank )
                         for ( k = j; k < col; k++ ) {        for ( k = j; k < col; k++ ) {
                                 t = bmat[i][k]; bmat[i][k] = bmat[rank][k]; bmat[rank][k] = t;          t = bmat[i][k]; bmat[i][k] = bmat[rank][k]; bmat[rank][k] = t;
                         }        }
                 for ( i = rank+1, v = bmat[rank][j]; i < row; i++ )      for ( i = rank+1, v = bmat[rank][j]; i < row; i++ )
                         for ( k = j, u = bmat[i][j]; k < col; k++ ) {        for ( k = j, u = bmat[i][j]; k < col; k++ ) {
                                 mulgz(bmat[i][k],v,&w); mulgz(bmat[rank][k],u,&x);          mulgz(bmat[i][k],v,&w); mulgz(bmat[rank][k],u,&x);
                                 subgz(w,x,&y); divsgz(y,d,&bmat[i][k]);          subgz(w,x,&y); divsgz(y,d,&bmat[i][k]);
                         }        }
                 d = v; rank++;      d = v; rank++;
         }    }
         *dn = d;    *dn = d;
         s = (GZ *)MALLOC(col*sizeof(GZ));    s = (GZ *)MALLOC(col*sizeof(GZ));
         for ( i = rank-1; i >= 0; i-- ) {    for ( i = rank-1; i >= 0; i-- ) {
                 for ( k = colpos[i]; k < col; k++ ) mulgz(bmat[i][k],d,&s[k]);      for ( k = colpos[i]; k < col; k++ ) mulgz(bmat[i][k],d,&s[k]);
                 for ( m = rank-1; m > i; m-- ) {      for ( m = rank-1; m > i; m-- ) {
                         for ( k = colpos[m], u = bmat[i][k]; k < col; k++ ) {        for ( k = colpos[m], u = bmat[i][k]; k < col; k++ ) {
                                 mulgz(bmat[m][k],u,&w); subgz(s[k],w,&x); s[k] = x;          mulgz(bmat[m][k],u,&w); subgz(s[k],w,&x); s[k] = x;
                         }        }
                 }      }
                 for ( k = colpos[i], u = bmat[i][k]; k < col; k++ )      for ( k = colpos[i], u = bmat[i][k]; k < col; k++ )
                         divgz(s[k],u,&bmat[i][k]);        divgz(s[k],u,&bmat[i][k]);
         }    }
         *cindp = cind = (int *)MALLOC_ATOMIC((col-rank)*sizeof(int));    *cindp = cind = (int *)MALLOC_ATOMIC((col-rank)*sizeof(int));
         MKMAT(r,rank,col-rank); *nm = r;    MKMAT(r,rank,col-rank); *nm = r;
         for ( j = 0, k = 0; j < col; j++ )    for ( j = 0, k = 0; j < col; j++ )
                 if ( !colstat[j] ) {      if ( !colstat[j] ) {
                         cind[k] = j;        cind[k] = j;
                         for ( i = 0; i < rank; i++ ) r->body[i][k] = bmat[i][j];        for ( i = 0; i < rank; i++ ) r->body[i][k] = bmat[i][j];
                         k++;        k++;
                 }      }
         return rank;    return rank;
 }  }
   
 int gz_intmtoratm(MAT mat,GZ md,MAT nm,GZ *dn)  int gz_intmtoratm(MAT mat,GZ md,MAT nm,GZ *dn)
 {  {
         GZ t,s,b,dn0,dn1,nm1,q,u,unm,udn,dmy;    GZ t,s,b,dn0,dn1,nm1,q,u,unm,udn,dmy;
         int i,j,k,l,row,col,sgn,ret;    int i,j,k,l,row,col,sgn,ret;
         GZ **rmat,**tmat,*tmi,*nmk;    GZ **rmat,**tmat,*tmi,*nmk;
   
         if ( UNIGZ(md) )    if ( UNIGZ(md) )
                 return 0;      return 0;
         row = mat->row; col = mat->col;    row = mat->row; col = mat->col;
         bshiftgz(md,1,&t);    bshiftgz(md,1,&t);
         isqrtgz(t,&s);    isqrtgz(t,&s);
         bshiftgz(s,64,&b);    bshiftgz(s,64,&b);
         if ( !b ) b = ONEGZ;    if ( !b ) b = ONEGZ;
         dn0 = ONEGZ;    dn0 = ONEGZ;
         tmat = (GZ **)mat->body;    tmat = (GZ **)mat->body;
         rmat = (GZ **)nm->body;    rmat = (GZ **)nm->body;
         for ( i = 0; i < row; i++ )    for ( i = 0; i < row; i++ )
                 for ( j = 0, tmi = tmat[i]; j < col; j++ )      for ( j = 0, tmi = tmat[i]; j < col; j++ )
                         if ( tmi[j] ) {        if ( tmi[j] ) {
                                 mulgz(tmi[j],dn0,&s);          mulgz(tmi[j],dn0,&s);
                                 divqrgz(s,md,&dmy,&u);          divqrgz(s,md,&dmy,&u);
                                 ret = gz_inttorat(u,md,b,&sgn,&unm,&udn);          ret = gz_inttorat(u,md,b,&sgn,&unm,&udn);
                                 if ( !ret ) return 0;          if ( !ret ) return 0;
                                 else {          else {
                                         if ( sgn < 0 ) chsgngz(unm,&nm1);            if ( sgn < 0 ) chsgngz(unm,&nm1);
                                         else nm1 = unm;            else nm1 = unm;
                                         dn1 = udn;            dn1 = udn;
                                         if ( !UNIGZ(dn1) ) {            if ( !UNIGZ(dn1) ) {
                                                 for ( k = 0; k < i; k++ )              for ( k = 0; k < i; k++ )
                                                         for ( l = 0, nmk = rmat[k]; l < col; l++ ) {                for ( l = 0, nmk = rmat[k]; l < col; l++ ) {
                                                                 mulgz(nmk[l],dn1,&q); nmk[l] = q;                  mulgz(nmk[l],dn1,&q); nmk[l] = q;
                                                         }                }
                                                 for ( l = 0, nmk = rmat[i]; l < j; l++ ) {              for ( l = 0, nmk = rmat[i]; l < j; l++ ) {
                                                         mulgz(nmk[l],dn1,&q); nmk[l] = q;                mulgz(nmk[l],dn1,&q); nmk[l] = q;
                                                 }              }
                                         }            }
                                         rmat[i][j] = nm1;            rmat[i][j] = nm1;
                                         mulgz(dn0,dn1,&q); dn0 = q;            mulgz(dn0,dn1,&q); dn0 = q;
                                 }          }
                         }        }
         *dn = dn0;    *dn = dn0;
         return 1;    return 1;
 }  }
   
 int gz_intmtoratm2(MAT mat,GZ md,MAT nm,GZ *dn,int *stat)  int gz_intmtoratm2(MAT mat,GZ md,MAT nm,GZ *dn,int *stat)
 {  {
         int row,col,i,j,ret;    int row,col,i,j,ret;
         GZ dn0,dn1,t,s,b;    GZ dn0,dn1,t,s,b;
         GZ *w,*tmi;    GZ *w,*tmi;
         GZ **tmat;    GZ **tmat;
   
         bshiftgz(md,1,&t);    bshiftgz(md,1,&t);
         isqrtgz(t,&s);    isqrtgz(t,&s);
         bshiftgz(s,64,&b);    bshiftgz(s,64,&b);
         tmat = (GZ **)mat->body;    tmat = (GZ **)mat->body;
         if ( UNIGZ(md) ) return 0;    if ( UNIGZ(md) ) return 0;
         row = mat->row; col = mat->col;    row = mat->row; col = mat->col;
         dn0 = ONEGZ;    dn0 = ONEGZ;
         for ( i = 0; i < row; i++ )    for ( i = 0; i < row; i++ )
                 if ( cmpgz(dn[i],dn0) > 0 ) dn0 = dn[i];      if ( cmpgz(dn[i],dn0) > 0 ) dn0 = dn[i];
         w = (GZ *)MALLOC(col*sizeof(GZ));    w = (GZ *)MALLOC(col*sizeof(GZ));
         for ( i = 0; i < row; i++ )    for ( i = 0; i < row; i++ )
                 if ( stat[i] == 0 ) {      if ( stat[i] == 0 ) {
                         for ( j = 0, tmi = tmat[i]; j < col; j++ )        for ( j = 0, tmi = tmat[i]; j < col; j++ )
                                         mulgz(tmi[j],dn0,&w[j]);            mulgz(tmi[j],dn0,&w[j]);
                         ret = gz_intvtoratv(w,col,md,b,BDY(nm)[i],&dn[i]);        ret = gz_intvtoratv(w,col,md,b,BDY(nm)[i],&dn[i]);
                         if ( ret ) {        if ( ret ) {
                                 stat[i] = 1;          stat[i] = 1;
                                 mulgz(dn0,dn[i],&t); dn[i] = t; dn0 = t;          mulgz(dn0,dn[i],&t); dn[i] = t; dn0 = t;
                         }        }
                 }      }
         for ( i = 0; i < row; i++ ) if ( !stat[i] ) break;    for ( i = 0; i < row; i++ ) if ( !stat[i] ) break;
         if ( i == row ) return 1;    if ( i == row ) return 1;
         else return 0;    else return 0;
 }  }
   
 int gz_intvtoratv(GZ *v,int n,GZ md,GZ b,GZ *nm,GZ *dn)  int gz_intvtoratv(GZ *v,int n,GZ md,GZ b,GZ *nm,GZ *dn)
 {  {
         GZ dn0,dn1,q,s,u,nm1,unm,udn,dmy;    GZ dn0,dn1,q,s,u,nm1,unm,udn,dmy;
         GZ *nmk;    GZ *nmk;
         int j,l,col,ret,sgn;    int j,l,col,ret,sgn;
   
         for ( j = 0; j < n; j++ ) nm[j] = 0;    for ( j = 0; j < n; j++ ) nm[j] = 0;
         dn0 = ONEGZ;    dn0 = ONEGZ;
         for ( j = 0; j < n; j++ ) {    for ( j = 0; j < n; j++ ) {
                 if ( !v[j] ) continue;      if ( !v[j] ) continue;
                 mulgz(v[j],dn0,&s);      mulgz(v[j],dn0,&s);
                 divqrgz(s,md,&dmy,&u);      divqrgz(s,md,&dmy,&u);
                 ret = gz_inttorat(u,md,b,&sgn,&unm,&udn);      ret = gz_inttorat(u,md,b,&sgn,&unm,&udn);
                 if ( !ret ) return 0;      if ( !ret ) return 0;
                 if ( sgn < 0 ) chsgngz(unm,&nm1);      if ( sgn < 0 ) chsgngz(unm,&nm1);
                 else nm1 = unm;      else nm1 = unm;
                 dn1 = udn;      dn1 = udn;
                 if ( !UNIGZ(dn1) )      if ( !UNIGZ(dn1) )
                         for ( l = 0; l < j; l++ ) {        for ( l = 0; l < j; l++ ) {
                                 mulgz(nm[l],dn1,&q); nm[l] = q;          mulgz(nm[l],dn1,&q); nm[l] = q;
                         }        }
                 nm[j] = nm1;      nm[j] = nm1;
                 mulgz(dn0,dn1,&q); dn0 = q;      mulgz(dn0,dn1,&q); dn0 = q;
         }    }
         *dn = dn0;    *dn = dn0;
         return 1;    return 1;
 }  }
   
 /* assuming 0 < c < m */  /* assuming 0 < c < m */
   
 int gz_inttorat(GZ c,GZ m,GZ b,int *sgnp,GZ *nmp,GZ *dnp)  int gz_inttorat(GZ c,GZ m,GZ b,int *sgnp,GZ *nmp,GZ *dnp)
 {  {
         GZ qq,t,u1,v1,r1;    GZ qq,t,u1,v1,r1;
         GZ q,u2,v2,r2;    GZ q,u2,v2,r2;
   
         u1 = 0; v1 = ONEGZ; u2 = m; v2 = c;    u1 = 0; v1 = ONEGZ; u2 = m; v2 = c;
         while ( cmpgz(v2,b) >= 0 ) {    while ( cmpgz(v2,b) >= 0 ) {
                 divqrgz(u2,v2,&q,&r2); u2 = v2; v2 = r2;      divqrgz(u2,v2,&q,&r2); u2 = v2; v2 = r2;
                 mulgz(q,v1,&t); subgz(u1,t,&r1); u1 = v1; v1 = r1;      mulgz(q,v1,&t); subgz(u1,t,&r1); u1 = v1; v1 = r1;
         }    }
         if ( cmpgz(v1,b) >= 0 ) return 0;    if ( cmpgz(v1,b) >= 0 ) return 0;
         else {    else {
                 *nmp = v2;      *nmp = v2;
                 if ( mpz_sgn(BDY(v1))<0  ) {      if ( mpz_sgn(BDY(v1))<0  ) {
                         *sgnp = -1; chsgngz(v1,dnp);        *sgnp = -1; chsgngz(v1,dnp);
                 } else {      } else {
                         *sgnp = 1; *dnp = v1;        *sgnp = 1; *dnp = v1;
                 }      }
                 return 1;      return 1;
         }    }
 }  }
   
 extern int f4_nocheck;  extern int f4_nocheck;
   
 int gz_gensolve_check(MAT mat,MAT nm,GZ dn,int *rind,int *cind)  int gz_gensolve_check(MAT mat,MAT nm,GZ dn,int *rind,int *cind)
 {  {
         int row,col,rank,clen,i,j,k,l;    int row,col,rank,clen,i,j,k,l;
         GZ s,t;    GZ s,t;
         GZ *w;    GZ *w;
         GZ *mati,*nmk;    GZ *mati,*nmk;
   
         if ( f4_nocheck ) return 1;    if ( f4_nocheck ) return 1;
         row = mat->row; col = mat->col; rank = nm->row; clen = nm->col;    row = mat->row; col = mat->col; rank = nm->row; clen = nm->col;
         w = (GZ *)MALLOC(clen*sizeof(GZ));    w = (GZ *)MALLOC(clen*sizeof(GZ));
         for ( i = 0; i < row; i++ ) {    for ( i = 0; i < row; i++ ) {
                 mati = (GZ *)mat->body[i];      mati = (GZ *)mat->body[i];
                 bzero(w,clen*sizeof(GZ));      bzero(w,clen*sizeof(GZ));
                 for ( k = 0; k < rank; k++ )      for ( k = 0; k < rank; k++ )
                         for ( l = 0, nmk = (GZ *)nm->body[k]; l < clen; l++ ) {        for ( l = 0, nmk = (GZ *)nm->body[k]; l < clen; l++ ) {
                                 mulgz(mati[rind[k]],nmk[l],&t); addgz(w[l],t,&s); w[l] = s;          mulgz(mati[rind[k]],nmk[l],&t); addgz(w[l],t,&s); w[l] = s;
                         }        }
                 for ( j = 0; j < clen; j++ ) {      for ( j = 0; j < clen; j++ ) {
                         mulgz(dn,mati[cind[j]],&t);        mulgz(dn,mati[cind[j]],&t);
                         if ( cmpgz(w[j],t) ) break;        if ( cmpgz(w[j],t) ) break;
                 }      }
                 if ( j != clen ) break;      if ( j != clen ) break;
         }    }
         if ( i != row ) return 0;    if ( i != row ) return 0;
         else return 1;    else return 1;
 }  }
   
 int gz_gensolve_check2(MAT mat,MAT nm,GZ *dn,int *rind,int *cind)  int gz_gensolve_check2(MAT mat,MAT nm,GZ *dn,int *rind,int *cind)
 {  {
         int row,col,rank,clen,i,j,k,l;    int row,col,rank,clen,i,j,k,l;
         GZ s,t,u,d;    GZ s,t,u,d;
         GZ *w,*m;    GZ *w,*m;
         GZ *mati,*nmk;    GZ *mati,*nmk;
   
         if ( f4_nocheck ) return 1;    if ( f4_nocheck ) return 1;
         row = mat->row; col = mat->col; rank = nm->row; clen = nm->col;    row = mat->row; col = mat->col; rank = nm->row; clen = nm->col;
         w = (GZ *)MALLOC(clen*sizeof(GZ));    w = (GZ *)MALLOC(clen*sizeof(GZ));
         m = (GZ *)MALLOC(clen*sizeof(GZ));    m = (GZ *)MALLOC(clen*sizeof(GZ));
         for ( d = dn[0], i = 1; i < rank; i++ ) {    for ( d = dn[0], i = 1; i < rank; i++ ) {
                 lcmgz(d,dn[i],&t); d = t;      lcmgz(d,dn[i],&t); d = t;
         }    }
         for ( i = 0; i < rank; i++ ) divsgz(d,dn[i],&m[i]);    for ( i = 0; i < rank; i++ ) divsgz(d,dn[i],&m[i]);
         for ( i = 0; i < row; i++ ) {    for ( i = 0; i < row; i++ ) {
                 mati = (GZ *)mat->body[i];      mati = (GZ *)mat->body[i];
                 bzero(w,clen*sizeof(GZ));      bzero(w,clen*sizeof(GZ));
                 for ( k = 0; k < rank; k++ ) {      for ( k = 0; k < rank; k++ ) {
                         mulgz(mati[rind[k]],m[k],&u);        mulgz(mati[rind[k]],m[k],&u);
                         for ( l = 0, nmk = (GZ *)nm->body[k]; l < clen; l++ ) {        for ( l = 0, nmk = (GZ *)nm->body[k]; l < clen; l++ ) {
                                 mulgz(u,nmk[l],&t); addgz(w[l],t,&s); w[l] = s;          mulgz(u,nmk[l],&t); addgz(w[l],t,&s); w[l] = s;
                         }        }
                 }      }
                 for ( j = 0; j < clen; j++ ) {      for ( j = 0; j < clen; j++ ) {
                         mulgz(d,mati[cind[j]],&t);        mulgz(d,mati[cind[j]],&t);
                         if ( cmpgz(w[j],t) ) break;        if ( cmpgz(w[j],t) ) break;
                 }      }
                 if ( j != clen ) break;      if ( j != clen ) break;
         }    }
         if ( i != row ) return 0;    if ( i != row ) return 0;
         else return 1;    else return 1;
 }  }
   
 void isqrtgz(GZ a,GZ *r)  void isqrtgz(GZ a,GZ *r)
 {  {
         int k;    int k;
         GZ x,t,x2,xh,quo,rem;    GZ x,t,x2,xh,quo,rem;
         Q two;    Q two;
   
         if ( !a ) *r = 0;    if ( !a ) *r = 0;
         else if ( UNIGZ(a) ) *r = ONEGZ;    else if ( UNIGZ(a) ) *r = ONEGZ;
         else {    else {
                 k = gz_bits((GQ)a); /* a <= 2^k-1 */      k = gz_bits((GQ)a); /* a <= 2^k-1 */
                 bshiftgz(ONEGZ,-((k>>1)+(k&1)),&x); /* a <= x^2 */      bshiftgz(ONEGZ,-((k>>1)+(k&1)),&x); /* a <= x^2 */
                 STOQ(2,two);      STOQ(2,two);
                 while ( 1 ) {      while ( 1 ) {
                         pwrgz(x,two,&t);        pwrgz(x,two,&t);
                         if ( cmpgz(t,a) <= 0 ) {        if ( cmpgz(t,a) <= 0 ) {
                                 *r = x; return;          *r = x; return;
                         } else {        } else {
                                 if ( mpz_tstbit(BDY(x),0) ) addgz(x,a,&t);          if ( mpz_tstbit(BDY(x),0) ) addgz(x,a,&t);
                                 else t = a;          else t = a;
                                 bshiftgz(x,-1,&x2); divqrgz(t,x2,&quo,&rem);          bshiftgz(x,-1,&x2); divqrgz(t,x2,&quo,&rem);
                                 bshiftgz(x,1,&xh); addgz(quo,xh,&x);          bshiftgz(x,1,&xh); addgz(quo,xh,&x);
                         }        }
                 }      }
         }    }
 }  }
   
 void bshiftgz(GZ a,int n,GZ *r)  void bshiftgz(GZ a,int n,GZ *r)
 {  {
         mpz_t t;    mpz_t t;
   
         if ( !a ) *r = 0;    if ( !a ) *r = 0;
         else if ( n == 0 ) *r = a;    else if ( n == 0 ) *r = a;
         else if ( n < 0 ) {    else if ( n < 0 ) {
                 mpz_init(t); mpz_mul_2exp(t,BDY(a),-n); MPZTOGZ(t,*r);      mpz_init(t); mpz_mul_2exp(t,BDY(a),-n); MPZTOGZ(t,*r);
         } else {    } else {
                 mpz_init(t); mpz_fdiv_q_2exp(t,BDY(a),n);      mpz_init(t); mpz_fdiv_q_2exp(t,BDY(a),n);
                 if ( !mpz_sgn(t) ) *r = 0;      if ( !mpz_sgn(t) ) *r = 0;
                 else MPZTOGZ(t,*r);      else MPZTOGZ(t,*r);
         }    }
 }  }
   
 void addlf(GZ a,GZ b,GZ *c)  void addlf(GZ a,GZ b,GZ *c)
Legend:
Removed from v.1.9  
changed lines
  Added in v.1.10
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>