OpenXM_contrib2/asir2018/engine/Q.c - diff

Return to Q.c CVS log

Up to [local] / OpenXM_contrib2 / asir2018 / engine

Diff for /OpenXM_contrib2/asir2018/engine/Q.c between version 1.2 and 1.7

version 1.2, 2018/09/21 07:06:51

version 1.7, 2018/10/01 05:54:09

Line 1

/* $OpenXM: OpenXM_contrib2/asir2018/engine/Q.c,v 1.6 2018/10/01 05:49:06 noro Exp $ */

#include "ca.h"

#include "gmp.h"

#include "base.h"

Line 10 Z current_mod_lf;

Line 11 Z current_mod_lf;

int current_mod_lf_size;

gmp_randstate_t GMP_RAND;

#define F4_INTRAT_PERIOD 8

extern int DP_Print;

void isqrtz(Z a,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()

{

mp_set_memory_functions(Risa_GC_malloc,gc_realloc,gc_free);

mp_set_memory_functions(Risa_GC_malloc_atomic,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;

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)

{

#if 1

#if 0

Z t;

if ( n1 && n2 ) {

Line 171 void muladdtoz(Z n1,Z n2,Z *nr)

Line 188 void muladdtoz(Z n1,Z n2,Z *nr)

mpz_addmul(BDY(*nr),BDY(n1),BDY(n2));

if ( !mpz_sgn(BDY(*nr)) )

*nr = 0;

}

#else

Z t,s;

Line 183 void muladdtoz(Z n1,Z n2,Z *nr)

Line 200 void muladdtoz(Z n1,Z n2,Z *nr)

void mul1addtoz(Z n1,long u,Z *nr)

{

#if 0

Z t;

if ( n1 && u ) {

Line 196 void mul1addtoz(Z n1,long u,Z *nr)

Line 214 void mul1addtoz(Z n1,long u,Z *nr)

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)

Line 337 void pwrz(Z n1,Z n,Z *nr)

Line 360 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;

Line 437 void gcdvz_estimate(VECT v,Z *q)

Line 460 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;

Line 445 void gcdvz_estimate(VECT v,Z *q)

Line 468 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;

Line 589 void pwrq(Q n1,Q n,Q *nr)

Line 655 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 ) {

Line 646 void mkbc(int n,Z *t)

Line 712 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);

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++ )

t[i] = t[n-i];

Line 825 unsigned int remqi(Q a,unsigned int mod)

Line 891 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;

Line 839 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;

int ret;

#if SIZEOF_LONG == 8

return generic_gauss_elim64(mat,nm,dn,rindp,cindp);

#endif

bmat = (Z **)mat->body;

row = mat->row; col = mat->col;

wmat = (int **)almat(row,col);

Line 1310 void isqrtz(Z a,Z *r)

Line 1375 void isqrtz(Z a,Z *r)

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 ) {

Line 1421 int generic_gauss_elim_hensel(MAT mat,MAT *nmmat,Z *dn

Line 1486 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;

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);

Line 1478 int generic_gauss_elim_hensel(MAT mat,MAT *nmmat,Z *dn

Line 1544 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 )

/* 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 */

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];

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]);

}

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;

Line 1567 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);

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);

Line 1593 int generic_gauss_elim_hensel_dalg(MAT mat,DP *mb,MAT

Line 1673 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];

Line 1623 int generic_gauss_elim_hensel_dalg(MAT mat,DP *mb,MAT

Line 1703 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));

Line 1634 int generic_gauss_elim_hensel_dalg(MAT mat,DP *mb,MAT

Line 1712 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);

/* 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];

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]);

}

count++;

Line 1679 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

FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>