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.8 and 1.17

version 1.8, 2018/10/02 09:06:15

version 1.17, 2020/01/09 01:47:41

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.16 2019/05/22 02:41:16 noro Exp $ */

#include "ca.h"

#include "gmp.h"

#include "base.h"

#include "inline.h"

mpz_t ONEMPZ;

Z ONE;

extern Z ONE;

int lf_lazy;

Z current_mod_lf;

int current_mod_lf_size;

gmp_randstate_t GMP_RAND;

#define F4_INTRAT_PERIOD 8

#define F4_INTRAT_PERIOD 4

extern int DP_Print;

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

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

}

Line 558 void subq(Q n1,Q n2,Q *nr)

if ( !n1 ) {

if ( !n2 ) *nr = 0;

else if ( n1->z ) chsgnz((Z)n1,(Z *)nr);

else if ( n2->z ) chsgnz((Z)n2,(Z *)nr);

else {

mpq_init(t); mpq_neg(t,BDY(n2)); MPQTOQ(t,*nr);

}

Line 762 void lgp(P p,Z *g,Z *l);

void ptozp(P p,int sgn,Q *c,P *pr)

{

Z nm,dn;

Z nm,dn,nm1;

if ( !p ) {

*c = 0; *pr = 0;

} else {

lgp(p,&nm,&dn);

if ( sgn < 0 ) {

chsgnz(nm,&nm1); nm = nm1;

}

divz(nm,dn,(Z *)c);

divsp(CO,p,(P)*c,pr);

}

Line 1229 int intmtoratm(MAT mat,Z md,MAT nm,Z *dn)

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

Line 1337 int mpz_inttorat(mpz_t c,mpz_t m,mpz_t b,mpz_t nm,mpz_

Line 1340 int mpz_inttorat(mpz_t c,mpz_t m,mpz_t b,mpz_t nm,mpz_

}

if ( mpz_cmp(v1,b) >= 0 ) return 0;

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

mpz_neg(dn,v1); mpz_neg(nm,v2);

} else {

Line 1348 int mpz_inttorat(mpz_t c,mpz_t m,mpz_t b,mpz_t nm,mpz_

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

{

Z qq,t,u1,v1,r1;

Z qq,t,s,r,u1,v1,r1;

Z q,u2,v2,r2;

u1 = 0; v1 = ONE; u2 = m; v2 = c;

Line 1358 int inttorat(Z c,Z m,Z b,Z *nmp,Z *dnp)

Line 1370 int inttorat(Z c,Z m,Z b,Z *nmp,Z *dnp)

}

if ( cmpz(v1,b) >= 0 ) return 0;

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 {

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

}

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 *w;

Z *mati;

mpz_t *nmk;

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

mpz_init(t);

for ( i = 0; i < clen; i++ ) mpz_init(w[i]);

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

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

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

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

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;

Line 1703 get_eg(&tmp2); add_eg(&eg_mul2,&tmp1,&tmp2);

Line 1728 get_eg(&tmp2); add_eg(&eg_mul2,&tmp1,&tmp2);

if ( !cinfo[j] )

cind[k++] = j;

return rank;

}

} else

goto reset;

} else {

reset:

period = period*3/2;

count = 0;

}

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

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

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

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

for ( i = 0; i < rank; i++ )

for ( j = 0, bi = b[i], wi = wc[i]; j < ri; j++ )

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

for ( i = 0; i < rank; i++ )

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 1856 int generic_gauss_elim_hensel_dalg(MAT mat,DP *mb,MAT

mpz_init_set(uz,BDY(b[i][j]));

else

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

}

MPZTOZ(uz,u);

divsz(u,mdq,&b[i][j]);

}

count++;

/* q = q*md */

mulz(q,mdq,&u); q = u;

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

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

}

Line 2014 RESET:

Line 2042 RESET:

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,rind,cind) ) {

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

Line 2127 RESET:

Line 2155 RESET:

}

#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

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;

}

#endif

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