[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.8 and 1.21

version 1.8, 2018/10/02 09:06:15 version 1.21, 2021/10/08 07:27:05
Line 1 
Line 1 
 /* $OpenXM: OpenXM_contrib2/asir2018/engine/Q.c,v 1.7 2018/10/01 05:54:09 noro Exp $ */  /* $OpenXM: OpenXM_contrib2/asir2018/engine/Q.c,v 1.20 2020/12/04 08:09:33 noro Exp $ */
 #include "ca.h"  #include "ca.h"
 #include "gmp.h"  #include "gmp.h"
 #include "base.h"  #include "base.h"
 #include "inline.h"  #include "inline.h"
   
 mpz_t ONEMPZ;  mpz_t ONEMPZ;
 Z ONE;  extern Z ONE;
 int lf_lazy;  int lf_lazy;
 Z current_mod_lf;  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  #define F4_INTRAT_PERIOD 4
   
 extern int DP_Print;  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);
   int mpz_inttorat(mpz_t c,mpz_t m,mpz_t b,mpz_t nm,mpz_t dn);
   int generic_gauss_elim_hensel64(MAT mat,MAT *nmmat,Z *dn,int **rindp,int **cindp,DP *mb);
   int find_lhs_and_lu_mod64(mp_limb_t **a,int row,int col,mp_limb_t md,int **rinfo,int **cinfo);
   void solve_by_lu_mod64(mp_limb_t **a,int n,mp_limb_t md,mp_limb_signed_t **b,int l,int normalize);
   
 void *gc_realloc(void *p,size_t osize,size_t nsize)  void *gc_realloc(void *p,size_t osize,size_t nsize)
 {  {
Line 30  void gc_free(void *p,size_t size)
Line 34  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);
Line 93  void dnq(Q q,Z *r)
Line 97  void dnq(Q q,Z *r)
 int sgnq(Q q)  int sgnq(Q q)
 {  {
   if ( !q ) return 0;    if ( !q ) return 0;
   else if ( q->z ) return mpz_sgn(BDY((Z)q));    else if ( q->z ) {
   else return mpz_sgn(mpq_numref(BDY(q)));     int sgn;
      sgn = mpz_sgn(BDY((Z)q));
      return sgn;
     } else return mpz_sgn(mpq_numref(BDY(q)));
 }  }
   
 Q mpqtozq(mpq_t a)  Q mpqtozq(mpq_t a)
Line 388  int cmpz(Z q1,Z q2)
Line 395  int cmpz(Z q1,Z q2)
     else if ( sgn < 0 ) return -1;      else if ( sgn < 0 ) return -1;
     else return 0;      else return 0;
   }    }
     /* XXX */
     return 0;
 }  }
   
 void gcdz(Z n1,Z n2,Z *nq)  void gcdz(Z n1,Z n2,Z *nq)
Line 515  void factorialz(unsigned int n,Z *nr)
Line 524  void factorialz(unsigned int n,Z *nr)
 {  {
   mpz_t a;    mpz_t a;
   mpz_init(a);    mpz_init(a);
   mpz_fac_ui(a,n);    mpz_fac_ui(a,(unsigned long)n);
   MPZTOZ(a,*nr);    MPZTOZ(a,*nr);
 }  }
   
