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

Diff for /OpenXM_contrib2/asir2018/engine/Q.c between version 1.1 and 1.10

version 1.1, 2018/09/19 05:45:07 version 1.10, 2018/12/10 22:24:42
Line 1 
Line 1 
   /* $OpenXM: OpenXM_contrib2/asir2018/engine/Q.c,v 1.9 2018/10/19 23:27:38 noro Exp $ */
 #include "ca.h"  #include "ca.h"
 #include "gmp.h"  #include "gmp.h"
 #include "base.h"  #include "base.h"
Line 10  Z current_mod_lf;
Line 11  Z current_mod_lf;
 int current_mod_lf_size;  int current_mod_lf_size;
 gmp_randstate_t GMP_RAND;  gmp_randstate_t GMP_RAND;
   
   #define F4_INTRAT_PERIOD 4
   
   extern int DP_Print;
   
 void isqrtz(Z a,Z *r);  void isqrtz(Z a,Z *r);
 void bshiftz(Z a,int n,Z *r);  void bshiftz(Z a,int n,Z *r);
   
Line 25  void gc_free(void *p,size_t size)
Line 30  void gc_free(void *p,size_t size)
   
 void init_gmpq()  void init_gmpq()
 {  {
   mp_set_memory_functions(Risa_GC_malloc_atomic,gc_realloc,gc_free);    mp_set_memory_functions(Risa_GC_malloc,gc_realloc,gc_free);
   
   mpz_init(ONEMPZ); mpz_set_ui(ONEMPZ,1); MPZTOZ(ONEMPZ,ONE);    mpz_init(ONEMPZ); mpz_set_ui(ONEMPZ,1); MPZTOZ(ONEMPZ,ONE);
   gmp_randinit_default(GMP_RAND);    gmp_randinit_default(GMP_RAND);
 }  }
   
   void printexpr(VL,Obj);
   
   void pmat(Z **a,int row,int col)
   {
     int i,j;
   
     for ( i = 0; i < row; i++, printf("\n") )
       for ( j = 0; j < col; j++, printf(" ") )
         printexpr(CO,(Obj)a[i][j]);
     printf("\n");
   }
   
 Z utoz(unsigned int u)  Z utoz(unsigned int u)
 {  {
   mpz_t z;    mpz_t z;
Line 161  void mulz(Z n1,Z n2,Z *nr)
Line 178  void mulz(Z n1,Z n2,Z *nr)
   
 void muladdtoz(Z n1,Z n2,Z *nr)  void muladdtoz(Z n1,Z n2,Z *nr)
 {  {
   #if 0
   Z t;    Z t;
   
   if ( n1 && n2 ) {    if ( n1 && n2 ) {
Line 168  void muladdtoz(Z n1,Z n2,Z *nr)
Line 186  void muladdtoz(Z n1,Z n2,Z *nr)
           NEWZ(t); mpz_init(BDY(t)); *nr = t;            NEWZ(t); mpz_init(BDY(t)); *nr = t;
         }          }
         mpz_addmul(BDY(*nr),BDY(n1),BDY(n2));          mpz_addmul(BDY(*nr),BDY(n1),BDY(n2));
     }          if ( !mpz_sgn(BDY(*nr)) )
             *nr = 0;
     }
   #else
     Z t,s;
   
     mulz(n1,n2,&t); addz(*nr,t,&s); *nr = s;
   #endif
 }  }
   
 /* nr += n1*u */  /* nr += n1*u */
   
 void mul1addtoz(Z n1,long u,Z *nr)  void mul1addtoz(Z n1,long u,Z *nr)
 {  {
   #if 0
   Z t;    Z t;
   
   if ( n1 && u ) {    if ( n1 && u ) {
Line 185  void mul1addtoz(Z n1,long u,Z *nr)
Line 211  void mul1addtoz(Z n1,long u,Z *nr)
           mpz_addmul_ui(BDY(*nr),BDY(n1),(unsigned long)u);            mpz_addmul_ui(BDY(*nr),BDY(n1),(unsigned long)u);
         else          else
           mpz_submul_ui(BDY(*nr),BDY(n1),(unsigned long)(-u));            mpz_submul_ui(BDY(*nr),BDY(n1),(unsigned long)(-u));
           if ( !mpz_sgn(BDY(*nr)) )
             *nr = 0;
     }      }
   #else
     Z t,s;
   
     mul1z(n1,u,&t); addz(*nr,t,&s); *nr = s;
   #endif
 }  }
   
 void mul1z(Z n1,long n2,Z *nr)  void mul1z(Z n1,long n2,Z *nr)
