[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.9

version 1.8, 2018/10/02 09:06:15 version 1.9, 2018/10/19 23:27:38
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.8 2018/10/02 09:06:15 noro Exp $ */
 #include "ca.h"  #include "ca.h"
 #include "gmp.h"  #include "gmp.h"
 #include "base.h"  #include "base.h"
Line 11  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  #define F4_INTRAT_PERIOD 4
   
 extern int DP_Print;  extern int DP_Print;
   
Line 1591  int generic_gauss_elim_hensel(MAT mat,MAT *nmmat,Z *dn
Line 1591  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);
   #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 2127  RESET:
Line 2130  RESET:
   }    }
 }  }
 #endif  #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  #endif
Legend:
Removed from v.1.8  
changed lines
  Added in v.1.9
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>