version 1.2, 2018/09/21 07:06:51 |
version 1.4, 2018/09/25 07:36:01 |
Line 25 void gc_free(void *p,size_t size) |
|
Line 25 void gc_free(void *p,size_t size) |
|
|
|
void init_gmpq() |
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); |
mpz_init(ONEMPZ); mpz_set_ui(ONEMPZ,1); MPZTOZ(ONEMPZ,ONE); |
gmp_randinit_default(GMP_RAND); |
gmp_randinit_default(GMP_RAND); |
} |
} |
|
|
|
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,a[i][j]); |
|
printf("\n"); |
|
} |
|
|
Z utoz(unsigned int u) |
Z utoz(unsigned int u) |
{ |
{ |
mpz_t z; |
mpz_t z; |
Line 161 void mulz(Z n1,Z n2,Z *nr) |
|
Line 171 void mulz(Z n1,Z n2,Z *nr) |
|
|
|
void muladdtoz(Z n1,Z n2,Z *nr) |
void muladdtoz(Z n1,Z n2,Z *nr) |
{ |
{ |
#if 1 |
#if 0 |
Z t; |
Z t; |
|
|
if ( n1 && n2 ) { |
if ( n1 && n2 ) { |
Line 171 void muladdtoz(Z n1,Z n2,Z *nr) |
|
Line 181 void muladdtoz(Z n1,Z n2,Z *nr) |
|
mpz_addmul(BDY(*nr),BDY(n1),BDY(n2)); |
mpz_addmul(BDY(*nr),BDY(n1),BDY(n2)); |
if ( !mpz_sgn(BDY(*nr)) ) |
if ( !mpz_sgn(BDY(*nr)) ) |
*nr = 0; |
*nr = 0; |
} |
} |
#else |
#else |
Z t,s; |
Z t,s; |
|
|
Line 183 void muladdtoz(Z n1,Z n2,Z *nr) |
|
Line 193 void muladdtoz(Z n1,Z n2,Z *nr) |
|
|
|
void mul1addtoz(Z n1,long u,Z *nr) |
void mul1addtoz(Z n1,long u,Z *nr) |
{ |
{ |
|
#if 0 |
Z t; |
Z t; |
|
|
if ( n1 && u ) { |
if ( n1 && u ) { |
Line 196 void mul1addtoz(Z n1,long u,Z *nr) |
|
Line 207 void mul1addtoz(Z n1,long u,Z *nr) |
|
if ( !mpz_sgn(BDY(*nr)) ) |
if ( !mpz_sgn(BDY(*nr)) ) |
*nr = 0; |
*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) |
void mul1z(Z n1,long n2,Z *nr) |
Line 437 void gcdvz_estimate(VECT v,Z *q) |
|
Line 453 void gcdvz_estimate(VECT v,Z *q) |
|
else addz(s,b[i],&u); |
else addz(s,b[i],&u); |
s = 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); |
if ( b[i] && mpz_sgn(BDY(b[i]))<0 ) subz(t,b[i],&u); |
else addz(t,b[i],&u); |
else addz(t,b[i],&u); |
t = u; |
t = u; |
Line 445 void gcdvz_estimate(VECT v,Z *q) |
|
Line 461 void gcdvz_estimate(VECT v,Z *q) |
|
gcdz(s,t,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) |
void factorialz(unsigned int n,Z *nr) |
{ |
{ |
mpz_t a; |
mpz_t a; |
Line 1421 int generic_gauss_elim_hensel(MAT mat,MAT *nmmat,Z *dn |
|
Line 1480 int generic_gauss_elim_hensel(MAT mat,MAT *nmmat,Z *dn |
|
int *rind,*cind; |
int *rind,*cind; |
int count; |
int count; |
int ret; |
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 period; |
int *wx,*ptr; |
int *wx,*ptr; |
int wxsize,nsize; |
int wxsize,nsize; |
Z wn; |
Z wn; |
Z wq; |
Z wq; |
|
|
|
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; |
w = (int **)almat(row,col); |
w = (int **)almat(row,col); |
Line 1478 int generic_gauss_elim_hensel(MAT mat,MAT *nmmat,Z *dn |
|
Line 1538 int generic_gauss_elim_hensel(MAT mat,MAT *nmmat,Z *dn |
|
MKMAT(xmat,rank,ri); x = (Z **)(xmat)->body; |
MKMAT(xmat,rank,ri); x = (Z **)(xmat)->body; |
MKMAT(*nmmat,rank,ri); nm = (Z **)(*nmmat)->body; |
MKMAT(*nmmat,rank,ri); nm = (Z **)(*nmmat)->body; |
wc = (int **)almat(rank,ri); |
wc = (int **)almat(rank,ri); |
for ( i = 0; i < rank; i++ ) |
|
wc[i] = w[i]+rank; |
|
*rindp = rind = (int *)MALLOC_ATOMIC(rank*sizeof(int)); |
*rindp = rind = (int *)MALLOC_ATOMIC(rank*sizeof(int)); |
*cindp = cind = (int *)MALLOC_ATOMIC((ri)*sizeof(int)); |
*cindp = cind = (int *)MALLOC_ATOMIC((ri)*sizeof(int)); |
|
|
period = F4_INTRAT_PERIOD; |
period = F4_INTRAT_PERIOD; |
for ( q = ONE, count = 0; ; ) { |
for ( q = ONE, count = 0; ; ) { |
|
/* check Ax=B mod q */ |
if ( DP_Print > 3 ) |
if ( DP_Print > 3 ) |
fprintf(stderr,"o"); |
fprintf(stderr,"o"); |
/* wc = b mod md */ |
/* wc = b mod md */ |
for ( i = 0; i < rank; i++ ) |
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); |
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 */ |
/* x += q*wc */ |
|
get_eg(&tmp0); |
for ( i = 0; i < rank; i++ ) |
for ( i = 0; i < rank; i++ ) |
for ( j = 0, wi = wc[i]; j < ri; j++ ) mul1addtoz(q,wi[j],&x[i][j]); |
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 ( i = 0; i < rank; i++ ) |
for ( j = 0; j < ri; j++ ) { |
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]); |
divsz(u,mdq,&b[i][j]); |
} |
} |
|
get_eg(&tmp2); add_eg(&eg_mul2,&tmp1,&tmp2); |
count++; |
count++; |
/* q = q*md */ |
/* q = q*md */ |
mulz(q,mdq,&u); q = u; |
mulz(q,mdq,&u); q = u; |
if ( count == period ) { |
if ( count == period ) { |
ret = intmtoratm(xmat,q,*nmmat,dn); |
ret = intmtoratm(xmat,q,*nmmat,dn); |
if ( ret ) { |
if ( ret ) { |
|
print_eg("MUL1",&eg_mul1); |
|
print_eg("MUL2",&eg_mul2); |
for ( j = k = l = 0; j < col; j++ ) |
for ( j = k = l = 0; j < col; j++ ) |
if ( cinfo[j] ) |
if ( cinfo[j] ) |
rind[k++] = j; |
rind[k++] = j; |
Line 1593 int generic_gauss_elim_hensel_dalg(MAT mat,DP *mb,MAT |
|
Line 1667 int generic_gauss_elim_hensel_dalg(MAT mat,DP *mb,MAT |
|
a = (Z **)almat_pointer(rank,rank); /* lhs mat */ |
a = (Z **)almat_pointer(rank,rank); /* lhs mat */ |
MKMAT(bmat,rank,col-rank); b = (Z **)bmat->body; /* lhs mat */ |
MKMAT(bmat,rank,col-rank); b = (Z **)bmat->body; /* lhs mat */ |
for ( j = li = ri = 0; j < col; j++ ) |
for ( j = li = ri = 0; j < col; j++ ) |
if ( cinfo[j] ) { |
if ( cinfo[j] > 0 ) { |
/* the column is in lhs */ |
/* the column is in lhs */ |
for ( i = 0; i < rank; i++ ) { |
for ( i = 0; i < rank; i++ ) { |
w[i][li] = w[i][j]; |
w[i][li] = w[i][j]; |
a[i][li] = a0[rinfo[i]][j]; |
a[i][li] = a0[rinfo[i]][j]; |
} |
} |
li++; |
li++; |
} else { |
} else if ( !cinfo[j] ) { |
/* the column is in rhs */ |
/* the column is in rhs */ |
for ( i = 0; i < rank; i++ ) |
for ( i = 0; i < rank; i++ ) |
b[i][ri] = a0[rinfo[i]][j]; |
b[i][ri] = a0[rinfo[i]][j]; |
Line 1623 int generic_gauss_elim_hensel_dalg(MAT mat,DP *mb,MAT |
|
Line 1697 int generic_gauss_elim_hensel_dalg(MAT mat,DP *mb,MAT |
|
MKMAT(xmat,rank,ri); x = (Z **)(xmat)->body; |
MKMAT(xmat,rank,ri); x = (Z **)(xmat)->body; |
MKMAT(*nmmat,rank,ri); nm = (Z **)(*nmmat)->body; |
MKMAT(*nmmat,rank,ri); nm = (Z **)(*nmmat)->body; |
wc = (int **)almat(rank,ri); |
wc = (int **)almat(rank,ri); |
for ( i = 0; i < rank; i++ ) |
|
wc[i] = w[i]+rank; |
|
*rindp = rind = (int *)MALLOC_ATOMIC(rank*sizeof(int)); |
*rindp = rind = (int *)MALLOC_ATOMIC(rank*sizeof(int)); |
*cindp = cind = (int *)MALLOC_ATOMIC((ri)*sizeof(int)); |
*cindp = cind = (int *)MALLOC_ATOMIC((ri)*sizeof(int)); |
|
|
Line 1634 int generic_gauss_elim_hensel_dalg(MAT mat,DP *mb,MAT |
|
Line 1706 int generic_gauss_elim_hensel_dalg(MAT mat,DP *mb,MAT |
|
fprintf(stderr,"o"); |
fprintf(stderr,"o"); |
/* wc = b mod md */ |
/* wc = b mod md */ |
for ( i = 0; i < rank; i++ ) |
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); |
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 */ |
/* wc = A^(-1)wc; wc is normalized */ |
solve_by_lu_mod(w,rank,md,wc,ri,1); |
solve_by_lu_mod(w,rank,md,wc,ri,1); |
/* x += q*wc */ |
/* x += q*wc */ |
for ( i = 0; i < rank; i++ ) |
for ( i = 0; i < rank; i++ ) |
for ( j = 0, wi = wc[i]; j < ri; j++ ) mul1addtoz(q,wi[j],&x[i][j]); |
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 ( i = 0; i < rank; i++ ) |
for ( j = 0; j < ri; j++ ) { |
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]); |
divsz(u,mdq,&b[i][j]); |
} |
} |
count++; |
count++; |