Line 327  void pwrz(Z n1,Z n,Z *nr)
Line 360  void pwrz(Z n1,Z n,Z *nr)
   } else if ( !smallz(n) ) {    } else if ( !smallz(n) ) {
     error("exponent too big."); *nr = 0;      error("exponent too big."); *nr = 0;
   } else if ( n1->z && mpz_sgn(BDY((Z)n))>0 ) {    } else if ( n1->z && mpz_sgn(BDY((Z)n))>0 ) {
     mpz_init(z); mpz_pow_ui(z,BDY(n1),QTOS(n)); MPZTOZ(z,*nr);      mpz_init(z); mpz_pow_ui(z,BDY(n1),ZTOS(n)); MPZTOZ(z,*nr);
   } else {    } else {
     MPZTOMPQ(BDY(n1),q); MPQTOQ(q,r);      MPZTOMPQ(BDY(n1),q); MPQTOQ(q,r);
     pwrq(r,(Q)n,&p); *nr = (Z)p;      pwrq(r,(Q)n,&p); *nr = (Z)p;
Line 427  void gcdvz_estimate(VECT v,Z *q)
Line 460  void gcdvz_estimate(VECT v,Z *q)
     else addz(s,b[i],&u);      else addz(s,b[i],&u);
     s = u;      s = u;
   }    }
   for ( i = 0, t = 0; i < n; i++ ) {    for ( t = 0; i < n; i++ ) {
     if ( b[i] && mpz_sgn(BDY(b[i]))<0 ) subz(t,b[i],&u);      if ( b[i] && mpz_sgn(BDY(b[i]))<0 ) subz(t,b[i],&u);
     else addz(t,b[i],&u);      else addz(t,b[i],&u);
     t = u;      t = u;
Line 435  void gcdvz_estimate(VECT v,Z *q)
Line 468  void gcdvz_estimate(VECT v,Z *q)
   gcdz(s,t,q);    gcdz(s,t,q);
 }  }
   
   void gcdv_mpz_estimate(mpz_t g,mpz_t *b,int n)
   {
     int m,m2,i,j;
     mpz_t s,t;
   
     mpz_init(g);
     for ( i = 0, m = 0; i < n; i++ )
       if ( mpz_sgn(b[i]) ) m++;
     if ( !m ) {
       mpz_set_ui(g,0);
       return;
     }
     if ( m == 1 ) {
       for ( i = 0, m = 0; i < n; i++ )
         if ( mpz_sgn(b[i]) ) break;
       if ( mpz_sgn(b[i])<0 ) mpz_neg(g,b[i]);
       else mpz_set(g,b[i]);
       return ;
     }
     m2 = m/2;
     mpz_init_set_ui(s,0);
     for ( i = j = 0; j < m2; i++ ) {
       if ( mpz_sgn(b[i]) ) {
         if ( mpz_sgn(b[i])<0 )
           mpz_sub(s,s,b[i]);
         else
           mpz_add(s,s,b[i]);
         j++;
       }
     }
     mpz_init_set_ui(t,0);
     for ( ; i < n; i++ ) {
       if ( mpz_sgn(b[i]) ) {
         if ( mpz_sgn(b[i])<0 )
           mpz_sub(t,t,b[i]);
         else
           mpz_add(t,t,b[i]);
       }
     }
     mpz_gcd(g,s,t);
   }
   
   
 void factorialz(unsigned int n,Z *nr)  void factorialz(unsigned int n,Z *nr)
 {  {
   mpz_t a;    mpz_t a;
Line 579  void pwrq(Q n1,Q n,Q *nr)
Line 655  void pwrq(Q n1,Q n,Q *nr)
   } else if ( !smallz((Z)n) ) {    } else if ( !smallz((Z)n) ) {
     error("exponent too big."); *nr = 0;      error("exponent too big."); *nr = 0;
   } else {    } else {
     e = QTOS(n);      e = ZTOS(n);
     if ( e < 0 ) {      if ( e < 0 ) {
       e = -e;        e = -e;
       if ( n1->z ) {        if ( n1->z ) {
Line 636  void mkbc(int n,Z *t)
Line 712  void mkbc(int n,Z *t)
   Z c,d,iq;    Z c,d,iq;
   
   for ( t[0] = ONE, i = 1; i <= n/2; i++ ) {    for ( t[0] = ONE, i = 1; i <= n/2; i++ ) {
     STOQ(n-i+1,c); mulz(t[i-1],c,&d);      STOZ(n-i+1,c); mulz(t[i-1],c,&d);
     STOQ(i,iq); divsz(d,iq,&t[i]);      STOZ(i,iq); divsz(d,iq,&t[i]);
   }    }
   for ( ; i <= n; i++ )    for ( ; i <= n; i++ )
     t[i] = t[n-i];      t[i] = t[n-i];
Line 815  unsigned int remqi(Q a,unsigned int mod)
Line 891  unsigned int remqi(Q a,unsigned int mod)
   return c;    return c;
 }  }
   
 extern int DP_Print;  
   
 #define F4_INTRAT_PERIOD 8  
   
 int generic_gauss_elim(MAT mat,MAT *nm,Z *dn,int **rindp,int **cindp)  int generic_gauss_elim(MAT mat,MAT *nm,Z *dn,int **rindp,int **cindp)
 {  {
   int **wmat;    int **wmat;
Line 828  int generic_gauss_elim(MAT mat,MAT *nm,Z *dn,int **rin
Line 900  int generic_gauss_elim(MAT mat,MAT *nm,Z *dn,int **rin
   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;
     MAT mat2,nm2;
     Z dn2;
     int *rind2,*cind2;
     int ret2;
   
   #if SIZEOF_LONG == 8
     ret = generic_gauss_elim64(mat,nm,dn,rindp,cindp);
     return ret;
   #endif
   bmat = (Z **)mat->body;    bmat = (Z **)mat->body;
   row = mat->row; col = mat->col;    row = mat->row; col = mat->col;
   wmat = (int **)almat(row,col);    wmat = (int **)almat(row,col);
Line 1105  int generic_gauss_elim_direct(MAT mat,MAT *nm,Z *dn,in
Line 1185  int generic_gauss_elim_direct(MAT mat,MAT *nm,Z *dn,in
   return rank;    return rank;
 }  }
   
   int mpz_intmtoratm(mpz_t **mat,int row,int col,mpz_t md,mpz_t **nm,mpz_t dn)
   {
     mpz_t t,s,b,u,nm1,dn1;
     int i,j,k,l,ret;
     mpz_t *mi,*nmk;
   
     if ( UNIMPZ(md) )
       return 0;
     mpz_init(t); mpz_init(s); mpz_init(b); mpz_init(u);
     mpz_init(nm1); mpz_init(dn1);
     mpz_fdiv_q_2exp(t,md,1); mpz_sqrt(s,t); mpz_fdiv_q_2exp(b,s,64);
     if ( !mpz_sgn(b) ) mpz_set_ui(b,1);
     mpz_set_ui(dn,1);
     for ( i = 0; i < row; i++ )
       for ( j = 0, mi = mat[i]; j < col; j++ )
         if ( mpz_sgn(mi[j]) ) {
           mpz_mul(s,mi[j],dn);
           mpz_mod(u,s,md);
           ret = mpz_inttorat(u,md,b,nm1,dn1);
           if ( !ret )
             return 0;
           else {
             if ( !UNIMPZ(dn1) ) {
               for ( k = 0; k < i; k++ )
                 for ( l = 0, nmk = nm[k]; l < col; l++ ) mpz_mul(nmk[l],nmk[l],dn1);
               for ( l = 0, nmk = nm[i]; l < j; l++ ) mpz_mul(nmk[l],nmk[l],dn1);
             }
             mpz_set(nm[i][j],nm1);
             mpz_mul(dn,dn,dn1);
           }
         }
     return 1;
   }
   
 int intmtoratm(MAT mat,Z md,MAT nm,Z *dn)  int intmtoratm(MAT mat,Z md,MAT nm,Z *dn)
 {  {
   Z t,s,b,dn0,dn1,nm1,q,u,unm,udn,dmy;    Z t,s,b,dn0,dn1,nm1,q,u,unm,udn,dmy;
Line 1205  int intvtoratv(Z *v,int n,Z md,Z b,Z *nm,Z *dn)
Line 1319  int intvtoratv(Z *v,int n,Z md,Z b,Z *nm,Z *dn)
   
 /* assuming 0 < c < m */  /* assuming 0 < c < m */
   
   int mpz_inttorat(mpz_t c,mpz_t m,mpz_t b,mpz_t nm,mpz_t dn)
   {
     mpz_t u1,v1,u2,v2,r1,r2;
     mpz_t q,t;
   
     mpz_init_set_ui(u1,0); mpz_init_set_ui(v1,1);
     mpz_init_set(u2,m); mpz_init_set(v2,c);
     mpz_init(q); mpz_init(t); mpz_init(r1); mpz_init(r2);
     while ( mpz_cmp(v2,b) >= 0 ) {
       /* r2 = u2-q*v2 */
       mpz_fdiv_qr(q,r2,u2,v2);
       mpz_set(u2,v2); mpz_set(v2,r2);
       /* r1 = u1-q*v1 */
       mpz_mul(t,q,v1); mpz_sub(r1,u1,t);
       mpz_set(u1,v1); mpz_set(v1,r1);
     }
     if ( mpz_cmp(v1,b) >= 0 ) return 0;
     else {
       if ( mpz_sgn(v1)<0  ) {
         mpz_neg(dn,v1); mpz_neg(nm,v2);
       } else {
         mpz_set(dn,v1); mpz_set(nm,v2);
       }
       return 1;
     }
   }
   
 int inttorat(Z c,Z m,Z b,Z *nmp,Z *dnp)  int inttorat(Z c,Z m,Z b,Z *nmp,Z *dnp)
 {  {
   Z qq,t,u1,v1,r1;    Z qq,t,u1,v1,r1;
Line 1228  int inttorat(Z c,Z m,Z b,Z *nmp,Z *dnp)
Line 1369  int inttorat(Z c,Z m,Z b,Z *nmp,Z *dnp)
   
 extern int f4_nocheck;  extern int f4_nocheck;
   
   int mpz_gensolve_check(MAT mat,mpz_t **nm,mpz_t dn,int rank,int *rind,int *cind)
   {
     int row,col,clen,i,j,k,l;
     mpz_t t;
     mpz_t *w;
     Z *mati;
     mpz_t *nmk;
   
     if ( f4_nocheck ) return 1;
     row = mat->row; col = mat->col; clen = col-rank;
     w = (mpz_t *)MALLOC(clen*sizeof(mpz_t));
     mpz_init(t);
     for ( i = 0; i < clen; i++ ) mpz_init(w[i]);
     for ( i = 0; i < row; i++ ) {
       mati = (Z *)mat->body[i];
       for ( l = 0; l < clen; l++ ) mpz_set_ui(w[l],0);
       for ( k = 0; k < rank; k++ )
         for ( l = 0, nmk = (mpz_t *)nm[k]; l < clen; l++ ) {
           /* w[l] += mati[rind[k]]*nmk[k] */
           if ( mati[rind[k]] ) mpz_addmul(w[l],BDY(mati[rind[k]]),nmk[l]);
         }
       for ( j = 0; j < clen; j++ ) {
         if ( mati[cind[j]] ) mpz_mul(t,dn,BDY(mati[cind[j]]));
         else mpz_set_ui(t,0);
         if ( mpz_cmp(w[j],t) ) break;
       }
       if ( j != clen ) break;
     }
     if ( i != row ) return 0;
     else return 1;
   }
   
 int gensolve_check(MAT mat,MAT nm,Z dn,int *rind,int *cind)  int gensolve_check(MAT mat,MAT nm,Z dn,int *rind,int *cind)
 {  {
   int row,col,rank,clen,i,j,k,l;    int row,col,rank,clen,i,j,k,l;
Line 1300  void isqrtz(Z a,Z *r)
Line 1473  void isqrtz(Z a,Z *r)
   else {    else {
     k = z_bits((Q)a); /* a <= 2^k-1 */      k = z_bits((Q)a); /* a <= 2^k-1 */
     bshiftz(ONE,-((k>>1)+(k&1)),&x); /* a <= x^2 */      bshiftz(ONE,-((k>>1)+(k&1)),&x); /* a <= x^2 */
     STOQ(2,two);      STOZ(2,two);
     while ( 1 ) {      while ( 1 ) {
       pwrz(x,two,&t);        pwrz(x,two,&t);
       if ( cmpz(t,a) <= 0 ) {        if ( cmpz(t,a) <= 0 ) {
Line 1411  int generic_gauss_elim_hensel(MAT mat,MAT *nmmat,Z *dn
Line 1584  int generic_gauss_elim_hensel(MAT mat,MAT *nmmat,Z *dn
   int *rind,*cind;    int *rind,*cind;
   int count;    int count;
   int ret;    int ret;
   struct oEGT eg_mul,eg_inv,eg_intrat,eg_check,tmp0,tmp1;    struct oEGT eg_mul1,eg_mul2,tmp0,tmp1,tmp2;
   int period;    int period;
   int *wx,*ptr;    int *wx,*ptr;
   int wxsize,nsize;    int wxsize,nsize;
   Z wn;    Z wn;
   Z wq;    Z wq;
   
   #if SIZEOF_LONG == 8
     return generic_gauss_elim_hensel64(mat,nmmat,dn,rindp,cindp);
   #endif
   init_eg(&eg_mul1); init_eg(&eg_mul2);
   a0 = (Z **)mat->body;    a0 = (Z **)mat->body;
   row = mat->row; col = mat->col;    row = mat->row; col = mat->col;
   w = (int **)almat(row,col);    w = (int **)almat(row,col);
   for ( ind = 0; ; ind++ ) {    for ( ind = 0; ; ind++ ) {
     md = get_lprime(ind);      md = get_lprime(ind);
     STOQ(md,mdq);      STOZ(md,mdq);
     for ( i = 0; i < row; i++ )      for ( i = 0; i < row; i++ )
       for ( j = 0, ai = a0[i], wi = w[i]; j < col; j++ )        for ( j = 0, ai = a0[i], wi = w[i]; j < col; j++ )
         wi[j] = remqi((Q)ai[j],md);          wi[j] = remqi((Q)ai[j],md);
Line 1468  int generic_gauss_elim_hensel(MAT mat,MAT *nmmat,Z *dn
Line 1645  int generic_gauss_elim_hensel(MAT mat,MAT *nmmat,Z *dn
       MKMAT(xmat,rank,ri); x = (Z **)(xmat)->body;        MKMAT(xmat,rank,ri); x = (Z **)(xmat)->body;
       MKMAT(*nmmat,rank,ri); nm = (Z **)(*nmmat)->body;        MKMAT(*nmmat,rank,ri); nm = (Z **)(*nmmat)->body;
       wc = (int **)almat(rank,ri);        wc = (int **)almat(rank,ri);
       for ( i = 0; i < rank; i++ )  
         wc[i] = w[i]+rank;  
       *rindp = rind = (int *)MALLOC_ATOMIC(rank*sizeof(int));        *rindp = rind = (int *)MALLOC_ATOMIC(rank*sizeof(int));
       *cindp = cind = (int *)MALLOC_ATOMIC((ri)*sizeof(int));        *cindp = cind = (int *)MALLOC_ATOMIC((ri)*sizeof(int));
   
       period = F4_INTRAT_PERIOD;        period = F4_INTRAT_PERIOD;
       for ( q = ONE, count = 0; ; ) {        for ( q = ONE, count = 0; ; ) {
           /* check Ax=B mod q */
         if ( DP_Print > 3 )          if ( DP_Print > 3 )
           fprintf(stderr,"o");            fprintf(stderr,"o");
         /* wc = b mod md */          /* wc = b mod md */
         for ( i = 0; i < rank; i++ )          for ( i = 0; i < rank; i++ )
           for ( j = 0, bi = b[i], wi = wc[i]; j < ri; j++ ) {            for ( j = 0, bi = b[i], wi = wc[i]; j < ri; j++ )
             wi[j] = remqi((Q)bi[j],md);              wi[j] = remqi((Q)bi[j],md);
             if ( wi[j] && sgnz(bi[j]) < 0 )          /* wc = A^(-1)wc; wc is not normalized */
               wi[j] = md-wi[j];          solve_by_lu_mod(w,rank,md,wc,ri,0);
           }  
         /* wc = A^(-1)wc; wc is normalized */  
         solve_by_lu_mod(w,rank,md,wc,ri,1);  
         /* x += q*wc */          /* x += q*wc */
   get_eg(&tmp0);
         for ( i = 0; i < rank; i++ )          for ( i = 0; i < rank; i++ )
           for ( j = 0, wi = wc[i]; j < ri; j++ ) mul1addtoz(q,wi[j],&x[i][j]);            for ( j = 0, wi = wc[i]; j < ri; j++ ) mul1addtoz(q,wi[j],&x[i][j]);
         /* b =(A*wc+b)/md */          /* b =(b-A*wc)/md */
   get_eg(&tmp1); add_eg(&eg_mul1,&tmp0,&tmp1);
         for ( i = 0; i < rank; i++ )          for ( i = 0; i < rank; i++ )
           for ( j = 0; j < ri; j++ ) {            for ( j = 0; j < ri; j++ ) {
             u = b[i][j];              mpz_t uz;
             for ( k = 0; k < rank; k++ ) mul1addtoz(a[i][k],wc[k][j],&u);  
               if ( b[i][j] )
                 mpz_init_set(uz,BDY(b[i][j]));
               else
                 mpz_init_set_ui(uz,0);
               for ( k = 0; k < rank; k++ ) {
                 if ( a[i][k] && wc[k][j] ) {
                   if ( wc[k][j] < 0 )
                     mpz_addmul_ui(uz,BDY(a[i][k]),-wc[k][j]);
                   else
                     mpz_submul_ui(uz,BDY(a[i][k]),wc[k][j]);
                 }
               }
               MPZTOZ(uz,u);
             divsz(u,mdq,&b[i][j]);              divsz(u,mdq,&b[i][j]);
           }            }
   get_eg(&tmp2); add_eg(&eg_mul2,&tmp1,&tmp2);
         count++;          count++;
         /* q = q*md */          /* q = q*md */
         mulz(q,mdq,&u); q = u;          mulz(q,mdq,&u); q = u;
         if ( count == period ) {          if ( count == period ) {
           ret = intmtoratm(xmat,q,*nmmat,dn);            ret = intmtoratm(xmat,q,*nmmat,dn);
           if ( ret ) {            if ( ret ) {
               print_eg("MUL1",&eg_mul1);
               print_eg("MUL2",&eg_mul2);
             for ( j = k = l = 0; j < col; j++ )              for ( j = k = l = 0; j < col; j++ )
               if ( cinfo[j] )                if ( cinfo[j] )
                 rind[k++] = j;                  rind[k++] = j;
Line 1557  int generic_gauss_elim_hensel_dalg(MAT mat,DP *mb,MAT 
Line 1748  int generic_gauss_elim_hensel_dalg(MAT mat,DP *mb,MAT 
   w = (int **)almat(row,col);    w = (int **)almat(row,col);
   for ( ind = 0; ; ind++ ) {    for ( ind = 0; ; ind++ ) {
     md = get_lprime(ind);      md = get_lprime(ind);
     STOQ(md,mdq);      STOZ(md,mdq);
     for ( i = 0; i < row; i++ )      for ( i = 0; i < row; i++ )
       for ( j = 0, ai = a0[i], wi = w[i]; j < col; j++ )        for ( j = 0, ai = a0[i], wi = w[i]; j < col; j++ )
         wi[j] = remqi((Q)ai[j],md);          wi[j] = remqi((Q)ai[j],md);
Line 1583  int generic_gauss_elim_hensel_dalg(MAT mat,DP *mb,MAT 
Line 1774  int generic_gauss_elim_hensel_dalg(MAT mat,DP *mb,MAT 
     a = (Z **)almat_pointer(rank,rank); /* lhs mat */      a = (Z **)almat_pointer(rank,rank); /* lhs mat */
     MKMAT(bmat,rank,col-rank); b = (Z **)bmat->body; /* lhs mat */      MKMAT(bmat,rank,col-rank); b = (Z **)bmat->body; /* lhs mat */
     for ( j = li = ri = 0; j < col; j++ )      for ( j = li = ri = 0; j < col; j++ )
       if ( cinfo[j] ) {        if ( cinfo[j] > 0 ) {
         /* the column is in lhs */          /* the column is in lhs */
         for ( i = 0; i < rank; i++ ) {          for ( i = 0; i < rank; i++ ) {
           w[i][li] = w[i][j];            w[i][li] = w[i][j];
           a[i][li] = a0[rinfo[i]][j];            a[i][li] = a0[rinfo[i]][j];
         }          }
         li++;          li++;
       } else {        } else if ( !cinfo[j] ) {
         /* the column is in rhs */          /* the column is in rhs */
         for ( i = 0; i < rank; i++ )          for ( i = 0; i < rank; i++ )
           b[i][ri] = a0[rinfo[i]][j];            b[i][ri] = a0[rinfo[i]][j];
Line 1613  int generic_gauss_elim_hensel_dalg(MAT mat,DP *mb,MAT 
Line 1804  int generic_gauss_elim_hensel_dalg(MAT mat,DP *mb,MAT 
       MKMAT(xmat,rank,ri); x = (Z **)(xmat)->body;        MKMAT(xmat,rank,ri); x = (Z **)(xmat)->body;
       MKMAT(*nmmat,rank,ri); nm = (Z **)(*nmmat)->body;        MKMAT(*nmmat,rank,ri); nm = (Z **)(*nmmat)->body;
       wc = (int **)almat(rank,ri);        wc = (int **)almat(rank,ri);
       for ( i = 0; i < rank; i++ )  
         wc[i] = w[i]+rank;  
       *rindp = rind = (int *)MALLOC_ATOMIC(rank*sizeof(int));        *rindp = rind = (int *)MALLOC_ATOMIC(rank*sizeof(int));
       *cindp = cind = (int *)MALLOC_ATOMIC((ri)*sizeof(int));        *cindp = cind = (int *)MALLOC_ATOMIC((ri)*sizeof(int));
   
Line 1624  int generic_gauss_elim_hensel_dalg(MAT mat,DP *mb,MAT 
Line 1813  int generic_gauss_elim_hensel_dalg(MAT mat,DP *mb,MAT 
           fprintf(stderr,"o");            fprintf(stderr,"o");
         /* wc = b mod md */          /* wc = b mod md */
         for ( i = 0; i < rank; i++ )          for ( i = 0; i < rank; i++ )
           for ( j = 0, bi = b[i], wi = wc[i]; j < ri; j++ ) {            for ( j = 0, bi = b[i], wi = wc[i]; j < ri; j++ )
             wi[j] = remqi((Q)bi[j],md);              wi[j] = remqi((Q)bi[j],md);
             if ( wi[j] && sgnz(bi[j]) < 0 )  
               wi[j] = md-wi[j];  
           }  
         /* wc = A^(-1)wc; wc is normalized */          /* wc = A^(-1)wc; wc is normalized */
         solve_by_lu_mod(w,rank,md,wc,ri,1);          solve_by_lu_mod(w,rank,md,wc,ri,1);
         /* x += q*wc */          /* x += q*wc */
         for ( i = 0; i < rank; i++ )          for ( i = 0; i < rank; i++ )
           for ( j = 0, wi = wc[i]; j < ri; j++ ) mul1addtoz(q,wi[j],&x[i][j]);            for ( j = 0, wi = wc[i]; j < ri; j++ ) mul1addtoz(q,wi[j],&x[i][j]);
         /* b =(A*wc+b)/md */          /* b =(b-A*wc)/md */
         for ( i = 0; i < rank; i++ )          for ( i = 0; i < rank; i++ )
           for ( j = 0; j < ri; j++ ) {            for ( j = 0; j < ri; j++ ) {
             u = b[i][j];              mpz_t uz;
             for ( k = 0; k < rank; k++ ) mul1addtoz(a[i][k],wc[k][j],&u);  
               if ( b[i][j] )
                 mpz_init_set(uz,BDY(b[i][j]));
               else
                 mpz_init_set_ui(uz,0);
               for ( k = 0; k < rank; k++ ) {
                 if ( a[i][k] && wc[k][j] ) {
                   if ( wc[k][j] < 0 )
                     mpz_addmul_ui(uz,BDY(a[i][k]),-wc[k][j]);
                   else
                     mpz_submul_ui(uz,BDY(a[i][k]),wc[k][j]);
                 }
               }
               MPZTOZ(uz,u);
             divsz(u,mdq,&b[i][j]);              divsz(u,mdq,&b[i][j]);
           }            }
         count++;          count++;
Line 1669  int generic_gauss_elim_hensel_dalg(MAT mat,DP *mb,MAT 
Line 1868  int generic_gauss_elim_hensel_dalg(MAT mat,DP *mb,MAT 
       }        }
   }    }
 }  }
   
   #if SIZEOF_LONG == 8
   mp_limb_t remqi64(Q a,mp_limb_t mod)
   {
     mp_limb_t c,nm,dn;
     mpz_t r;
   
     if ( !a ) return 0;
     else if ( a->z ) {
       mpz_init(r);
       c = mpz_fdiv_r_ui(r,BDY((Z)a),mod);
     } else {
       mpz_init(r);
       nm = mpz_fdiv_r_ui(r,mpq_numref(BDY(a)),mod);
       dn = mpz_fdiv_r_ui(r,mpq_denref(BDY(a)),mod);
       dn = invmod64(dn,mod);
       c = mulmod64(nm,dn,mod);
     }
     return c;
   }
   
   int generic_gauss_elim_mod64(mp_limb_t **mat,int row,int col,mp_limb_t md,int *colstat);
   mp_limb_t get_lprime64(int ind);
   
   void mpz_print(mpz_t a)
   {
     mpz_out_str(stdout,10,a); printf("\n");
   }
   
   void mpz_printmat(mpz_t **a,int row,int col)
   {
     int i,j;
     for ( i = 0; i < row; i++ ) {
       for ( j = 0; j < col; j++ ) {
         mpz_out_str(stdout,10,a[i][j]); printf(" ");
       }
       printf("\n");
     }
   }
   
   mpz_t **mpz_allocmat(int row,int col)
   {
     mpz_t **p;
     int i,j;
   
     p = (mpz_t **)MALLOC(row*sizeof(mpz_t *));
     for ( i = 0; i < row; i++ ) {
       p[i] = (mpz_t *)MALLOC(col*sizeof(mpz_t));
       for ( j = 0; j < col; j++ ) mpz_init(p[i][j]);
     }
     return p;
   }
   
   #if 1
   int generic_gauss_elim64(MAT mat,MAT *nm,Z *dn,int **rindp,int **cindp)
   {
     mp_limb_t **wmat;
     mp_limb_t *wmi;
     mp_limb_t md,inv,t,t1;
     Z z;
     Z **bmat,*bmi;
     mpz_t **tmat,**num;
     mpz_t *tmi;
     mpz_t den;
     mpz_t q,m1,m3,s,u;
     int *colstat,*wcolstat,*rind,*cind;
     int row,col,ind,i,j,k,l,rank,rank0;
     MAT r;
     int ret;
   
     bmat = (Z **)mat->body;
     row = mat->row; col = mat->col;
     wmat = (mp_limb_t **)almat64(row,col);
     colstat = (int *)MALLOC_ATOMIC(col*sizeof(int));
     wcolstat = (int *)MALLOC_ATOMIC(col*sizeof(int));
     mpz_init(m1); mpz_init(m3); mpz_init(den);
     for ( ind = 0; ; ind++ ) {
       if ( DP_Print ) {
         fprintf(asir_out,"."); fflush(asir_out);
       }
       md = get_lprime64(ind);
       for ( i = 0; i < row; i++ )
         for ( j = 0, bmi = bmat[i], wmi = wmat[i]; j < col; j++ )
           wmi[j] = bmi[j]==0?0:mpz_fdiv_ui(BDY(bmi[j]),md);
       rank = generic_gauss_elim_mod64(wmat,row,col,md,wcolstat);
       if ( !ind ) {
   RESET:
         mpz_set_ui(m1,md);
         rank0 = rank;
         bcopy(wcolstat,colstat,col*sizeof(int));
         // crmat
         tmat = mpz_allocmat(rank,col-rank);
         //
         num = mpz_allocmat(rank,col-rank);
         for ( i = 0; i < rank; i++ )
           for ( j = k = 0, tmi = tmat[i], wmi = wmat[i]; j < col; j++ )
             if ( !colstat[j] ) { mpz_set_ui(tmi[k],wmi[j]); k++; }
       } else {
         if ( rank < rank0 ) {
           if ( DP_Print ) {
             fprintf(asir_out,"lower rank matrix; continuing...\n");
             fflush(asir_out);
           }
           continue;
         } else if ( rank > rank0 ) {
           if ( DP_Print ) {
             fprintf(asir_out,"higher rank matrix; resetting...\n");
             fflush(asir_out);
           }
           goto RESET;
         } else {
           for ( j = 0; (j<col) && (colstat[j]==wcolstat[j]); j++ );
           if ( j < col ) {
             if ( DP_Print ) {
               fprintf(asir_out,"inconsitent colstat; resetting...\n");
               fflush(asir_out);
             }
             goto RESET;
           }
         }
   
         inv = invmod64(mpz_fdiv_ui(m1,md),md);
         mpz_mul_ui(m3,m1,md);
         for ( i = 0; i < rank; i++ )
           for ( j = k = 0, tmi = tmat[i], wmi = wmat[i]; j < col; j++ )
             if ( !colstat[j] ) {
               if ( mpz_sgn(tmi[k]) ) {
               /* f3 = f1+m1*(m1 mod md)^(-1)*(f2 - f1 mod md) */
                 t = mpz_fdiv_ui(tmi[k],md);
                 if ( wmi[j] >= t ) t = wmi[j]-t;
                 else t = md-(t-wmi[j]);
                 mpz_addmul_ui(tmi[k],m1,mulmod64(t,inv,md));
               } else if ( wmi[j] ) {
               /* f3 = m1*(m1 mod m2)^(-1)*f2 */
                 mpz_mul_ui(tmi[k],m1,mulmod64(wmi[j],inv,md));
               }
               k++;
             }
         mpz_set(m1,m3);
         if ( ind % F4_INTRAT_PERIOD )
           ret = 0;
         else
           ret = mpz_intmtoratm(tmat,rank,col-rank,m1,num,den);
         if ( ret ) {
           *rindp = rind = (int *)MALLOC_ATOMIC(rank*sizeof(int));
           *cindp = cind = (int *)MALLOC_ATOMIC((col-rank)*sizeof(int));
           for ( j = k = l = 0; j < col; j++ )
             if ( colstat[j] ) rind[k++] = j;
             else cind[l++] = j;
           if ( mpz_gensolve_check(mat,num,den,rank,rind,cind) ) {
             MKMAT(r,rank,col-rank); *nm = r;
             for ( i = 0; i < rank; i++ )
               for ( j = 0; j < col-rank; j++ ) {
                 MPZTOZ(num[i][j],z); BDY(r)[i][j] = z;
               }
             MPZTOZ(den,*dn);
             return rank;
           }
         }
       }
     }
   }
   #else
   int generic_gauss_elim64(MAT mat,MAT *nm,Z *dn,int **rindp,int **cindp)
   {
     mp_limb_t **wmat;
     mp_limb_t *wmi;
     mp_limb_t md,inv,t,t1;
     Z **bmat,**tmat,*bmi,*tmi;
     Z q,m1,m2,m3,s,u;
     int *colstat,*wcolstat,*rind,*cind;
     int row,col,ind,i,j,k,l,rank,rank0;
     MAT r,crmat;
     int ret;
   
     bmat = (Z **)mat->body;
     row = mat->row; col = mat->col;
     wmat = (mp_limb_t **)almat64(row,col);
     colstat = (int *)MALLOC_ATOMIC(col*sizeof(int));
     wcolstat = (int *)MALLOC_ATOMIC(col*sizeof(int));
     for ( ind = 0; ; ind++ ) {
       if ( DP_Print ) {
         fprintf(asir_out,"."); fflush(asir_out);
       }
       md = get_lprime64(ind);
       for ( i = 0; i < row; i++ )
         for ( j = 0, bmi = bmat[i], wmi = wmat[i]; j < col; j++ )
           wmi[j] = remqi64((Q)bmi[j],md);
       rank = generic_gauss_elim_mod64(wmat,row,col,md,wcolstat);
       if ( !ind ) {
   RESET:
         UTOZ(md,m1);
         rank0 = rank;
         bcopy(wcolstat,colstat,col*sizeof(int));
         MKMAT(crmat,rank,col-rank);
         MKMAT(r,rank,col-rank); *nm = r;
         tmat = (Z **)crmat->body;
         for ( i = 0; i < rank; i++ )
           for ( j = k = 0, tmi = tmat[i], wmi = wmat[i]; j < col; j++ )
             if ( !colstat[j] ) { UTOZ(wmi[j],tmi[k]); k++; }
       } else {
         if ( rank < rank0 ) {
           if ( DP_Print ) {
             fprintf(asir_out,"lower rank matrix; continuing...\n");
             fflush(asir_out);
           }
           continue;
         } else if ( rank > rank0 ) {
           if ( DP_Print ) {
             fprintf(asir_out,"higher rank matrix; resetting...\n");
             fflush(asir_out);
           }
           goto RESET;
         } else {
           for ( j = 0; (j<col) && (colstat[j]==wcolstat[j]); j++ );
           if ( j < col ) {
             if ( DP_Print ) {
               fprintf(asir_out,"inconsitent colstat; resetting...\n");
               fflush(asir_out);
             }
             goto RESET;
           }
         }
   
         inv = invmod64(remqi64((Q)m1,md),md);
         UTOZ(md,m2); mulz(m1,m2,&m3);
         for ( i = 0; i < rank; i++ )
           for ( j = k = 0, tmi = tmat[i], wmi = wmat[i]; j < col; j++ )
             if ( !colstat[j] ) {
               if ( tmi[k] ) {
               /* f3 = f1+m1*(m1 mod m2)^(-1)*(f2 - f1 mod m2) */
                 t = remqi64((Q)tmi[k],md);
                 if ( wmi[j] >= t ) t = wmi[j]-t;
                 else t = md-(t-wmi[j]);
                 t1 = mulmod64(t,inv,md);
                 UTOZ(t1,u); mulz(m1,u,&s);
                 addz(tmi[k],s,&u); tmi[k] = u;
               } else if ( wmi[j] ) {
               /* f3 = m1*(m1 mod m2)^(-1)*f2 */
                 t = mulmod64(wmi[j],inv,md);
                 UTOZ(t,u); mulz(m1,u,&s); tmi[k] = s;
               }
               k++;
             }
         m1 = m3;
         if ( ind % F4_INTRAT_PERIOD )
           ret = 0;
         else
           ret = intmtoratm(crmat,m1,*nm,dn);
         if ( ret ) {
           *rindp = rind = (int *)MALLOC_ATOMIC(rank*sizeof(int));
           *cindp = cind = (int *)MALLOC_ATOMIC((col-rank)*sizeof(int));
           for ( j = k = l = 0; j < col; j++ )
             if ( colstat[j] ) rind[k++] = j;
             else cind[l++] = j;
           if ( gensolve_check(mat,*nm,*dn,rind,cind) )
             return rank;
         }
       }
     }
   }
   #endif
   
   int generic_gauss_elim_hensel64(MAT mat,MAT *nmmat,Z *dn,int **rindp,int **cindp)
   {
     MAT r;
     Z z;
     Z **a0;
     Z *ai;
     mpz_t **a,**b,**x,**nm;
     mpz_t *bi,*xi;
     mpz_t q,u,den;
     mp_limb_t **w;
     mp_limb_t *wi;
     mp_limb_t **wc;
     mp_limb_t md;
     int row,col;
     int ind,i,j,k,l,li,ri,rank;
     int *cinfo,*rinfo;
     int *rind,*cind;
     int count;
     int ret;
     int period;
   
     a0 = (Z **)mat->body;
     row = mat->row; col = mat->col;
     w = (mp_limb_t **)almat64(row,col);
     for ( ind = 0; ; ind++ ) {
       md = get_lprime64(ind);
       for ( i = 0; i < row; i++ )
         for ( j = 0, ai = a0[i], wi = w[i]; j < col; j++ )
           wi[j] = remqi64((Q)ai[j],md);
   
       if ( DP_Print > 3 ) {
         fprintf(asir_out,"LU decomposition.."); fflush(asir_out);
       }
       rank = find_lhs_and_lu_mod64(w,row,col,md,&rinfo,&cinfo);
       if ( DP_Print > 3 ) {
         fprintf(asir_out,"done.\n"); fflush(asir_out);
       }
       a = (mpz_t **)mpz_allocmat(rank,rank); /* lhs mat */
       b = (mpz_t **)mpz_allocmat(rank,col-rank);
       for ( j = li = ri = 0; j < col; j++ )
         if ( cinfo[j] ) {
           /* the column is in lhs */
           for ( i = 0; i < rank; i++ ) {
             w[i][li] = w[i][j];
             if ( a0[rinfo[i]][j] )
               mpz_set(a[i][li],BDY(a0[rinfo[i]][j]));
             else
               mpz_set_ui(a[i][li],0);
           }
           li++;
         } else {
           /* the column is in rhs */
           for ( i = 0; i < rank; i++ ) {
             if ( a0[rinfo[i]][j] )
               mpz_set(b[i][ri],BDY(a0[rinfo[i]][j]));
             else
               mpz_set_ui(b[i][ri],0);
           }
           ri++;
         }
   
         /* solve Ax=B; A: rank x rank, B: rank x ri */
         /* algorithm
            c <- B
            x <- 0
            q <- 1
            do
              t <- A^(-1)c mod p
              x <- x+qt
              c <- (c-At)/p
              q <- qp
            end do
            then Ax-B=0 mod q and b=(B-Ax)/q hold after "do".
         */
         x = (mpz_t **)mpz_allocmat(rank,ri);
         nm = (mpz_t **)mpz_allocmat(rank,ri);
         wc = (mp_limb_t **)almat64(rank,ri);
         *rindp = rind = (int *)MALLOC_ATOMIC(rank*sizeof(int));
         *cindp = cind = (int *)MALLOC_ATOMIC((ri)*sizeof(int));
   
         period = F4_INTRAT_PERIOD;
         mpz_init_set_ui(q,1);
         mpz_init(u);
         mpz_init(den);
         for ( count = 0; ; ) {
           /* check Ax=B mod q */
           if ( DP_Print > 3 )
             fprintf(stderr,"o");
           /* wc = b mod md */
           for ( i = 0; i < rank; i++ )
             for ( j = 0, bi = b[i], wi = wc[i]; j < ri; j++ )
               wi[j] = mpz_fdiv_ui(bi[j],md);
           /* wc = A^(-1)wc; wc is not normalized */
           solve_by_lu_mod64(w,rank,md,wc,ri,0);
           /* x += q*wc */
           for ( i = 0; i < rank; i++ )
             for ( j = 0, wi = wc[i]; j < ri; j++ )
               if ( wi[j] > 0 )
                 mpz_addmul_ui(x[i][j],q,wi[j]);
               else if ( wi[j] < 0 )
                 mpz_submul_ui(x[i][j],q,-wi[j]);
           /* b =(b-A*wc)/md */
           for ( i = 0; i < rank; i++ )
             for ( j = 0; j < ri; j++ ) {
               mpz_set(u,b[i][j]);
               for ( k = 0; k < rank; k++ ) {
                 if ( a[i][k] && wc[k][j] ) {
                   if ( wc[k][j] < 0 )
                     mpz_addmul_ui(u,a[i][k],-wc[k][j]);
                   else
                     mpz_submul_ui(u,a[i][k],wc[k][j]);
                 }
               }
               mpz_divexact_ui(b[i][j],u,md);
             }
           count++;
           /* q = q*md */
           mpz_mul_ui(q,q,md);
           fprintf(stderr,".");
           if ( count == period ) {
             ret = mpz_intmtoratm(x,rank,ri,q,nm,den);
             if ( ret ) {
               for ( j = k = l = 0; j < col; j++ )
                 if ( cinfo[j] )
                   rind[k++] = j;
                 else
                   cind[l++] = j;
               ret = mpz_gensolve_check(mat,nm,den,rank,rind,cind);
               if ( ret ) {
                 *rindp = rind;
                 *cindp = cind;
                 for ( j = k = 0; j < col; j++ )
                   if ( !cinfo[j] )
                     cind[k++] = j;
                 MKMAT(r,rank,ri); *nmmat = r;
                 for ( i = 0; i < rank; i++ )
                   for ( j = 0; j < ri; j++ ) {
                     MPZTOZ(nm[i][j],z); BDY(r)[i][j] = z;
                   }
                 MPZTOZ(den,*dn);
                 return rank;
               }
             } else {
               fprintf(stderr,"F");
               period = period*3/2;
               count = 0;
             }
           }
         }
     }
   }
   
   #endif
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.10
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>