version 1.4, 2018/09/25 07:36:01 |
version 1.6, 2018/10/01 05:49:06 |
|
|
|
/* $OpenXM$ */ |
#include "ca.h" |
#include "ca.h" |
#include "gmp.h" |
#include "gmp.h" |
#include "base.h" |
#include "base.h" |
Line 10 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 |
|
|
|
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); |
|
|
Line 353 void pwrz(Z n1,Z n,Z *nr) |
|
Line 358 void pwrz(Z n1,Z n,Z *nr) |
|
} else if ( !smallz(n) ) { |
} else if ( !smallz(n) ) { |
error("exponent too big."); *nr = 0; |
error("exponent too big."); *nr = 0; |
} else if ( n1->z && mpz_sgn(BDY((Z)n))>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 { |
} else { |
MPZTOMPQ(BDY(n1),q); MPQTOQ(q,r); |
MPZTOMPQ(BDY(n1),q); MPQTOQ(q,r); |
pwrq(r,(Q)n,&p); *nr = (Z)p; |
pwrq(r,(Q)n,&p); *nr = (Z)p; |
Line 648 void pwrq(Q n1,Q n,Q *nr) |
|
Line 653 void pwrq(Q n1,Q n,Q *nr) |
|
} else if ( !smallz((Z)n) ) { |
} else if ( !smallz((Z)n) ) { |
error("exponent too big."); *nr = 0; |
error("exponent too big."); *nr = 0; |
} else { |
} else { |
e = QTOS(n); |
e = ZTOS(n); |
if ( e < 0 ) { |
if ( e < 0 ) { |
e = -e; |
e = -e; |
if ( n1->z ) { |
if ( n1->z ) { |
Line 705 void mkbc(int n,Z *t) |
|
Line 710 void mkbc(int n,Z *t) |
|
Z c,d,iq; |
Z c,d,iq; |
|
|
for ( t[0] = ONE, i = 1; i <= n/2; i++ ) { |
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++ ) |
for ( ; i <= n; i++ ) |
t[i] = t[n-i]; |
t[i] = t[n-i]; |
Line 884 unsigned int remqi(Q a,unsigned int mod) |
|
Line 889 unsigned int remqi(Q a,unsigned int mod) |
|
return c; |
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 generic_gauss_elim(MAT mat,MAT *nm,Z *dn,int **rindp,int **cindp) |
{ |
{ |
int **wmat; |
int **wmat; |
Line 898 int generic_gauss_elim(MAT mat,MAT *nm,Z *dn,int **rin |
|
Line 899 int generic_gauss_elim(MAT mat,MAT *nm,Z *dn,int **rin |
|
MAT r,crmat; |
MAT r,crmat; |
int ret; |
int ret; |
|
|
|
#if SIZEOF_LONG == 8 |
|
return generic_gauss_elim64(mat,nm,dn,rindp,cindp); |
|
#endif |
bmat = (Z **)mat->body; |
bmat = (Z **)mat->body; |
row = mat->row; col = mat->col; |
row = mat->row; col = mat->col; |
wmat = (int **)almat(row,col); |
wmat = (int **)almat(row,col); |
Line 1369 void isqrtz(Z a,Z *r) |
|
Line 1373 void isqrtz(Z a,Z *r) |
|
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 */ |
STOQ(2,two); |
STOZ(2,two); |
while ( 1 ) { |
while ( 1 ) { |
pwrz(x,two,&t); |
pwrz(x,two,&t); |
if ( cmpz(t,a) <= 0 ) { |
if ( cmpz(t,a) <= 0 ) { |
Line 1493 init_eg(&eg_mul1); init_eg(&eg_mul2); |
|
Line 1497 init_eg(&eg_mul1); init_eg(&eg_mul2); |
|
w = (int **)almat(row,col); |
w = (int **)almat(row,col); |
for ( ind = 0; ; ind++ ) { |
for ( ind = 0; ; ind++ ) { |
md = get_lprime(ind); |
md = get_lprime(ind); |
STOQ(md,mdq); |
STOZ(md,mdq); |
for ( i = 0; i < row; i++ ) |
for ( i = 0; i < row; i++ ) |
for ( j = 0, ai = a0[i], wi = w[i]; j < col; j++ ) |
for ( j = 0, ai = a0[i], wi = w[i]; j < col; j++ ) |
wi[j] = remqi((Q)ai[j],md); |
wi[j] = remqi((Q)ai[j],md); |
Line 1641 int generic_gauss_elim_hensel_dalg(MAT mat,DP *mb,MAT |
|
Line 1645 int generic_gauss_elim_hensel_dalg(MAT mat,DP *mb,MAT |
|
w = (int **)almat(row,col); |
w = (int **)almat(row,col); |
for ( ind = 0; ; ind++ ) { |
for ( ind = 0; ; ind++ ) { |
md = get_lprime(ind); |
md = get_lprime(ind); |
STOQ(md,mdq); |
STOZ(md,mdq); |
for ( i = 0; i < row; i++ ) |
for ( i = 0; i < row; i++ ) |
for ( j = 0, ai = a0[i], wi = w[i]; j < col; j++ ) |
for ( j = 0, ai = a0[i], wi = w[i]; j < col; j++ ) |
wi[j] = remqi((Q)ai[j],md); |
wi[j] = remqi((Q)ai[j],md); |
Line 1761 int generic_gauss_elim_hensel_dalg(MAT mat,DP *mb,MAT |
|
Line 1765 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 |