[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.4 and 1.7

version 1.4, 2018/09/25 07:36:01 version 1.7, 2018/10/01 05:54:09
Line 1 
Line 1 
   /* $OpenXM: OpenXM_contrib2/asir2018/engine/Q.c,v 1.6 2018/10/01 05:49:06 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 8
   
   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 31  void init_gmpq()
Line 36  void init_gmpq()
   gmp_randinit_default(GMP_RAND);    gmp_randinit_default(GMP_RAND);
 }  }
   
   void printexpr(VL,Obj);
   
 void pmat(Z **a,int row,int col)  void pmat(Z **a,int row,int col)
 {  {
   int i,j;    int i,j;
   
   for ( i = 0; i < row; i++, printf("\n") )    for ( i = 0; i < row; i++, printf("\n") )
     for ( j = 0; j < col; j++, printf(" ") )      for ( j = 0; j < col; j++, printf(" ") )
       printexpr(CO,a[i][j]);        printexpr(CO,(Obj)a[i][j]);
   printf("\n");    printf("\n");
 }  }
   
Line 353  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 648  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 705  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 884  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 898  int generic_gauss_elim(MAT mat,MAT *nm,Z *dn,int **rin
Line 901  int generic_gauss_elim(MAT mat,MAT *nm,Z *dn,int **rin
   MAT r,crmat;    MAT r,crmat;
   int ret;    int ret;
   
   #if SIZEOF_LONG == 8
     return generic_gauss_elim64(mat,nm,dn,rindp,cindp);
   #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 1369  void isqrtz(Z a,Z *r)
Line 1375  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 1493  init_eg(&eg_mul1); init_eg(&eg_mul2);
Line 1499  init_eg(&eg_mul1); init_eg(&eg_mul2);
   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 1641  int generic_gauss_elim_hensel_dalg(MAT mat,DP *mb,MAT 
Line 1647  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 1761  int generic_gauss_elim_hensel_dalg(MAT mat,DP *mb,MAT 
Line 1767  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);
   
   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
Legend:
Removed from v.1.4  
changed lines
  Added in v.1.7
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>