=================================================================== RCS file: /home/cvs/OpenXM_contrib2/asir2018/engine/Q.c,v retrieving revision 1.1 retrieving revision 1.13 diff -u -p -r1.1 -r1.13 --- OpenXM_contrib2/asir2018/engine/Q.c 2018/09/19 05:45:07 1.1 +++ OpenXM_contrib2/asir2018/engine/Q.c 2019/03/18 07:00:33 1.13 @@ -1,3 +1,4 @@ +/* $OpenXM: OpenXM_contrib2/asir2018/engine/Q.c,v 1.12 2019/03/17 02:48:13 noro Exp $ */ #include "ca.h" #include "gmp.h" #include "base.h" @@ -10,6 +11,10 @@ Z current_mod_lf; int current_mod_lf_size; gmp_randstate_t GMP_RAND; +#define F4_INTRAT_PERIOD 4 + +extern int DP_Print; + void isqrtz(Z a,Z *r); void bshiftz(Z a,int n,Z *r); @@ -25,12 +30,24 @@ void gc_free(void *p,size_t size) 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); 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) { mpz_t z; @@ -161,6 +178,7 @@ void mulz(Z n1,Z n2,Z *nr) void muladdtoz(Z n1,Z n2,Z *nr) { +#if 0 Z t; if ( n1 && n2 ) { @@ -168,13 +186,21 @@ void muladdtoz(Z n1,Z n2,Z *nr) NEWZ(t); mpz_init(BDY(t)); *nr = t; } 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 */ void mul1addtoz(Z n1,long u,Z *nr) { +#if 0 Z t; if ( n1 && u ) { @@ -185,7 +211,14 @@ void mul1addtoz(Z n1,long u,Z *nr) mpz_addmul_ui(BDY(*nr),BDY(n1),(unsigned long)u); else 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) @@ -327,7 +360,7 @@ void pwrz(Z n1,Z n,Z *nr) } else if ( !smallz(n) ) { error("exponent too big."); *nr = 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 { MPZTOMPQ(BDY(n1),q); MPQTOQ(q,r); pwrq(r,(Q)n,&p); *nr = (Z)p; @@ -427,7 +460,7 @@ void gcdvz_estimate(VECT v,Z *q) else addz(s,b[i],&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); else addz(t,b[i],&u); t = u; @@ -435,11 +468,54 @@ void gcdvz_estimate(VECT v,Z *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) { mpz_t a; mpz_init(a); - mpz_fac_ui(a,n); + mpz_fac_ui(a,(unsigned long)n); MPZTOZ(a,*nr); } @@ -579,7 +655,7 @@ void pwrq(Q n1,Q n,Q *nr) } else if ( !smallz((Z)n) ) { error("exponent too big."); *nr = 0; } else { - e = QTOS(n); + e = ZTOS(n); if ( e < 0 ) { e = -e; if ( n1->z ) { @@ -636,8 +712,8 @@ void mkbc(int n,Z *t) Z c,d,iq; for ( t[0] = ONE, i = 1; i <= n/2; i++ ) { - STOQ(n-i+1,c); mulz(t[i-1],c,&d); - STOQ(i,iq); divsz(d,iq,&t[i]); + STOZ(n-i+1,c); mulz(t[i-1],c,&d); + STOZ(i,iq); divsz(d,iq,&t[i]); } for ( ; i <= n; i++ ) t[i] = t[n-i]; @@ -815,10 +891,6 @@ unsigned int remqi(Q a,unsigned int mod) 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 **wmat; @@ -828,7 +900,15 @@ 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; MAT r,crmat; 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; row = mat->row; col = mat->col; wmat = (int **)almat(row,col); @@ -1105,6 +1185,40 @@ int generic_gauss_elim_direct(MAT mat,MAT *nm,Z *dn,in 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) { Z t,s,b,dn0,dn1,nm1,q,u,unm,udn,dmy; @@ -1115,7 +1229,7 @@ int intmtoratm(MAT mat,Z md,MAT nm,Z *dn) return 0; row = mat->row; col = mat->col; bshiftz(md,1,&t); - isqrtz(t,&s); + isqrt(t,&s); bshiftz(s,64,&b); if ( !b ) b = ONE; dn0 = ONE; @@ -1205,6 +1319,33 @@ int intvtoratv(Z *v,int n,Z md,Z b,Z *nm,Z *dn) /* 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) { Z qq,t,u1,v1,r1; @@ -1228,6 +1369,38 @@ int inttorat(Z c,Z m,Z b,Z *nmp,Z *dnp) extern int f4_nocheck; +int mpz_gensolve_check(MAT mat,mpz_t **nm,mpz_t dn,int rank,int clen,int *rind,int *cind) +{ + int row,col,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; + 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 row,col,rank,clen,i,j,k,l; @@ -1296,11 +1469,11 @@ void isqrtz(Z a,Z *r) Z two; if ( !a ) *r = 0; - else if ( UNIQ(a) ) *r = ONE; + else if ( UNIZ(a) ) *r = ONE; else { k = z_bits((Q)a); /* a <= 2^k-1 */ bshiftz(ONE,-((k>>1)+(k&1)),&x); /* a <= x^2 */ - STOQ(2,two); + STOZ(2,two); while ( 1 ) { pwrz(x,two,&t); if ( cmpz(t,a) <= 0 ) { @@ -1411,19 +1584,23 @@ int generic_gauss_elim_hensel(MAT mat,MAT *nmmat,Z *dn int *rind,*cind; int count; 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 *wx,*ptr; int wxsize,nsize; Z wn; 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); a0 = (Z **)mat->body; row = mat->row; col = mat->col; w = (int **)almat(row,col); for ( ind = 0; ; ind++ ) { md = get_lprime(ind); - STOQ(md,mdq); + STOZ(md,mdq); for ( i = 0; i < row; i++ ) for ( j = 0, ai = a0[i], wi = w[i]; j < col; j++ ) wi[j] = remqi((Q)ai[j],md); @@ -1468,40 +1645,54 @@ int generic_gauss_elim_hensel(MAT mat,MAT *nmmat,Z *dn MKMAT(xmat,rank,ri); x = (Z **)(xmat)->body; MKMAT(*nmmat,rank,ri); nm = (Z **)(*nmmat)->body; wc = (int **)almat(rank,ri); - for ( i = 0; i < rank; i++ ) - wc[i] = w[i]+rank; *rindp = rind = (int *)MALLOC_ATOMIC(rank*sizeof(int)); *cindp = cind = (int *)MALLOC_ATOMIC((ri)*sizeof(int)); period = F4_INTRAT_PERIOD; for ( q = ONE, 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++ ) { + for ( j = 0, bi = b[i], wi = wc[i]; j < ri; j++ ) 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 */ - solve_by_lu_mod(w,rank,md,wc,ri,1); + /* wc = A^(-1)wc; wc is not normalized */ + solve_by_lu_mod(w,rank,md,wc,ri,0); /* x += q*wc */ +get_eg(&tmp0); for ( i = 0; i < rank; i++ ) 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 ( j = 0; j < ri; j++ ) { - u = b[i][j]; - for ( k = 0; k < rank; k++ ) mul1addtoz(a[i][k],wc[k][j],&u); + mpz_t uz; + + 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]); } +get_eg(&tmp2); add_eg(&eg_mul2,&tmp1,&tmp2); count++; /* q = q*md */ mulz(q,mdq,&u); q = u; if ( count == period ) { ret = intmtoratm(xmat,q,*nmmat,dn); if ( ret ) { + print_eg("MUL1",&eg_mul1); + print_eg("MUL2",&eg_mul2); for ( j = k = l = 0; j < col; j++ ) if ( cinfo[j] ) rind[k++] = j; @@ -1515,8 +1706,10 @@ int generic_gauss_elim_hensel(MAT mat,MAT *nmmat,Z *dn if ( !cinfo[j] ) cind[k++] = j; return rank; - } + } else + goto reset; } else { +reset: period = period*3/2; count = 0; } @@ -1552,12 +1745,15 @@ int generic_gauss_elim_hensel_dalg(MAT mat,DP *mb,MAT Z wq; DP m; +#if SIZEOF_LONG == 8 + return generic_gauss_elim_hensel64(mat,nmmat,dn,rindp,cindp,mb); +#endif a0 = (Z **)mat->body; row = mat->row; col = mat->col; w = (int **)almat(row,col); for ( ind = 0; ; ind++ ) { md = get_lprime(ind); - STOQ(md,mdq); + STOZ(md,mdq); for ( i = 0; i < row; i++ ) for ( j = 0, ai = a0[i], wi = w[i]; j < col; j++ ) wi[j] = remqi((Q)ai[j],md); @@ -1583,14 +1779,14 @@ int generic_gauss_elim_hensel_dalg(MAT mat,DP *mb,MAT a = (Z **)almat_pointer(rank,rank); /* lhs mat */ MKMAT(bmat,rank,col-rank); b = (Z **)bmat->body; /* lhs mat */ for ( j = li = ri = 0; j < col; j++ ) - if ( cinfo[j] ) { + if ( cinfo[j] > 0 ) { /* the column is in lhs */ for ( i = 0; i < rank; i++ ) { w[i][li] = w[i][j]; a[i][li] = a0[rinfo[i]][j]; } li++; - } else { + } else if ( !cinfo[j] ) { /* the column is in rhs */ for ( i = 0; i < rank; i++ ) b[i][ri] = a0[rinfo[i]][j]; @@ -1613,8 +1809,6 @@ int generic_gauss_elim_hensel_dalg(MAT mat,DP *mb,MAT MKMAT(xmat,rank,ri); x = (Z **)(xmat)->body; MKMAT(*nmmat,rank,ri); nm = (Z **)(*nmmat)->body; wc = (int **)almat(rank,ri); - for ( i = 0; i < rank; i++ ) - wc[i] = w[i]+rank; *rindp = rind = (int *)MALLOC_ATOMIC(rank*sizeof(int)); *cindp = cind = (int *)MALLOC_ATOMIC((ri)*sizeof(int)); @@ -1624,23 +1818,29 @@ int generic_gauss_elim_hensel_dalg(MAT mat,DP *mb,MAT 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++ ) { + for ( j = 0, bi = b[i], wi = wc[i]; j < ri; j++ ) 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 */ - solve_by_lu_mod(w,rank,md,wc,ri,1); + /* wc = A^(-1)wc; wc is not normalized */ + solve_by_lu_mod(w,rank,md,wc,ri,0); /* x += q*wc */ for ( i = 0; i < rank; i++ ) 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 ( j = 0; j < ri; j++ ) { - u = b[i][j]; - for ( k = 0; k < rank; k++ ) mul1addtoz(a[i][k],wc[k][j],&u); + mpz_t uz; + + 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] ) + mpz_submul_ui(uz,BDY(a[i][k]),wc[k][j]); + MPZTOZ(uz,u); divsz(u,mdq,&b[i][j]); } + count++; /* q = q*md */ mulz(q,mdq,&u); q = u; @@ -1660,8 +1860,439 @@ int generic_gauss_elim_hensel_dalg(MAT mat,DP *mb,MAT if ( !cinfo[j] ) cind[k++] = j; return rank; + } else + goto reset; + } else { +reset: + period = period*3/2; + count = 0; + } + } + } + } +} + +#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= 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,col-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= 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,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,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; } @@ -1669,3 +2300,5 @@ int generic_gauss_elim_hensel_dalg(MAT mat,DP *mb,MAT } } } + +#endif