version 1.10, 2018/12/10 22:24:42 |
version 1.11, 2019/03/03 05:21:17 |
|
|
/* $OpenXM: OpenXM_contrib2/asir2018/engine/Q.c,v 1.9 2018/10/19 23:27:38 noro Exp $ */ |
/* $OpenXM: OpenXM_contrib2/asir2018/engine/Q.c,v 1.10 2018/12/10 22:24:42 noro Exp $ */ |
#include "ca.h" |
#include "ca.h" |
#include "gmp.h" |
#include "gmp.h" |
#include "base.h" |
#include "base.h" |
Line 1229 int intmtoratm(MAT mat,Z md,MAT nm,Z *dn) |
|
Line 1229 int intmtoratm(MAT mat,Z md,MAT nm,Z *dn) |
|
return 0; |
return 0; |
row = mat->row; col = mat->col; |
row = mat->row; col = mat->col; |
bshiftz(md,1,&t); |
bshiftz(md,1,&t); |
isqrtz(t,&s); |
isqrt(t,&s); |
bshiftz(s,64,&b); |
bshiftz(s,64,&b); |
if ( !b ) b = ONE; |
if ( !b ) b = ONE; |
dn0 = ONE; |
dn0 = ONE; |
Line 1469 void isqrtz(Z a,Z *r) |
|
Line 1469 void isqrtz(Z a,Z *r) |
|
Z two; |
Z two; |
|
|
if ( !a ) *r = 0; |
if ( !a ) *r = 0; |
else if ( UNIQ(a) ) *r = ONE; |
else if ( UNIZ(a) ) *r = ONE; |
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 */ |
Line 1592 int generic_gauss_elim_hensel(MAT mat,MAT *nmmat,Z *dn |
|
Line 1592 int generic_gauss_elim_hensel(MAT mat,MAT *nmmat,Z *dn |
|
Z wq; |
Z wq; |
|
|
#if SIZEOF_LONG == 8 |
#if SIZEOF_LONG == 8 |
return generic_gauss_elim_hensel64(mat,nmmat,dn,rindp,cindp); |
return generic_gauss_elim_hensel64(mat,nmmat,dn,rindp,cindp,0); |
#endif |
#endif |
init_eg(&eg_mul1); init_eg(&eg_mul2); |
init_eg(&eg_mul1); init_eg(&eg_mul2); |
a0 = (Z **)mat->body; |
a0 = (Z **)mat->body; |
Line 1706 get_eg(&tmp2); add_eg(&eg_mul2,&tmp1,&tmp2); |
|
Line 1706 get_eg(&tmp2); add_eg(&eg_mul2,&tmp1,&tmp2); |
|
if ( !cinfo[j] ) |
if ( !cinfo[j] ) |
cind[k++] = j; |
cind[k++] = j; |
return rank; |
return rank; |
} |
} else |
|
goto reset; |
} else { |
} else { |
|
reset: |
period = period*3/2; |
period = period*3/2; |
count = 0; |
count = 0; |
} |
} |
Line 1743 int generic_gauss_elim_hensel_dalg(MAT mat,DP *mb,MAT |
|
Line 1745 int generic_gauss_elim_hensel_dalg(MAT mat,DP *mb,MAT |
|
Z wq; |
Z wq; |
DP m; |
DP m; |
|
|
|
#if SIZEOF_LONG == 8 |
|
return generic_gauss_elim_hensel64(mat,nmmat,dn,rindp,cindp,mb); |
|
#endif |
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 1815 int generic_gauss_elim_hensel_dalg(MAT mat,DP *mb,MAT |
|
Line 1820 int generic_gauss_elim_hensel_dalg(MAT mat,DP *mb,MAT |
|
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); |
/* 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 */ |
/* 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]); |
Line 1829 int generic_gauss_elim_hensel_dalg(MAT mat,DP *mb,MAT |
|
Line 1834 int generic_gauss_elim_hensel_dalg(MAT mat,DP *mb,MAT |
|
mpz_init_set(uz,BDY(b[i][j])); |
mpz_init_set(uz,BDY(b[i][j])); |
else |
else |
mpz_init_set_ui(uz,0); |
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]); |
mpz_submul_ui(uz,BDY(a[i][k]),wc[k][j]); |
} |
|
} |
|
MPZTOZ(uz,u); |
MPZTOZ(uz,u); |
divsz(u,mdq,&b[i][j]); |
divsz(u,mdq,&b[i][j]); |
} |
} |
|
|
count++; |
count++; |
/* q = q*md */ |
/* q = q*md */ |
mulz(q,mdq,&u); q = u; |
mulz(q,mdq,&u); q = u; |
Line 1859 int generic_gauss_elim_hensel_dalg(MAT mat,DP *mb,MAT |
|
Line 1860 int generic_gauss_elim_hensel_dalg(MAT mat,DP *mb,MAT |
|
if ( !cinfo[j] ) |
if ( !cinfo[j] ) |
cind[k++] = j; |
cind[k++] = j; |
return rank; |
return rank; |
} |
} else |
|
goto reset; |
} else { |
} else { |
|
reset: |
period = period*3/2; |
period = period*3/2; |
count = 0; |
count = 0; |
} |
} |
|
|
} |
} |
#endif |
#endif |
|
|
int generic_gauss_elim_hensel64(MAT mat,MAT *nmmat,Z *dn,int **rindp,int **cindp) |
int generic_gauss_elim_hensel64(MAT mat,MAT *nmmat,Z *dn,int **rindp,int **cindp,DP *mb) |
{ |
{ |
MAT r; |
MAT r; |
Z z; |
Z z; |
Line 2151 int generic_gauss_elim_hensel64(MAT mat,MAT *nmmat,Z * |
|
Line 2154 int generic_gauss_elim_hensel64(MAT mat,MAT *nmmat,Z * |
|
int count; |
int count; |
int ret; |
int ret; |
int period; |
int period; |
|
DP m; |
|
|
a0 = (Z **)mat->body; |
a0 = (Z **)mat->body; |
row = mat->row; col = mat->col; |
row = mat->row; col = mat->col; |
Line 2168 int generic_gauss_elim_hensel64(MAT mat,MAT *nmmat,Z * |
|
Line 2172 int generic_gauss_elim_hensel64(MAT mat,MAT *nmmat,Z * |
|
if ( DP_Print > 3 ) { |
if ( DP_Print > 3 ) { |
fprintf(asir_out,"done.\n"); fflush(asir_out); |
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 */ |
a = (mpz_t **)mpz_allocmat(rank,rank); /* lhs mat */ |
b = (mpz_t **)mpz_allocmat(rank,col-rank); |
b = (mpz_t **)mpz_allocmat(rank,col-rank); |
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]; |
Line 2181 int generic_gauss_elim_hensel64(MAT mat,MAT *nmmat,Z * |
|
Line 2198 int generic_gauss_elim_hensel64(MAT mat,MAT *nmmat,Z * |
|
mpz_set_ui(a[i][li],0); |
mpz_set_ui(a[i][li],0); |
} |
} |
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++ ) { |
if ( a0[rinfo[i]][j] ) |
if ( a0[rinfo[i]][j] ) |
Line 2249 int generic_gauss_elim_hensel64(MAT mat,MAT *nmmat,Z * |
|
Line 2266 int generic_gauss_elim_hensel64(MAT mat,MAT *nmmat,Z * |
|
count++; |
count++; |
/* q = q*md */ |
/* q = q*md */ |
mpz_mul_ui(q,q,md); |
mpz_mul_ui(q,q,md); |
fprintf(stderr,"."); |
|
if ( count == period ) { |
if ( count == period ) { |
ret = mpz_intmtoratm(x,rank,ri,q,nm,den); |
ret = mpz_intmtoratm(x,rank,ri,q,nm,den); |
if ( ret ) { |
if ( ret ) { |
for ( j = k = l = 0; j < col; j++ ) |
for ( j = k = l = 0; j < col; j++ ) |
if ( cinfo[j] ) |
if ( cinfo[j] > 0 ) |
rind[k++] = j; |
rind[k++] = j; |
else |
else if ( !cinfo[j] ) |
cind[l++] = j; |
cind[l++] = j; |
ret = mpz_gensolve_check(mat,nm,den,rank,rind,cind); |
ret = mpz_gensolve_check(mat,nm,den,rank,rind,cind); |
if ( ret ) { |
if ( ret ) { |
Line 2272 int generic_gauss_elim_hensel64(MAT mat,MAT *nmmat,Z * |
|
Line 2288 int generic_gauss_elim_hensel64(MAT mat,MAT *nmmat,Z * |
|
} |
} |
MPZTOZ(den,*dn); |
MPZTOZ(den,*dn); |
return rank; |
return rank; |
} |
} else |
|
goto reset; |
} else { |
} else { |
|
reset: |
fprintf(stderr,"F"); |
fprintf(stderr,"F"); |
period = period*3/2; |
period = period*3/2; |
count = 0; |
count = 0; |