Line 558  void subq(Q n1,Q n2,Q *nr)
Line 567  void subq(Q n1,Q n2,Q *nr)
   
   if ( !n1 ) {    if ( !n1 ) {
     if ( !n2 ) *nr = 0;      if ( !n2 ) *nr = 0;
     else if ( n1->z ) chsgnz((Z)n1,(Z *)nr);      else if ( n2->z ) chsgnz((Z)n2,(Z *)nr);
     else {      else {
         mpq_init(t); mpq_neg(t,BDY(n2)); MPQTOQ(t,*nr);          mpq_init(t); mpq_neg(t,BDY(n2)); MPQTOQ(t,*nr);
       }        }
Line 665  void pwrq(Q n1,Q n,Q *nr)
Line 674  void pwrq(Q n1,Q n,Q *nr)
         nm[0] = mpq_denref(BDY(n1))[0];          nm[0] = mpq_denref(BDY(n1))[0];
         dn[0] = mpq_numref(BDY(n1))[0];          dn[0] = mpq_numref(BDY(n1))[0];
       }        }
         if ( mpz_sgn(dn)<0 ) {
           mpz_neg(nm,nm);
           mpz_neg(dn,dn);
         }
     } else {      } else {
       if ( n1->z ) {        if ( n1->z ) {
         nm[0] = BDY((Z)n1)[0];          nm[0] = BDY((Z)n1)[0];
Line 762  void lgp(P p,Z *g,Z *l);
Line 775  void lgp(P p,Z *g,Z *l);
   
 void ptozp(P p,int sgn,Q *c,P *pr)  void ptozp(P p,int sgn,Q *c,P *pr)
 {  {
   Z nm,dn;    Z nm,dn,nm1;
   
   if ( !p ) {    if ( !p ) {
     *c = 0; *pr = 0;      *c = 0; *pr = 0;
   } else {    } else {
     lgp(p,&nm,&dn);      lgp(p,&nm,&dn);
       if ( sgn < 0 ) {
         chsgnz(nm,&nm1); nm = nm1;
       }
     divz(nm,dn,(Z *)c);      divz(nm,dn,(Z *)c);
     divsp(CO,p,(P)*c,pr);      divsp(CO,p,(P)*c,pr);
   }    }
Line 780  void lgp(P p,Z *g,Z *l)
Line 796  void lgp(P p,Z *g,Z *l)
   
   if ( NUM(p) ) {    if ( NUM(p) ) {
     if ( ((Q)p)->z ) {      if ( ((Q)p)->z ) {
       MPZTOZ(BDY((Z)p),*g);        absz((Z)p,g);
       *l = ONE;        *l = ONE;
     } else {      } else {
       MPZTOZ(mpq_numref(BDY((Q)p)),*g);        MPZTOZ(mpq_numref(BDY((Q)p)),g1); absz(g1,g);
       MPZTOZ(mpq_denref(BDY((Q)p)),*l);        MPZTOZ(mpq_denref(BDY((Q)p)),*l);
     }      }
   } else {    } else {
Line 1337  int mpz_inttorat(mpz_t c,mpz_t m,mpz_t b,mpz_t nm,mpz_
Line 1353  int mpz_inttorat(mpz_t c,mpz_t m,mpz_t b,mpz_t nm,mpz_
   }    }
   if ( mpz_cmp(v1,b) >= 0 ) return 0;    if ( mpz_cmp(v1,b) >= 0 ) return 0;
   else {    else {
       mpz_gcd(t,v1,v2);
       if ( UNIMPZ(t) )
         mpz_set_ui(r1,0);
       else {
         /* v1 /= t, v2 /= t, t=c*v1-v2, r1=t%m */
         mpz_divexact(v1,v1,t); mpz_divexact(v2,v2,t);
         mpz_mul(t,c,v1); mpz_sub(t,t,v2); mpz_mod(r1,t,m);
       }
       if ( mpz_sgn(r1) ) return 0;
     if ( mpz_sgn(v1)<0  ) {      if ( mpz_sgn(v1)<0  ) {
       mpz_neg(dn,v1); mpz_neg(nm,v2);        mpz_neg(dn,v1); mpz_neg(nm,v2);
     } else {      } else {
Line 1348  int mpz_inttorat(mpz_t c,mpz_t m,mpz_t b,mpz_t nm,mpz_
Line 1373  int mpz_inttorat(mpz_t c,mpz_t m,mpz_t b,mpz_t nm,mpz_
   
 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,s,r,u1,v1,r1;
   Z q,u2,v2,r2;    Z q,u2,v2,r2;
   
   u1 = 0; v1 = ONE; u2 = m; v2 = c;    u1 = 0; v1 = ONE; u2 = m; v2 = c;
Line 1358  int inttorat(Z c,Z m,Z b,Z *nmp,Z *dnp)
Line 1383  int inttorat(Z c,Z m,Z b,Z *nmp,Z *dnp)
   }    }
   if ( cmpz(v1,b) >= 0 ) return 0;    if ( cmpz(v1,b) >= 0 ) return 0;
   else {    else {
     if ( mpz_sgn(BDY(v1))<0  ) {      /* reduction and check */
       chsgnz(v1,dnp); chsgnz(v2,nmp);      /* v2/v1 = u2/u1, c*u1-u2 = 0 mod m? */
       gcdz(v1,v2,&t);
       if ( UNIZ(t) ) {
         u1 = v1; u2 = v2; r = 0;
     } else {      } else {
       *dnp = v1; *nmp = v2;        divsz(v1,t,&u1); divsz(v2,t,&u2);
         mulz(c,u1,&t); subz(t,u2,&s); remz(s,m,&r);
     }      }
       if ( r ) return 0;
       if ( mpz_sgn(BDY(u1))<0  ) {
         chsgnz(u1,dnp); chsgnz(u2,nmp);
       } else {
         *dnp = u1; *nmp = u2;
       }
     return 1;      return 1;
   }    }
 }  }
   
 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 mpz_gensolve_check(MAT mat,mpz_t **nm,mpz_t dn,int rank,int clen,int *rind,int *cind)
 {  {
   int row,col,clen,i,j,k,l;    int row,col,i,j,k,l;
   mpz_t t;    mpz_t t;
   mpz_t *w;    mpz_t *w;
   Z *mati;    Z *mati;
   mpz_t *nmk;    mpz_t *nmk;
   
   if ( f4_nocheck ) return 1;    if ( f4_nocheck ) return 1;
   row = mat->row; col = mat->col; clen = col-rank;    row = mat->row; col = mat->col;
   w = (mpz_t *)MALLOC(clen*sizeof(mpz_t));    w = (mpz_t *)MALLOC(clen*sizeof(mpz_t));
   mpz_init(t);    mpz_init(t);
   for ( i = 0; i < clen; i++ ) mpz_init(w[i]);    for ( i = 0; i < clen; i++ ) mpz_init(w[i]);
Line 1469  void isqrtz(Z a,Z *r)
Line 1504  void isqrtz(Z a,Z *r)
   Z two;    Z two;
   
   if ( !a ) *r = 0;    if ( !a ) *r = 0;
   else if ( UNIQ(a) ) *r = ONE;    else if ( UNIZ(a) ) *r = ONE;
   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 */
Line 1591  int generic_gauss_elim_hensel(MAT mat,MAT *nmmat,Z *dn
Line 1626  int generic_gauss_elim_hensel(MAT mat,MAT *nmmat,Z *dn
   Z wn;    Z wn;
   Z wq;    Z wq;
   
   #if SIZEOF_LONG == 8
     return generic_gauss_elim_hensel64(mat,nmmat,dn,rindp,cindp,0);
   #endif
 init_eg(&eg_mul1); init_eg(&eg_mul2);  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;
Line 1703  get_eg(&tmp2); add_eg(&eg_mul2,&tmp1,&tmp2);
Line 1741  get_eg(&tmp2); add_eg(&eg_mul2,&tmp1,&tmp2);
                 if ( !cinfo[j] )                  if ( !cinfo[j] )
                   cind[k++] = j;                    cind[k++] = j;
               return rank;                return rank;
             }              } else
                 goto reset;
           } else {            } else {
   reset:
             period = period*3/2;              period = period*3/2;
             count = 0;              count = 0;
           }            }
Line 1740  int generic_gauss_elim_hensel_dalg(MAT mat,DP *mb,MAT 
Line 1780  int generic_gauss_elim_hensel_dalg(MAT mat,DP *mb,MAT 
   Z wq;    Z wq;
   DP m;    DP m;
   
   #if SIZEOF_LONG == 8
     return generic_gauss_elim_hensel64(mat,nmmat,dn,rindp,cindp,mb);
   #endif
   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);
Line 1812  int generic_gauss_elim_hensel_dalg(MAT mat,DP *mb,MAT 
Line 1855  int generic_gauss_elim_hensel_dalg(MAT mat,DP *mb,MAT 
         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);
         /* wc = A^(-1)wc; wc is normalized */          /* wc = A^(-1)wc; wc is not normalized */
         solve_by_lu_mod(w,rank,md,wc,ri,1);          solve_by_lu_mod(w,rank,md,wc,ri,0);
         /* 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]);
Line 1826  int generic_gauss_elim_hensel_dalg(MAT mat,DP *mb,MAT 
Line 1869  int generic_gauss_elim_hensel_dalg(MAT mat,DP *mb,MAT 
               mpz_init_set(uz,BDY(b[i][j]));                mpz_init_set(uz,BDY(b[i][j]));
             else              else
               mpz_init_set_ui(uz,0);                mpz_init_set_ui(uz,0);
             for ( k = 0; k < rank; k++ ) {              for ( k = 0; k < rank; k++ )
               if ( a[i][k] && wc[k][j] ) {                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]);                    mpz_submul_ui(uz,BDY(a[i][k]),wc[k][j]);
               }  
             }  
             MPZTOZ(uz,u);              MPZTOZ(uz,u);
             divsz(u,mdq,&b[i][j]);              divsz(u,mdq,&b[i][j]);
           }            }
   
         count++;          count++;
         /* q = q*md */          /* q = q*md */
         mulz(q,mdq,&u); q = u;          mulz(q,mdq,&u); q = u;
Line 1856  int generic_gauss_elim_hensel_dalg(MAT mat,DP *mb,MAT 
Line 1895  int generic_gauss_elim_hensel_dalg(MAT mat,DP *mb,MAT 
                 if ( !cinfo[j] )                  if ( !cinfo[j] )
                   cind[k++] = j;                    cind[k++] = j;
               return rank;                return rank;
             }              } else
                 goto reset;
           } else {            } else {
   reset:
             period = period*3/2;              period = period*3/2;
             count = 0;              count = 0;
           }            }
Line 2014  RESET:
Line 2055  RESET:
         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 ( mpz_gensolve_check(mat,num,den,rank,rind,cind) ) {          if ( mpz_gensolve_check(mat,num,den,rank,col-rank,rind,cind) ) {
           MKMAT(r,rank,col-rank); *nm = r;            MKMAT(r,rank,col-rank); *nm = r;
           for ( i = 0; i < rank; i++ )            for ( i = 0; i < rank; i++ )
             for ( j = 0; j < col-rank; j++ ) {              for ( j = 0; j < col-rank; j++ ) {
Line 2127  RESET:
Line 2168  RESET:
   }    }
 }  }
 #endif  #endif
   
   int generic_gauss_elim_hensel64(MAT mat,MAT *nmmat,Z *dn,int **rindp,int **cindp,DP *mb)
   {
     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;
     DP m;
   
     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);
       }
   
       if ( mb ) {
         /* this part is added for inv_or_split_dalg */
         for ( i = 0; i < col-1; i++ ) {
           if ( !cinfo[i] ) {
             m = mb[i];
             for ( j = i+1; j < col-1; j++ )
               if ( dp_redble(mb[j],m) )
                 cinfo[j] = -1;
           }
         }
       }
   
       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] > 0 ) {
           /* 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 if ( !cinfo[j] ) {
           /* 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,(mp_limb_signed_t **)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);
           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] > 0 )
                   rind[k++] = j;
                 else if ( !cinfo[j] )
                   cind[l++] = j;
               ret = mpz_gensolve_check(mat,nm,den,rank,ri,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
                 goto reset;
             } else {
   reset:
               fprintf(stderr,"F");
               period = period*3/2;
               count = 0;
             }
           }
         }
     }
   }
   
 #endif  #endif
Legend:
Removed from v.1.8  
changed lines
  Added in v.1.21
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>