version 1.5, 2017/01/08 03:05:39 |
version 1.11, 2018/07/28 00:45:55 |
|
|
GZ ONEGZ; |
GZ ONEGZ; |
int lf_lazy; |
int lf_lazy; |
GZ current_mod_lf; |
GZ current_mod_lf; |
|
int current_mod_lf_size; |
|
|
void isqrtgz(GZ a,GZ *r); |
void isqrtgz(GZ a,GZ *r); |
void bshiftgz(GZ a,int n,GZ *r); |
void bshiftgz(GZ a,int n,GZ *r); |
|
|
void *gc_realloc(void *p,size_t osize,size_t nsize) |
void *gc_realloc(void *p,size_t osize,size_t nsize) |
{ |
{ |
return (void *)Risa_GC_realloc(p,nsize); |
return (void *)Risa_GC_realloc(p,nsize); |
} |
} |
|
|
void gc_free(void *p,size_t size) |
void gc_free(void *p,size_t size) |
{ |
{ |
Risa_GC_free(p); |
Risa_GC_free(p); |
} |
} |
|
|
void init_gmpq() |
void init_gmpq() |
{ |
{ |
mp_set_memory_functions(Risa_GC_malloc,gc_realloc,gc_free); |
mp_set_memory_functions(Risa_GC_malloc,gc_realloc,gc_free); |
|
|
mpz_init(ONEMPZ); mpz_set_ui(ONEMPZ,1); MPZTOGZ(ONEMPZ,ONEGZ); |
mpz_init(ONEMPZ); mpz_set_ui(ONEMPZ,1); MPZTOGZ(ONEMPZ,ONEGZ); |
} |
} |
|
|
GZ utogz(unsigned int u) |
GZ utogz(unsigned int u) |
{ |
{ |
mpz_t z; |
mpz_t z; |
GZ r; |
GZ r; |
|
|
if ( !u ) return 0; |
if ( !u ) return 0; |
mpz_init(z); mpz_set_ui(z,u); MPZTOGZ(z,r); return r; |
mpz_init(z); mpz_set_ui(z,u); MPZTOGZ(z,r); return r; |
} |
} |
|
|
GZ stogz(int s) |
GZ stogz(int s) |
{ |
{ |
mpz_t z; |
mpz_t z; |
GZ r; |
GZ r; |
|
|
if ( !s ) return 0; |
if ( !s ) return 0; |
mpz_init(z); mpz_set_si(z,s); MPZTOGZ(z,r); return r; |
mpz_init(z); mpz_set_si(z,s); MPZTOGZ(z,r); return r; |
} |
} |
|
|
GQ mpqtogzq(mpq_t a) |
GQ mpqtogzq(mpq_t a) |
{ |
{ |
GZ z; |
GZ z; |
GQ q; |
GQ q; |
|
|
if ( INTMPQ(a) ) { |
if ( INTMPQ(a) ) { |
MPZTOGZ(mpq_numref(a),z); return (GQ)z; |
MPZTOGZ(mpq_numref(a),z); return (GQ)z; |
} else { |
} else { |
MPQTOGQ(a,q); return q; |
MPQTOGQ(a,q); return q; |
} |
} |
} |
} |
|
|
GZ ztogz(Q a) |
GZ ztogz(Q a) |
{ |
{ |
mpz_t z; |
mpz_t z; |
mpq_t b; |
mpq_t b; |
N nm; |
N nm; |
GZ s; |
GZ s; |
|
|
if ( !a ) return 0; |
if ( !a ) return 0; |
nm = NM(a); |
nm = NM(a); |
mpz_init(z); |
mpz_init(z); |
mpz_import(z,PL(nm),-1,sizeof(BD(nm)[0]),0,0,BD(nm)); |
mpz_import(z,PL(nm),-1,sizeof(BD(nm)[0]),0,0,BD(nm)); |
if ( SGN(a)<0 ) mpz_neg(z,z); |
if ( SGN(a)<0 ) mpz_neg(z,z); |
MPZTOGZ(z,s); |
MPZTOGZ(z,s); |
return s; |
return s; |
} |
} |
|
|
Q gztoz(GZ a) |
Q gztoz(GZ a) |
{ |
{ |
N nm; |
N nm; |
Q q; |
Q q; |
int sgn; |
int sgn; |
size_t len; |
size_t len; |
|
|
if ( !a ) return 0; |
if ( !a ) return 0; |
len = WORDSIZE_IN_N(BDY((GZ)a)); nm = NALLOC(len); |
len = WORDSIZE_IN_N(BDY((GZ)a)); nm = NALLOC(len); |
mpz_export(BD(nm),&len,-1,sizeof(int),0,0,BDY((GZ)a)); |
mpz_export(BD(nm),&len,-1,sizeof(int),0,0,BDY((GZ)a)); |
PL(nm) = len; |
PL(nm) = len; |
sgn = mpz_sgn(BDY((GZ)a)); NTOQ(nm,sgn,q); |
sgn = mpz_sgn(BDY((GZ)a)); NTOQ(nm,sgn,q); |
return q; |
return q; |
} |
} |
|
|
void dupgz(GZ a,GZ *b) |
void dupgz(GZ a,GZ *b) |
Line 101 void dupgz(GZ a,GZ *b) |
|
Line 102 void dupgz(GZ a,GZ *b) |
|
|
|
int n_bits_gz(GZ a) |
int n_bits_gz(GZ a) |
{ |
{ |
return a ? mpz_sizeinbase(BDY(a),2) : 0; |
return a ? mpz_sizeinbase(BDY(a),2) : 0; |
} |
} |
|
|
GQ qtogq(Q a) |
GQ qtogq(Q a) |
{ |
{ |
mpz_t z; |
mpz_t z; |
mpq_t b; |
mpq_t b; |
N nm,dn; |
N nm,dn; |
GZ s; |
GZ s; |
GQ r; |
GQ r; |
|
|
if ( !a ) return 0; |
if ( !a ) return 0; |
if ( INT(a) ) { |
if ( INT(a) ) { |
nm = NM(a); |
nm = NM(a); |
mpz_init(z); |
mpz_init(z); |
mpz_import(z,PL(nm),-1,sizeof(BD(nm)[0]),0,0,BD(nm)); |
mpz_import(z,PL(nm),-1,sizeof(BD(nm)[0]),0,0,BD(nm)); |
if ( SGN(a)<0 ) mpz_neg(z,z); |
if ( SGN(a)<0 ) mpz_neg(z,z); |
MPZTOGZ(z,s); |
MPZTOGZ(z,s); |
return (GQ)s; |
return (GQ)s; |
} else { |
} else { |
nm = NM(a); dn = DN(a); |
nm = NM(a); dn = DN(a); |
mpq_init(b); |
mpq_init(b); |
mpz_import(mpq_numref(b),PL(nm),-1,sizeof(BD(nm)[0]),0,0,BD(nm)); |
mpz_import(mpq_numref(b),PL(nm),-1,sizeof(BD(nm)[0]),0,0,BD(nm)); |
mpz_import(mpq_denref(b),PL(dn),-1,sizeof(BD(dn)[0]),0,0,BD(dn)); |
mpz_import(mpq_denref(b),PL(dn),-1,sizeof(BD(dn)[0]),0,0,BD(dn)); |
if ( SGN(a)<0 ) mpq_neg(b,b); |
if ( SGN(a)<0 ) mpq_neg(b,b); |
MPQTOGQ(b,r); |
MPQTOGQ(b,r); |
return r; |
return r; |
} |
} |
} |
} |
|
|
Q gqtoq(GQ a) |
Q gqtoq(GQ a) |
{ |
{ |
N nm,dn; |
N nm,dn; |
Q q; |
Q q; |
int sgn; |
int sgn; |
size_t len; |
size_t len; |
|
|
if ( !a ) return 0; |
if ( !a ) return 0; |
if ( NID(a) == N_GZ ) { |
if ( NID(a) == N_GZ ) { |
len = WORDSIZE_IN_N(BDY((GZ)a)); nm = NALLOC(len); |
len = WORDSIZE_IN_N(BDY((GZ)a)); nm = NALLOC(len); |
mpz_export(BD(nm),&len,-1,sizeof(int),0,0,BDY((GZ)a)); |
mpz_export(BD(nm),&len,-1,sizeof(int),0,0,BDY((GZ)a)); |
PL(nm) = len; |
PL(nm) = len; |
sgn = mpz_sgn(BDY((GZ)a)); NTOQ(nm,sgn,q); |
sgn = mpz_sgn(BDY((GZ)a)); NTOQ(nm,sgn,q); |
} else { |
} else { |
len = WORDSIZE_IN_N(mpq_numref(BDY(a))); nm = NALLOC(len); |
len = WORDSIZE_IN_N(mpq_numref(BDY(a))); nm = NALLOC(len); |
mpz_export(BD(nm),&len,-1,sizeof(int),0,0,mpq_numref(BDY(a))); |
mpz_export(BD(nm),&len,-1,sizeof(int),0,0,mpq_numref(BDY(a))); |
PL(nm) = len; |
PL(nm) = len; |
len = WORDSIZE_IN_N(mpq_denref(BDY(a))); dn = NALLOC(len); |
len = WORDSIZE_IN_N(mpq_denref(BDY(a))); dn = NALLOC(len); |
mpz_export(BD(dn),&len,-1,sizeof(int),0,0,mpq_denref(BDY(a))); |
mpz_export(BD(dn),&len,-1,sizeof(int),0,0,mpq_denref(BDY(a))); |
PL(dn) = len; |
PL(dn) = len; |
sgn = mpz_sgn(mpq_numref(BDY(a))); NDTOQ(nm,dn,sgn,q); |
sgn = mpz_sgn(mpq_numref(BDY(a))); NDTOQ(nm,dn,sgn,q); |
} |
} |
return q; |
return q; |
} |
} |
|
|
P ptogp(P a) |
P ptogp(P a) |
{ |
{ |
DCP dc,dcr,dcr0; |
DCP dc,dcr,dcr0; |
P b; |
P b; |
|
|
if ( !a ) return 0; |
if ( !a ) return 0; |
if ( NUM(a) ) return (P)qtogq((Q)a); |
if ( NUM(a) ) return (P)qtogq((Q)a); |
for ( dc = DC(a), dcr0 = 0; dc; dc = NEXT(dc) ) { |
for ( dc = DC(a), dcr0 = 0; dc; dc = NEXT(dc) ) { |
NEXTDC(dcr0,dcr); DEG(dcr) = DEG(dc); COEF(dcr) = (P)ptogp(COEF(dc)); |
NEXTDC(dcr0,dcr); DEG(dcr) = DEG(dc); COEF(dcr) = (P)ptogp(COEF(dc)); |
} |
} |
NEXT(dcr) = 0; MKP(VR(a),dcr0,b); |
NEXT(dcr) = 0; MKP(VR(a),dcr0,b); |
return b; |
return b; |
} |
} |
|
|
P gptop(P a) |
P gptop(P a) |
{ |
{ |
DCP dc,dcr,dcr0; |
DCP dc,dcr,dcr0; |
P b; |
P b; |
|
|
if ( !a ) return 0; |
if ( !a ) return 0; |
if ( NUM(a) ) b = (P)gqtoq((GQ)a); |
if ( NUM(a) ) b = (P)gqtoq((GQ)a); |
else { |
else { |
for ( dc = DC(a), dcr0 = 0; dc; dc = NEXT(dc) ) { |
for ( dc = DC(a), dcr0 = 0; dc; dc = NEXT(dc) ) { |
NEXTDC(dcr0,dcr); DEG(dcr) = DEG(dc); |
NEXTDC(dcr0,dcr); DEG(dcr) = DEG(dc); |
COEF(dcr) = (P)gptop(COEF(dc)); |
COEF(dcr) = (P)gptop(COEF(dc)); |
} |
} |
NEXT(dcr) = 0; MKP(VR(a),dcr0,b); |
NEXT(dcr) = 0; MKP(VR(a),dcr0,b); |
} |
} |
return b; |
return b; |
} |
} |
|
|
void addgz(GZ n1,GZ n2,GZ *nr) |
void addgz(GZ n1,GZ n2,GZ *nr) |
{ |
{ |
mpz_t t; |
mpz_t t; |
int s1,s2; |
int s1,s2; |
|
|
if ( !n1 ) *nr = n2; |
if ( !n1 ) *nr = n2; |
else if ( !n2 ) *nr = n1; |
else if ( !n2 ) *nr = n1; |
else { |
else { |
mpz_init(t); mpz_add(t,BDY(n1),BDY(n2)); MPZTOGZ(t,*nr); |
mpz_init(t); mpz_add(t,BDY(n1),BDY(n2)); MPZTOGZ(t,*nr); |
} |
} |
} |
} |
|
|
void subgz(GZ n1,GZ n2,GZ *nr) |
void subgz(GZ n1,GZ n2,GZ *nr) |
{ |
{ |
mpz_t t; |
mpz_t t; |
|
|
if ( !n1 ) |
if ( !n1 ) |
if ( !n2 ) |
if ( !n2 ) |
*nr = 0; |
*nr = 0; |
else { |
else { |
t[0] = BDY(n2)[0]; mpz_neg(t,t); MPZTOGZ(t,*nr); |
t[0] = BDY(n2)[0]; mpz_neg(t,t); MPZTOGZ(t,*nr); |
} |
} |
else if ( !n2 ) |
else if ( !n2 ) |
*nr = n1; |
*nr = n1; |
else if ( n1 == n2 ) |
else if ( n1 == n2 ) |
*nr = 0; |
*nr = 0; |
else { |
else { |
mpz_init(t); mpz_sub(t,BDY(n1),BDY(n2)); MPZTOGZ(t,*nr); |
mpz_init(t); mpz_sub(t,BDY(n1),BDY(n2)); MPZTOGZ(t,*nr); |
} |
} |
} |
} |
|
|
void mulgz(GZ n1,GZ n2,GZ *nr) |
void mulgz(GZ n1,GZ n2,GZ *nr) |
{ |
{ |
mpz_t t; |
mpz_t t; |
|
|
if ( !n1 || !n2 ) *nr = 0; |
if ( !n1 || !n2 ) *nr = 0; |
#if 1 |
#if 1 |
else if ( UNIGZ(n1) ) *nr = n2; |
else if ( UNIGZ(n1) ) *nr = n2; |
else if ( UNIGZ(n2) ) *nr = n1; |
else if ( UNIGZ(n2) ) *nr = n1; |
else if ( MUNIGZ(n1) ) chsgngz(n2,nr); |
else if ( MUNIGZ(n1) ) chsgngz(n2,nr); |
else if ( MUNIGZ(n2) ) chsgngz(n1,nr); |
else if ( MUNIGZ(n2) ) chsgngz(n1,nr); |
#endif |
#endif |
else { |
else { |
mpz_init(t); mpz_mul(t,BDY(n1),BDY(n2)); MPZTOGZ(t,*nr); |
mpz_init(t); mpz_mul(t,BDY(n1),BDY(n2)); MPZTOGZ(t,*nr); |
} |
} |
} |
} |
|
|
/* nr += n1*n2 */ |
/* nr += n1*n2 */ |
Line 240 void muladdtogz(GZ n1,GZ n2,GZ *nr) |
|
Line 241 void muladdtogz(GZ n1,GZ n2,GZ *nr) |
|
{ |
{ |
GZ t; |
GZ t; |
|
|
if ( n1 && n2 ) { |
if ( n1 && n2 ) { |
if ( !(*nr) ) { |
if ( !(*nr) ) { |
NEWGZ(t); mpz_init(BDY(t)); *nr = t; |
NEWGZ(t); mpz_init(BDY(t)); *nr = t; |
} |
} |
Line 250 void muladdtogz(GZ n1,GZ n2,GZ *nr) |
|
Line 251 void muladdtogz(GZ n1,GZ n2,GZ *nr) |
|
|
|
void mul1gz(GZ n1,int n2,GZ *nr) |
void mul1gz(GZ n1,int n2,GZ *nr) |
{ |
{ |
mpz_t t; |
mpz_t t; |
|
|
if ( !n1 || !n2 ) *nr = 0; |
if ( !n1 || !n2 ) *nr = 0; |
else { |
else { |
mpz_init(t); mpz_mul_ui(t,BDY(n1),(long)n2); MPZTOGZ(t,*nr); |
mpz_init(t); mpz_mul_ui(t,BDY(n1),(long)n2); MPZTOGZ(t,*nr); |
} |
} |
} |
} |
|
|
void divgz(GZ n1,GZ n2,GZ *nq) |
void divgz(GZ n1,GZ n2,GZ *nq) |
{ |
{ |
mpz_t t; |
mpz_t t; |
mpq_t a, b, q; |
mpq_t a, b, q; |
|
|
if ( !n2 ) { |
if ( !n2 ) { |
error("division by 0"); |
error("division by 0"); |
*nq = 0; |
*nq = 0; |
} else if ( !n1 ) |
} else if ( !n1 ) |
*nq = 0; |
*nq = 0; |
else if ( n1 == n2 ) { |
else if ( n1 == n2 ) { |
mpz_init(t); mpz_set_ui(t,1); MPZTOGZ(t,*nq); |
mpz_init(t); mpz_set_ui(t,1); MPZTOGZ(t,*nq); |
} else { |
} else { |
MPZTOMPQ(BDY(n1),a); MPZTOMPQ(BDY(n2),b); |
MPZTOMPQ(BDY(n1),a); MPZTOMPQ(BDY(n2),b); |
mpq_init(q); mpq_div(q,a,b); *nq = (GZ)mpqtogzq(q); |
mpq_init(q); mpq_div(q,a,b); *nq = (GZ)mpqtogzq(q); |
} |
} |
} |
} |
|
|
void remgz(GZ n1,GZ n2,GZ *nr) |
void remgz(GZ n1,GZ n2,GZ *nr) |
{ |
{ |
mpz_t r; |
mpz_t r; |
|
|
if ( !n2 ) { |
if ( !n2 ) { |
error("division by 0"); |
error("division by 0"); |
*nr = 0; |
*nr = 0; |
} else if ( !n1 || n1 == n2 ) |
} else if ( !n1 || n1 == n2 ) |
*nr = 0; |
*nr = 0; |
else { |
else { |
mpz_init(r); |
mpz_init(r); |
mpz_mod(r,BDY(n1),BDY(n2)); |
mpz_mod(r,BDY(n1),BDY(n2)); |
if ( !mpz_sgn(r) ) *nr = 0; |
if ( !mpz_sgn(r) ) *nr = 0; |
else MPZTOGZ(r,*nr); |
else MPZTOGZ(r,*nr); |
} |
} |
} |
} |
|
|
void divqrgz(GZ n1,GZ n2,GZ *nq,GZ *nr) |
void divqrgz(GZ n1,GZ n2,GZ *nq,GZ *nr) |
{ |
{ |
mpz_t t, a, b, q, r; |
mpz_t t, a, b, q, r; |
|
|
if ( !n2 ) { |
if ( !n2 ) { |
error("division by 0"); |
error("division by 0"); |
*nq = 0; *nr = 0; |
*nq = 0; *nr = 0; |
} else if ( !n1 ) { |
} else if ( !n1 ) { |
*nq = 0; *nr = 0; |
*nq = 0; *nr = 0; |
} else if ( n1 == n2 ) { |
} else if ( n1 == n2 ) { |
mpz_init(t); mpz_set_ui(t,1); MPZTOGZ(t,*nq); *nr = 0; |
mpz_init(t); mpz_set_ui(t,1); MPZTOGZ(t,*nq); *nr = 0; |
} else { |
} else { |
mpz_init(q); mpz_init(r); |
mpz_init(q); mpz_init(r); |
mpz_fdiv_qr(q,r,BDY(n1),BDY(n2)); |
mpz_fdiv_qr(q,r,BDY(n1),BDY(n2)); |
if ( !mpz_sgn(q) ) *nq = 0; |
if ( !mpz_sgn(q) ) *nq = 0; |
else MPZTOGZ(q,*nq); |
else MPZTOGZ(q,*nq); |
if ( !mpz_sgn(r) ) *nr = 0; |
if ( !mpz_sgn(r) ) *nr = 0; |
else MPZTOGZ(r,*nr); |
else MPZTOGZ(r,*nr); |
} |
} |
} |
} |
|
|
void divsgz(GZ n1,GZ n2,GZ *nq) |
void divsgz(GZ n1,GZ n2,GZ *nq) |
{ |
{ |
mpz_t t; |
mpz_t t; |
mpq_t a, b, q; |
mpq_t a, b, q; |
|
|
if ( !n2 ) { |
if ( !n2 ) { |
error("division by 0"); |
error("division by 0"); |
*nq = 0; |
*nq = 0; |
} else if ( !n1 ) |
} else if ( !n1 ) |
*nq = 0; |
*nq = 0; |
else if ( n1 == n2 ) { |
else if ( n1 == n2 ) { |
mpz_init(t); mpz_set_ui(t,1); MPZTOGZ(t,*nq); |
mpz_init(t); mpz_set_ui(t,1); MPZTOGZ(t,*nq); |
} else { |
} else { |
mpz_init(t); mpz_divexact(t,BDY(n1),BDY(n2)); MPZTOGZ(t,*nq); |
mpz_init(t); mpz_divexact(t,BDY(n1),BDY(n2)); MPZTOGZ(t,*nq); |
} |
} |
} |
} |
|
|
void chsgngz(GZ n,GZ *nr) |
void chsgngz(GZ n,GZ *nr) |
{ |
{ |
mpz_t t; |
mpz_t t; |
|
|
if ( !n ) |
if ( !n ) |
*nr = 0; |
*nr = 0; |
else { |
else { |
t[0] = BDY(n)[0]; mpz_neg(t,t); MPZTOGZ(t,*nr); |
t[0] = BDY(n)[0]; mpz_neg(t,t); MPZTOGZ(t,*nr); |
} |
} |
} |
} |
|
|
void pwrgz(GZ n1,Q n,GZ *nr) |
void pwrgz(GZ n1,Q n,GZ *nr) |
{ |
{ |
mpq_t t,q; |
mpq_t t,q; |
mpz_t z; |
mpz_t z; |
GQ p,r; |
GQ p,r; |
|
|
if ( !n || UNIGZ(n1) ) *nr = ONEGZ; |
if ( !n || UNIGZ(n1) ) *nr = ONEGZ; |
else if ( !n1 ) *nr = 0; |
else if ( !n1 ) *nr = 0; |
else if ( !INT(n) ) { |
else if ( !INT(n) ) { |
error("can't calculate fractional power."); *nr = 0; |
error("can't calculate fractional power."); *nr = 0; |
} else if ( MUNIGZ(n1) ) *nr = BD(NM(n))[0]%2 ? n1 : ONEGZ; |
} else if ( MUNIGZ(n1) ) *nr = BD(NM(n))[0]%2 ? n1 : ONEGZ; |
else if ( PL(NM(n)) > 1 ) { |
else if ( PL(NM(n)) > 1 ) { |
error("exponent too big."); *nr = 0; |
error("exponent too big."); *nr = 0; |
} else if ( NID(n1)==N_GZ && SGN(n)>0 ) { |
} else if ( NID(n1)==N_GZ && SGN(n)>0 ) { |
mpz_init(z); mpz_pow_ui(z,BDY(n1),QTOS(n)); MPZTOGZ(z,*nr); |
mpz_init(z); mpz_pow_ui(z,BDY(n1),QTOS(n)); MPZTOGZ(z,*nr); |
} else { |
} else { |
MPZTOMPQ(BDY(n1),q); MPQTOGQ(q,r); |
MPZTOMPQ(BDY(n1),q); MPQTOGQ(q,r); |
pwrgq(r,n,&p); *nr = (GZ)p; |
pwrgq(r,n,&p); *nr = (GZ)p; |
} |
} |
} |
} |
|
|
int cmpgz(GZ q1,GZ q2) |
int cmpgz(GZ q1,GZ q2) |
{ |
{ |
int sgn; |
int sgn; |
|
|
if ( !q1 ) |
if ( !q1 ) |
if ( !q2 ) |
if ( !q2 ) |
return 0; |
return 0; |
else |
else |
return -mpz_sgn(BDY(q2)); |
return -mpz_sgn(BDY(q2)); |
else if ( !q2 ) |
else if ( !q2 ) |
return mpz_sgn(BDY(q1)); |
return mpz_sgn(BDY(q1)); |
else if ( (sgn = mpz_sgn(BDY(q1))) != mpz_sgn(BDY(q2)) ) |
else if ( (sgn = mpz_sgn(BDY(q1))) != mpz_sgn(BDY(q2)) ) |
return sgn; |
return sgn; |
else { |
else { |
sgn = mpz_cmp(BDY(q1),BDY(q2)); |
sgn = mpz_cmp(BDY(q1),BDY(q2)); |
if ( sgn > 0 ) return 1; |
if ( sgn > 0 ) return 1; |
else if ( sgn < 0 ) return -1; |
else if ( sgn < 0 ) return -1; |
else return 0; |
else return 0; |
} |
} |
} |
} |
|
|
void gcdgz(GZ n1,GZ n2,GZ *nq) |
void gcdgz(GZ n1,GZ n2,GZ *nq) |
{ |
{ |
mpz_t t; |
mpz_t t; |
|
|
if ( !n1 ) *nq = n2; |
if ( !n1 ) *nq = n2; |
else if ( !n2 ) *nq = n1; |
else if ( !n2 ) *nq = n1; |
else { |
else { |
mpz_init(t); mpz_gcd(t,BDY(n1),BDY(n2)); |
mpz_init(t); mpz_gcd(t,BDY(n1),BDY(n2)); |
MPZTOGZ(t,*nq); |
MPZTOGZ(t,*nq); |
} |
} |
} |
} |
|
|
void invgz(GZ n1,GZ *nq) |
void invgz(GZ n1,GZ *nq) |
{ |
{ |
mpz_t t; |
mpz_t t; |
|
|
mpz_init(t); mpz_invert(t,BDY(n1),BDY(current_mod_lf)); |
mpz_init(t); mpz_invert(t,BDY(n1),BDY(current_mod_lf)); |
MPZTOGZ(t,*nq); |
MPZTOGZ(t,*nq); |
} |
} |
|
|
void lcmgz(GZ n1,GZ n2,GZ *nq) |
void lcmgz(GZ n1,GZ n2,GZ *nq) |
{ |
{ |
GZ g,t; |
GZ g,t; |
|
|
if ( !n1 || !n2 ) *nq = 0; |
if ( !n1 || !n2 ) *nq = 0; |
else { |
else { |
gcdgz(n1,n2,&g); divsgz(n1,g,&t); |
gcdgz(n1,n2,&g); divsgz(n1,g,&t); |
mulgz(n2,t,nq); |
mulgz(n2,t,nq); |
} |
} |
} |
} |
|
|
void gcdvgz(VECT v,GZ *q) |
void gcdvgz(VECT v,GZ *q) |
{ |
{ |
int n,i; |
int n,i; |
GZ *b; |
GZ *b; |
GZ g,g1; |
GZ g,g1; |
|
|
n = v->len; |
n = v->len; |
b = (GZ *)v->body; |
b = (GZ *)v->body; |
g = b[0]; |
g = b[0]; |
for ( i = 1; i < n; i++ ) { |
for ( i = 1; i < n; i++ ) { |
gcdgz(g,b[i],&g1); g = g1; |
gcdgz(g,b[i],&g1); g = g1; |
} |
} |
*q = g; |
*q = g; |
} |
} |
|
|
void gcdvgz_estimate(VECT v,GZ *q) |
void gcdvgz_estimate(VECT v,GZ *q) |
{ |
{ |
int n,m,i; |
int n,m,i; |
GZ s,t,u; |
GZ s,t,u; |
GZ *b; |
GZ *b; |
|
|
n = v->len; |
n = v->len; |
b = (GZ *)v->body; |
b = (GZ *)v->body; |
if ( n == 1 ) { |
if ( n == 1 ) { |
if ( mpz_sgn(BDY(b[0]))<0 ) chsgngz(b[0],q); |
if ( mpz_sgn(BDY(b[0]))<0 ) chsgngz(b[0],q); |
else *q = b[0]; |
else *q = b[0]; |
} |
} |
m = n/2; |
m = n/2; |
for ( i = 0, s = 0; i < m; i++ ) { |
for ( i = 0, s = 0; i < m; i++ ) { |
if ( b[i] && mpz_sgn(BDY(b[i]))<0 ) subgz(s,b[i],&u); |
if ( b[i] && mpz_sgn(BDY(b[i]))<0 ) subgz(s,b[i],&u); |
else addgz(s,b[i],&u); |
else addgz(s,b[i],&u); |
s = u; |
s = u; |
} |
} |
for ( i = 0, t = 0; i < n; i++ ) { |
for ( i = 0, t = 0; i < n; i++ ) { |
if ( b[i] && mpz_sgn(BDY(b[i]))<0 ) subgz(t,b[i],&u); |
if ( b[i] && mpz_sgn(BDY(b[i]))<0 ) subgz(t,b[i],&u); |
else addgz(t,b[i],&u); |
else addgz(t,b[i],&u); |
t = u; |
t = u; |
} |
} |
gcdgz(s,t,q); |
gcdgz(s,t,q); |
} |
} |
|
|
void addgq(GQ n1,GQ n2,GQ *nr) |
void addgq(GQ n1,GQ n2,GQ *nr) |
{ |
{ |
mpq_t q1,q2,t; |
mpq_t q1,q2,t; |
|
|
if ( !n1 ) *nr = n2; |
if ( !n1 ) *nr = n2; |
else if ( !n2 ) *nr = n1; |
else if ( !n2 ) *nr = n1; |
else { |
else { |
if ( NID(n1) == N_GZ ) MPZTOMPQ(BDY((GZ)n1),q1); |
if ( NID(n1) == N_GZ ) MPZTOMPQ(BDY((GZ)n1),q1); |
else q1[0] = BDY(n1)[0]; |
else q1[0] = BDY(n1)[0]; |
if ( NID(n2) == N_GZ ) MPZTOMPQ(BDY((GZ)n2),q2); |
if ( NID(n2) == N_GZ ) MPZTOMPQ(BDY((GZ)n2),q2); |
else q2[0] = BDY(n2)[0]; |
else q2[0] = BDY(n2)[0]; |
mpq_init(t); mpq_add(t,q1,q2); *nr = mpqtogzq(t); |
mpq_init(t); mpq_add(t,q1,q2); *nr = mpqtogzq(t); |
} |
} |
} |
} |
|
|
void subgq(GQ n1,GQ n2,GQ *nr) |
void subgq(GQ n1,GQ n2,GQ *nr) |
{ |
{ |
mpq_t q1,q2,t; |
mpq_t q1,q2,t; |
|
|
if ( !n1 ) |
if ( !n1 ) |
if ( !n2 ) *nr = 0; |
if ( !n2 ) *nr = 0; |
else { |
else { |
if ( NID(n1) == N_GZ ) chsgngz((GZ)n1,(GZ *)nr); |
if ( NID(n1) == N_GZ ) chsgngz((GZ)n1,(GZ *)nr); |
else { |
else { |
mpq_init(t); mpq_neg(t,BDY(n2)); MPQTOGQ(t,*nr); |
mpq_init(t); mpq_neg(t,BDY(n2)); MPQTOGQ(t,*nr); |
} |
} |
} |
} |
else if ( !n2 ) *nr = n1; |
else if ( !n2 ) *nr = n1; |
else if ( n1 == n2 ) *nr = 0; |
else if ( n1 == n2 ) *nr = 0; |
else { |
else { |
if ( NID(n1) == N_GZ ) MPZTOMPQ(BDY((GZ)n1),q1); |
if ( NID(n1) == N_GZ ) MPZTOMPQ(BDY((GZ)n1),q1); |
else q1[0] = BDY(n1)[0]; |
else q1[0] = BDY(n1)[0]; |
if ( NID(n2) == N_GZ ) MPZTOMPQ(BDY((GZ)n2),q2); |
if ( NID(n2) == N_GZ ) MPZTOMPQ(BDY((GZ)n2),q2); |
else q2[0] = BDY(n2)[0]; |
else q2[0] = BDY(n2)[0]; |
mpq_init(t); mpq_sub(t,q1,q2); *nr = mpqtogzq(t); |
mpq_init(t); mpq_sub(t,q1,q2); *nr = mpqtogzq(t); |
} |
} |
} |
} |
|
|
void mulgq(GQ n1,GQ n2,GQ *nr) |
void mulgq(GQ n1,GQ n2,GQ *nr) |
{ |
{ |
mpq_t t,q1,q2; |
mpq_t t,q1,q2; |
|
|
if ( !n1 || !n2 ) *nr = 0; |
if ( !n1 || !n2 ) *nr = 0; |
else { |
else { |
if ( NID(n1) == N_GZ ) MPZTOMPQ(BDY((GZ)n1),q1); |
if ( NID(n1) == N_GZ ) MPZTOMPQ(BDY((GZ)n1),q1); |
else q1[0] = BDY(n1)[0]; |
else q1[0] = BDY(n1)[0]; |
if ( NID(n2) == N_GZ ) MPZTOMPQ(BDY((GZ)n2),q2); |
if ( NID(n2) == N_GZ ) MPZTOMPQ(BDY((GZ)n2),q2); |
else q2[0] = BDY(n2)[0]; |
else q2[0] = BDY(n2)[0]; |
mpq_init(t); mpq_mul(t,q1,q2); *nr = mpqtogzq(t); |
mpq_init(t); mpq_mul(t,q1,q2); *nr = mpqtogzq(t); |
} |
} |
} |
} |
|
|
void divgq(GQ n1,GQ n2,GQ *nq) |
void divgq(GQ n1,GQ n2,GQ *nq) |
{ |
{ |
mpq_t t,q1,q2; |
mpq_t t,q1,q2; |
|
|
if ( !n2 ) { |
if ( !n2 ) { |
error("division by 0"); |
error("division by 0"); |
*nq = 0; |
*nq = 0; |
return; |
return; |
} else if ( !n1 ) *nq = 0; |
} else if ( !n1 ) *nq = 0; |
else if ( n1 == n2 ) *nq = (GQ)ONEGZ; |
else if ( n1 == n2 ) *nq = (GQ)ONEGZ; |
else { |
else { |
if ( NID(n1) == N_GZ ) MPZTOMPQ(BDY((GZ)n1),q1); |
if ( NID(n1) == N_GZ ) MPZTOMPQ(BDY((GZ)n1),q1); |
else q1[0] = BDY(n1)[0]; |
else q1[0] = BDY(n1)[0]; |
if ( NID(n2) == N_GZ ) MPZTOMPQ(BDY((GZ)n2),q2); |
if ( NID(n2) == N_GZ ) MPZTOMPQ(BDY((GZ)n2),q2); |
else q2[0] = BDY(n2)[0]; |
else q2[0] = BDY(n2)[0]; |
mpq_init(t); mpq_div(t,q1,q2); *nq = mpqtogzq(t); |
mpq_init(t); mpq_div(t,q1,q2); *nq = mpqtogzq(t); |
} |
} |
} |
} |
|
|
void chsgngq(GQ n,GQ *nr) |
void chsgngq(GQ n,GQ *nr) |
{ |
{ |
mpq_t t; |
mpq_t t; |
|
|
if ( !n ) *nr = 0; |
if ( !n ) *nr = 0; |
else if ( NID(n) == N_GZ ) chsgngz((GZ)n,(GZ *)nr); |
else if ( NID(n) == N_GZ ) chsgngz((GZ)n,(GZ *)nr); |
else { |
else { |
mpq_init(t); mpq_neg(t,BDY(n)); MPQTOGQ(t,*nr); |
mpq_init(t); mpq_neg(t,BDY(n)); MPQTOGQ(t,*nr); |
} |
} |
} |
} |
|
|
void pwrgq(GQ n1,Q n,GQ *nr) |
void pwrgq(GQ n1,Q n,GQ *nr) |
{ |
{ |
int e; |
int e; |
mpz_t nm,dn; |
mpz_t nm,dn; |
mpq_t t; |
mpq_t t; |
|
|
if ( !n || UNIGZ((GZ)n1) || UNIGQ(n1) ) *nr = (GQ)ONEGZ; |
if ( !n || UNIGZ((GZ)n1) || UNIGQ(n1) ) *nr = (GQ)ONEGZ; |
else if ( !n1 ) *nr = 0; |
else if ( !n1 ) *nr = 0; |
else if ( !INT(n) ) { |
else if ( !INT(n) ) { |
error("can't calculate fractional power."); *nr = 0; |
error("can't calculate fractional power."); *nr = 0; |
} else if ( PL(NM(n)) > 1 ) { |
} else if ( PL(NM(n)) > 1 ) { |
error("exponent too big."); *nr = 0; |
error("exponent too big."); *nr = 0; |
} else { |
} else { |
e = QTOS(n); |
e = QTOS(n); |
if ( e < 0 ) { |
if ( e < 0 ) { |
e = -e; |
e = -e; |
if ( NID(n1)==N_GZ ) { |
if ( NID(n1)==N_GZ ) { |
nm[0] = ONEMPZ[0]; |
nm[0] = ONEMPZ[0]; |
dn[0] = BDY((GZ)n1)[0]; |
dn[0] = BDY((GZ)n1)[0]; |
} else { |
} else { |
nm[0] = mpq_denref(BDY(n1))[0]; dn[0] = mpq_numref(BDY(n1))[0]; |
nm[0] = mpq_denref(BDY(n1))[0]; dn[0] = mpq_numref(BDY(n1))[0]; |
} |
} |
} else { |
} else { |
if ( NID(n1)==N_GZ ) { |
if ( NID(n1)==N_GZ ) { |
nm[0] = BDY((GZ)n1)[0]; dn[0] = ONEMPZ[0]; |
nm[0] = BDY((GZ)n1)[0]; dn[0] = ONEMPZ[0]; |
} else { |
} else { |
nm[0] = mpq_numref(BDY(n1))[0]; dn[0] = mpq_denref(BDY(n1))[0]; |
nm[0] = mpq_numref(BDY(n1))[0]; dn[0] = mpq_denref(BDY(n1))[0]; |
} |
} |
} |
} |
mpq_init(t); |
mpq_init(t); |
mpz_pow_ui(mpq_numref(t),nm,e); mpz_pow_ui(mpq_denref(t),dn,e); |
mpz_pow_ui(mpq_numref(t),nm,e); mpz_pow_ui(mpq_denref(t),dn,e); |
*nr = mpqtogzq(t); |
*nr = mpqtogzq(t); |
} |
} |
} |
} |
|
|
int cmpgq(GQ n1,GQ n2) |
int cmpgq(GQ n1,GQ n2) |
{ |
{ |
mpq_t q1,q2; |
mpq_t q1,q2; |
int sgn; |
int sgn; |
|
|
if ( !n1 ) |
if ( !n1 ) |
if ( !n2 ) return 0; |
if ( !n2 ) return 0; |
else return (NID(n2)==N_GZ) ? -mpz_sgn(BDY((GZ)n2)) : -mpq_sgn(BDY(n2)); |
else return (NID(n2)==N_GZ) ? -mpz_sgn(BDY((GZ)n2)) : -mpq_sgn(BDY(n2)); |
if ( !n2 ) return (NID(n1)==N_GZ) ? mpz_sgn(BDY((GZ)n1)) : mpq_sgn(BDY(n1)); |
if ( !n2 ) return (NID(n1)==N_GZ) ? mpz_sgn(BDY((GZ)n1)) : mpq_sgn(BDY(n1)); |
else if ( (sgn = mpq_sgn(BDY(n1))) != mpq_sgn(BDY(n2)) ) return sgn; |
else if ( (sgn = mpq_sgn(BDY(n1))) != mpq_sgn(BDY(n2)) ) return sgn; |
else { |
else { |
if ( NID(n1) == N_GZ ) MPZTOMPQ(BDY((GZ)n1),q1); |
if ( NID(n1) == N_GZ ) MPZTOMPQ(BDY((GZ)n1),q1); |
else q1[0] = BDY(n1)[0]; |
else q1[0] = BDY(n1)[0]; |
if ( NID(n2) == N_GZ ) MPZTOMPQ(BDY((GZ)n2),q2); |
if ( NID(n2) == N_GZ ) MPZTOMPQ(BDY((GZ)n2),q2); |
else q2[0] = BDY(n2)[0]; |
else q2[0] = BDY(n2)[0]; |
sgn = mpq_cmp(q1,q2); |
sgn = mpq_cmp(q1,q2); |
if ( sgn > 0 ) return 1; |
if ( sgn > 0 ) return 1; |
else if ( sgn < 0 ) return -1; |
else if ( sgn < 0 ) return -1; |
else return 0; |
else return 0; |
} |
} |
} |
} |
|
|
void mkgwc(int k,int l,GZ *t) |
void mkgwc(int k,int l,GZ *t) |
{ |
{ |
mpz_t a,b,q,nm,z,u; |
mpz_t a,b,q,nm,z,u; |
int i,n; |
int i,n; |
|
|
n = MIN(k,l); |
n = MIN(k,l); |
mpz_init_set_ui(z,1); |
mpz_init_set_ui(z,1); |
mpz_init(u); mpz_set(u,z); MPZTOGZ(u,t[0]); |
mpz_init(u); mpz_set(u,z); MPZTOGZ(u,t[0]); |
mpz_init(a); mpz_init(b); mpz_init(nm); |
mpz_init(a); mpz_init(b); mpz_init(nm); |
for ( i = 1; i <= n; i++ ) { |
for ( i = 1; i <= n; i++ ) { |
mpz_set_ui(a,k-i+1); mpz_set_ui(b,l-i+1); mpz_mul(nm,a,b); |
mpz_set_ui(a,k-i+1); mpz_set_ui(b,l-i+1); mpz_mul(nm,a,b); |
mpz_mul(z,BDY(t[i-1]),nm); mpz_fdiv_q_ui(z,z,i); |
mpz_mul(z,BDY(t[i-1]),nm); mpz_fdiv_q_ui(z,z,i); |
mpz_init(u); mpz_set(u,z); MPZTOGZ(u,t[i]); |
mpz_init(u); mpz_set(u,z); MPZTOGZ(u,t[i]); |
} |
} |
} |
} |
|
|
void gz_lgp(P p,GZ *g,GZ *l); |
void gz_lgp(P p,GZ *g,GZ *l); |
|
|
void gz_ptozp(P p,int sgn,GQ *c,P *pr) |
void gz_ptozp(P p,int sgn,GQ *c,P *pr) |
{ |
{ |
GZ nm,dn; |
GZ nm,dn; |
|
|
if ( !p ) { |
if ( !p ) { |
*c = 0; *pr = 0; |
*c = 0; *pr = 0; |
} else { |
} else { |
gz_lgp(p,&nm,&dn); |
gz_lgp(p,&nm,&dn); |
divgz(nm,dn,(GZ *)c); |
divgz(nm,dn,(GZ *)c); |
divsp(CO,p,(P)*c,pr); |
divsp(CO,p,(P)*c,pr); |
} |
} |
} |
} |
|
|
void gz_lgp(P p,GZ *g,GZ *l) |
void gz_lgp(P p,GZ *g,GZ *l) |
{ |
{ |
DCP dc; |
DCP dc; |
GZ g1,g2,l1,l2,l3,l4; |
GZ g1,g2,l1,l2,l3,l4; |
|
|
if ( NUM(p) ) { |
if ( NUM(p) ) { |
if ( NID((GZ)p)==N_GZ ) { |
if ( NID((GZ)p)==N_GZ ) { |
MPZTOGZ(BDY((GZ)p),*g); |
MPZTOGZ(BDY((GZ)p),*g); |
*l = ONEGZ; |
*l = ONEGZ; |
} else { |
} else { |
MPZTOGZ(mpq_numref(BDY((GQ)p)),*g); |
MPZTOGZ(mpq_numref(BDY((GQ)p)),*g); |
MPZTOGZ(mpq_denref(BDY((GQ)p)),*l); |
MPZTOGZ(mpq_denref(BDY((GQ)p)),*l); |
} |
} |
} else { |
} else { |
dc = DC(p); gz_lgp(COEF(dc),g,l); |
dc = DC(p); gz_lgp(COEF(dc),g,l); |
for ( dc = NEXT(dc); dc; dc = NEXT(dc) ) { |
for ( dc = NEXT(dc); dc; dc = NEXT(dc) ) { |
gz_lgp(COEF(dc),&g1,&l1); gcdgz(*g,g1,&g2); *g = g2; |
gz_lgp(COEF(dc),&g1,&l1); gcdgz(*g,g1,&g2); *g = g2; |
gcdgz(*l,l1,&l2); mulgz(*l,l1,&l3); divgz(l3,l2,l); |
gcdgz(*l,l1,&l2); mulgz(*l,l1,&l3); divgz(l3,l2,l); |
} |
} |
} |
} |
} |
} |
|
|
void gz_qltozl(GQ *w,int n,GZ *dvr) |
void gz_qltozl(GQ *w,int n,GZ *dvr) |
{ |
{ |
GZ nm,dn; |
GZ nm,dn; |
GZ g,g1,l1,l2,l3; |
GZ g,g1,l1,l2,l3; |
GQ c; |
GQ c; |
int i; |
int i; |
struct oVECT v; |
struct oVECT v; |
|
|
for ( i = 0; i < n; i++ ) |
for ( i = 0; i < n; i++ ) |
if ( w[i] && NID(w[i])==N_GQ ) |
if ( w[i] && NID(w[i])==N_GQ ) |
break; |
break; |
if ( i == n ) { |
if ( i == n ) { |
v.id = O_VECT; v.len = n; v.body = (pointer *)w; |
v.id = O_VECT; v.len = n; v.body = (pointer *)w; |
gcdvgz(&v,dvr); return; |
gcdvgz(&v,dvr); return; |
} |
} |
for ( i = 0; !w[i]; i++ ); |
for ( i = 0; !w[i]; i++ ); |
c = w[i]; |
c = w[i]; |
if ( NID(c)==N_GQ ) { |
if ( NID(c)==N_GQ ) { |
MPZTOGZ(mpq_numref(BDY(c)),nm); MPZTOGZ(mpq_denref(BDY(c)),dn); |
MPZTOGZ(mpq_numref(BDY(c)),nm); MPZTOGZ(mpq_denref(BDY(c)),dn); |
} else { |
} else { |
MPZTOGZ(BDY((GZ)c),nm); dn = ONEGZ; |
MPZTOGZ(BDY((GZ)c),nm); dn = ONEGZ; |
} |
} |
for ( i++; i < n; i++ ) { |
for ( i++; i < n; i++ ) { |
c = w[i]; |
c = w[i]; |
if ( !c ) continue; |
if ( !c ) continue; |
if ( NID(c)==N_GQ ) { |
if ( NID(c)==N_GQ ) { |
MPZTOGZ(mpq_numref(BDY(c)),g1); MPZTOGZ(mpq_denref(BDY(c)),l1); |
MPZTOGZ(mpq_numref(BDY(c)),g1); MPZTOGZ(mpq_denref(BDY(c)),l1); |
} else { |
} else { |
MPZTOGZ(BDY((GZ)c),g1); l1 = ONEGZ; |
MPZTOGZ(BDY((GZ)c),g1); l1 = ONEGZ; |
} |
} |
gcdgz(nm,g1,&g); nm = g; |
gcdgz(nm,g1,&g); nm = g; |
gcdgz(dn,l1,&l2); mulgz(dn,l1,&l3); divgz(l3,l2,&dn); |
gcdgz(dn,l1,&l2); mulgz(dn,l1,&l3); divgz(l3,l2,&dn); |
} |
} |
divgz(nm,dn,dvr); |
divgz(nm,dn,dvr); |
} |
} |
|
|
int gz_bits(GQ q) |
int gz_bits(GQ q) |
{ |
{ |
if ( !q ) return 0; |
if ( !q ) return 0; |
else if ( NID(q)==N_Q ) |
else if ( NID(q)==N_Q ) |
return n_bits(NM((Q)q))+(INT((Q)q)?0:n_bits(DN((Q)q))); |
return n_bits(NM((Q)q))+(INT((Q)q)?0:n_bits(DN((Q)q))); |
else if ( NID(q)==N_GZ ) return mpz_sizeinbase(BDY((GZ)q),2); |
else if ( NID(q)==N_GZ ) return mpz_sizeinbase(BDY((GZ)q),2); |
else |
else |
return mpz_sizeinbase(mpq_numref(BDY(q)),2) |
return mpz_sizeinbase(mpq_numref(BDY(q)),2) |
+ mpz_sizeinbase(mpq_denref(BDY(q)),2); |
+ mpz_sizeinbase(mpq_denref(BDY(q)),2); |
} |
} |
|
|
int gzp_mag(P p) |
int gzp_mag(P p) |
{ |
{ |
int s; |
int s; |
DCP dc; |
DCP dc; |
|
|
if ( !p ) return 0; |
if ( !p ) return 0; |
else if ( OID(p) == O_N ) return gz_bits((GQ)p); |
else if ( OID(p) == O_N ) return gz_bits((GQ)p); |
else { |
else { |
for ( dc = DC(p), s = 0; dc; dc = NEXT(dc) ) s += gzp_mag(COEF(dc)); |
for ( dc = DC(p), s = 0; dc; dc = NEXT(dc) ) s += gzp_mag(COEF(dc)); |
return s; |
return s; |
} |
} |
} |
} |
|
|
void makesubstgz(VL v,NODE *s) |
void makesubstgz(VL v,NODE *s) |
{ |
{ |
NODE r,r0; |
NODE r,r0; |
GZ q; |
GZ q; |
unsigned int n; |
unsigned int n; |
|
|
for ( r0 = 0; v; v = NEXT(v) ) { |
for ( r0 = 0; v; v = NEXT(v) ) { |
NEXTNODE(r0,r); BDY(r) = (pointer)v->v; |
NEXTNODE(r0,r); BDY(r) = (pointer)v->v; |
#if defined(_PA_RISC1_1) |
#if defined(_PA_RISC1_1) |
n = mrand48()&BMASK; q = utogz(n); |
n = mrand48()&BMASK; q = utogz(n); |
#else |
#else |
n = random(); q = utogz(n); |
n = random(); q = utogz(n); |
#endif |
#endif |
NEXTNODE(r0,r); BDY(r) = (pointer)q; |
NEXTNODE(r0,r); BDY(r) = (pointer)q; |
} |
} |
if ( r0 ) NEXT(r) = 0; |
if ( r0 ) NEXT(r) = 0; |
*s = r0; |
*s = r0; |
} |
} |
|
|
unsigned int remgq(GQ a,unsigned int mod) |
unsigned int remgq(GQ a,unsigned int mod) |
{ |
{ |
unsigned int c,nm,dn; |
unsigned int c,nm,dn; |
mpz_t r; |
mpz_t r; |
|
|
if ( !a ) return 0; |
if ( !a ) return 0; |
else if ( NID(a)==N_GZ ) { |
else if ( NID(a)==N_GZ ) { |
mpz_init(r); |
mpz_init(r); |
c = mpz_fdiv_r_ui(r,BDY((GZ)a),mod); |
c = mpz_fdiv_r_ui(r,BDY((GZ)a),mod); |
} else { |
} else { |
mpz_init(r); |
mpz_init(r); |
nm = mpz_fdiv_r_ui(r,mpq_numref(BDY(a)),mod); |
nm = mpz_fdiv_r_ui(r,mpq_numref(BDY(a)),mod); |
dn = mpz_fdiv_r_ui(r,mpq_denref(BDY(a)),mod); |
dn = mpz_fdiv_r_ui(r,mpq_denref(BDY(a)),mod); |
dn = invm(dn,mod); |
dn = invm(dn,mod); |
DMAR(nm,dn,0,mod,c); |
DMAR(nm,dn,0,mod,c); |
} |
} |
return c; |
return c; |
} |
} |
|
|
extern int DP_Print; |
extern int DP_Print; |
Line 752 extern int DP_Print; |
|
Line 753 extern int DP_Print; |
|
|
|
int gz_generic_gauss_elim(MAT mat,MAT *nm,GZ *dn,int **rindp,int **cindp) |
int gz_generic_gauss_elim(MAT mat,MAT *nm,GZ *dn,int **rindp,int **cindp) |
{ |
{ |
int **wmat; |
int **wmat; |
GZ **bmat,**tmat,*bmi,*tmi; |
GZ **bmat,**tmat,*bmi,*tmi; |
GZ q,m1,m2,m3,s,u; |
GZ q,m1,m2,m3,s,u; |
int *wmi,*colstat,*wcolstat,*rind,*cind; |
int *wmi,*colstat,*wcolstat,*rind,*cind; |
int row,col,ind,md,i,j,k,l,t,t1,rank,rank0,inv; |
int row,col,ind,md,i,j,k,l,t,t1,rank,rank0,inv; |
MAT r,crmat; |
MAT r,crmat; |
int ret; |
int ret; |
|
|
bmat = (GZ **)mat->body; |
bmat = (GZ **)mat->body; |
row = mat->row; col = mat->col; |
row = mat->row; col = mat->col; |
wmat = (int **)almat(row,col); |
wmat = (int **)almat(row,col); |
colstat = (int *)MALLOC_ATOMIC(col*sizeof(int)); |
colstat = (int *)MALLOC_ATOMIC(col*sizeof(int)); |
wcolstat = (int *)MALLOC_ATOMIC(col*sizeof(int)); |
wcolstat = (int *)MALLOC_ATOMIC(col*sizeof(int)); |
for ( ind = 0; ; ind++ ) { |
for ( ind = 0; ; ind++ ) { |
if ( DP_Print ) { |
if ( DP_Print ) { |
fprintf(asir_out,"."); fflush(asir_out); |
fprintf(asir_out,"."); fflush(asir_out); |
} |
} |
md = get_lprime(ind); |
md = get_lprime(ind); |
for ( i = 0; i < row; i++ ) |
for ( i = 0; i < row; i++ ) |
for ( j = 0, bmi = bmat[i], wmi = wmat[i]; j < col; j++ ) |
for ( j = 0, bmi = bmat[i], wmi = wmat[i]; j < col; j++ ) |
wmi[j] = remgq((GQ)bmi[j],md); |
wmi[j] = remgq((GQ)bmi[j],md); |
rank = generic_gauss_elim_mod(wmat,row,col,md,wcolstat); |
rank = generic_gauss_elim_mod(wmat,row,col,md,wcolstat); |
if ( !ind ) { |
if ( !ind ) { |
RESET: |
RESET: |
m1 = utogz(md); |
m1 = utogz(md); |
rank0 = rank; |
rank0 = rank; |
bcopy(wcolstat,colstat,col*sizeof(int)); |
bcopy(wcolstat,colstat,col*sizeof(int)); |
MKMAT(crmat,rank,col-rank); |
MKMAT(crmat,rank,col-rank); |
MKMAT(r,rank,col-rank); *nm = r; |
MKMAT(r,rank,col-rank); *nm = r; |
tmat = (GZ **)crmat->body; |
tmat = (GZ **)crmat->body; |
for ( i = 0; i < rank; i++ ) |
for ( i = 0; i < rank; i++ ) |
for ( j = k = 0, tmi = tmat[i], wmi = wmat[i]; j < col; j++ ) |
for ( j = k = 0, tmi = tmat[i], wmi = wmat[i]; j < col; j++ ) |
if ( !colstat[j] ) tmi[k++] = utogz(wmi[j]); |
if ( !colstat[j] ) tmi[k++] = utogz(wmi[j]); |
} else { |
} else { |
if ( rank < rank0 ) { |
if ( rank < rank0 ) { |
if ( DP_Print ) { |
if ( DP_Print ) { |
fprintf(asir_out,"lower rank matrix; continuing...\n"); |
fprintf(asir_out,"lower rank matrix; continuing...\n"); |
fflush(asir_out); |
fflush(asir_out); |
} |
} |
continue; |
continue; |
} else if ( rank > rank0 ) { |
} else if ( rank > rank0 ) { |
if ( DP_Print ) { |
if ( DP_Print ) { |
fprintf(asir_out,"higher rank matrix; resetting...\n"); |
fprintf(asir_out,"higher rank matrix; resetting...\n"); |
fflush(asir_out); |
fflush(asir_out); |
} |
} |
goto RESET; |
goto RESET; |
} else { |
} else { |
for ( j = 0; (j<col) && (colstat[j]==wcolstat[j]); j++ ); |
for ( j = 0; (j<col) && (colstat[j]==wcolstat[j]); j++ ); |
if ( j < col ) { |
if ( j < col ) { |
if ( DP_Print ) { |
if ( DP_Print ) { |
fprintf(asir_out,"inconsitent colstat; resetting...\n"); |
fprintf(asir_out,"inconsitent colstat; resetting...\n"); |
fflush(asir_out); |
fflush(asir_out); |
} |
} |
goto RESET; |
goto RESET; |
} |
} |
} |
} |
|
|
inv = invm(remgq((GQ)m1,md),md); |
inv = invm(remgq((GQ)m1,md),md); |
m2 = utogz(md); mulgz(m1,m2,&m3); |
m2 = utogz(md); mulgz(m1,m2,&m3); |
for ( i = 0; i < rank; i++ ) |
for ( i = 0; i < rank; i++ ) |
for ( j = k = 0, tmi = tmat[i], wmi = wmat[i]; j < col; j++ ) |
for ( j = k = 0, tmi = tmat[i], wmi = wmat[i]; j < col; j++ ) |
if ( !colstat[j] ) { |
if ( !colstat[j] ) { |
if ( tmi[k] ) { |
if ( tmi[k] ) { |
/* f3 = f1+m1*(m1 mod m2)^(-1)*(f2 - f1 mod m2) */ |
/* f3 = f1+m1*(m1 mod m2)^(-1)*(f2 - f1 mod m2) */ |
t = remgq((GQ)tmi[k],md); |
t = remgq((GQ)tmi[k],md); |
if ( wmi[j] >= t ) t = wmi[j]-t; |
if ( wmi[j] >= t ) t = wmi[j]-t; |
else t = md-(t-wmi[j]); |
else t = md-(t-wmi[j]); |
DMAR(t,inv,0,md,t1) |
DMAR(t,inv,0,md,t1) |
u = utogz(t1); mulgz(m1,u,&s); |
u = utogz(t1); mulgz(m1,u,&s); |
addgz(tmi[k],s,&u); tmi[k] = u; |
addgz(tmi[k],s,&u); tmi[k] = u; |
} else if ( wmi[j] ) { |
} else if ( wmi[j] ) { |
/* f3 = m1*(m1 mod m2)^(-1)*f2 */ |
/* f3 = m1*(m1 mod m2)^(-1)*f2 */ |
DMAR(wmi[j],inv,0,md,t) |
DMAR(wmi[j],inv,0,md,t) |
u = utogz(t); mulgz(m1,u,&s); tmi[k] = s; |
u = utogz(t); mulgz(m1,u,&s); tmi[k] = s; |
} |
} |
k++; |
k++; |
} |
} |
m1 = m3; |
m1 = m3; |
if ( ind % GZ_F4_INTRAT_PERIOD ) |
if ( ind % GZ_F4_INTRAT_PERIOD ) |
ret = 0; |
ret = 0; |
else |
else |
ret = gz_intmtoratm(crmat,m1,*nm,dn); |
ret = gz_intmtoratm(crmat,m1,*nm,dn); |
if ( ret ) { |
if ( ret ) { |
*rindp = rind = (int *)MALLOC_ATOMIC(rank*sizeof(int)); |
*rindp = rind = (int *)MALLOC_ATOMIC(rank*sizeof(int)); |
*cindp = cind = (int *)MALLOC_ATOMIC((col-rank)*sizeof(int)); |
*cindp = cind = (int *)MALLOC_ATOMIC((col-rank)*sizeof(int)); |
for ( j = k = l = 0; j < col; j++ ) |
for ( j = k = l = 0; j < col; j++ ) |
if ( colstat[j] ) rind[k++] = j; |
if ( colstat[j] ) rind[k++] = j; |
else cind[l++] = j; |
else cind[l++] = j; |
if ( gz_gensolve_check(mat,*nm,*dn,rind,cind) ) |
if ( gz_gensolve_check(mat,*nm,*dn,rind,cind) ) |
return rank; |
return rank; |
} |
} |
} |
} |
} |
} |
} |
} |
|
|
int gz_generic_gauss_elim2(MAT mat,MAT *nm,GZ *dn,int **rindp,int **cindp) |
int gz_generic_gauss_elim2(MAT mat,MAT *nm,GZ *dn,int **rindp,int **cindp) |
{ |
{ |
|
|
MAT full; |
MAT full; |
GZ **bmat,**b; |
GZ **bmat,**b; |
GZ *bmi; |
GZ *bmi; |
GZ dn0; |
GZ dn0; |
int row,col,md,i,j,rank,ret; |
int row,col,md,i,j,rank,ret; |
int **wmat; |
int **wmat; |
int *wmi; |
int *wmi; |
int *colstat,*rowstat; |
int *colstat,*rowstat; |
|
|
bmat = (GZ **)mat->body; |
bmat = (GZ **)mat->body; |
row = mat->row; col = mat->col; |
row = mat->row; col = mat->col; |
wmat = (int **)almat(row,col); |
wmat = (int **)almat(row,col); |
colstat = (int *)MALLOC_ATOMIC(col*sizeof(int)); |
colstat = (int *)MALLOC_ATOMIC(col*sizeof(int)); |
rowstat = (int *)MALLOC_ATOMIC(row*sizeof(int)); |
rowstat = (int *)MALLOC_ATOMIC(row*sizeof(int)); |
/* XXX */ |
/* XXX */ |
md = get_lprime(0); |
md = get_lprime(0); |
for ( i = 0; i < row; i++ ) |
for ( i = 0; i < row; i++ ) |
for ( j = 0, bmi = bmat[i], wmi = wmat[i]; j < col; j++ ) |
for ( j = 0, bmi = bmat[i], wmi = wmat[i]; j < col; j++ ) |
wmi[j] = remgq((GQ)bmi[j],md); |
wmi[j] = remgq((GQ)bmi[j],md); |
rank = generic_gauss_elim_mod2(wmat,row,col,md,colstat,rowstat); |
rank = generic_gauss_elim_mod2(wmat,row,col,md,colstat,rowstat); |
b = (GZ **)MALLOC(rank*sizeof(GZ)); |
b = (GZ **)MALLOC(rank*sizeof(GZ)); |
for ( i = 0; i < rank; i++ ) b[i] = bmat[rowstat[i]]; |
for ( i = 0; i < rank; i++ ) b[i] = bmat[rowstat[i]]; |
NEWMAT(full); full->row = rank; full->col = col; full->body = (pointer **)b; |
NEWMAT(full); full->row = rank; full->col = col; full->body = (pointer **)b; |
ret = gz_generic_gauss_elim_full(full,nm,dn,rindp,cindp); |
ret = gz_generic_gauss_elim_full(full,nm,dn,rindp,cindp); |
if ( !ret ) { |
if ( !ret ) { |
rank = gz_generic_gauss_elim(mat,nm,&dn0,rindp,cindp); |
rank = gz_generic_gauss_elim(mat,nm,&dn0,rindp,cindp); |
for ( i = 0; i < rank; i++ ) dn[i] = dn0; |
for ( i = 0; i < rank; i++ ) dn[i] = dn0; |
} |
} |
return rank; |
return rank; |
} |
} |
|
|
int gz_generic_gauss_elim_full(MAT mat,MAT *nm,GZ *dn,int **rindp,int **cindp) |
int gz_generic_gauss_elim_full(MAT mat,MAT *nm,GZ *dn,int **rindp,int **cindp) |
{ |
{ |
int **wmat; |
int **wmat; |
int *stat; |
int *stat; |
GZ **bmat,**tmat,*bmi,*tmi,*ri; |
GZ **bmat,**tmat,*bmi,*tmi,*ri; |
GZ q,m1,m2,m3,s,u; |
GZ q,m1,m2,m3,s,u; |
int *wmi,*colstat,*wcolstat,*rind,*cind; |
int *wmi,*colstat,*wcolstat,*rind,*cind; |
int row,col,ind,md,i,j,k,l,t,t1,rank,rank0,inv,h; |
int row,col,ind,md,i,j,k,l,t,t1,rank,rank0,inv,h; |
MAT r,crmat; |
MAT r,crmat; |
int ret,initialized,done; |
int ret,initialized,done; |
|
|
initialized = 0; |
initialized = 0; |
bmat = (GZ **)mat->body; |
bmat = (GZ **)mat->body; |
row = mat->row; col = mat->col; |
row = mat->row; col = mat->col; |
wmat = (int **)almat(row,col); |
wmat = (int **)almat(row,col); |
stat = (int *)MALLOC_ATOMIC(row*sizeof(int)); |
stat = (int *)MALLOC_ATOMIC(row*sizeof(int)); |
for ( i = 0; i < row; i++ ) stat[i] = 0; |
for ( i = 0; i < row; i++ ) stat[i] = 0; |
colstat = (int *)MALLOC_ATOMIC(col*sizeof(int)); |
colstat = (int *)MALLOC_ATOMIC(col*sizeof(int)); |
wcolstat = (int *)MALLOC_ATOMIC(col*sizeof(int)); |
wcolstat = (int *)MALLOC_ATOMIC(col*sizeof(int)); |
for ( ind = 0; ; ind++ ) { |
for ( ind = 0; ; ind++ ) { |
if ( DP_Print ) { |
if ( DP_Print ) { |
fprintf(asir_out,"."); fflush(asir_out); |
fprintf(asir_out,"."); fflush(asir_out); |
} |
} |
md = get_lprime(ind); |
md = get_lprime(ind); |
for ( i = 0; i < row; i++ ) |
for ( i = 0; i < row; i++ ) |
for ( j = 0, bmi = bmat[i], wmi = wmat[i]; j < col; j++ ) |
for ( j = 0, bmi = bmat[i], wmi = wmat[i]; j < col; j++ ) |
wmi[j] = remgq((GQ)bmi[j],md); |
wmi[j] = remgq((GQ)bmi[j],md); |
rank = generic_gauss_elim_mod(wmat,row,col,md,wcolstat); |
rank = generic_gauss_elim_mod(wmat,row,col,md,wcolstat); |
if ( rank < row ) continue; |
if ( rank < row ) continue; |
if ( !initialized ) { |
if ( !initialized ) { |
m1 = utogz(md); |
m1 = utogz(md); |
bcopy(wcolstat,colstat,col*sizeof(int)); |
bcopy(wcolstat,colstat,col*sizeof(int)); |
MKMAT(crmat,row,col-row); |
MKMAT(crmat,row,col-row); |
MKMAT(r,row,col-row); *nm = r; |
MKMAT(r,row,col-row); *nm = r; |
tmat = (GZ **)crmat->body; |
tmat = (GZ **)crmat->body; |
for ( i = 0; i < row; i++ ) |
for ( i = 0; i < row; i++ ) |
for ( j = k = 0, tmi = tmat[i], wmi = wmat[i]; j < col; j++ ) |
for ( j = k = 0, tmi = tmat[i], wmi = wmat[i]; j < col; j++ ) |
if ( !colstat[j] ) tmi[k++] = utogz(wmi[j]); |
if ( !colstat[j] ) tmi[k++] = utogz(wmi[j]); |
initialized = 1; |
initialized = 1; |
} else { |
} else { |
for ( j = 0; (j<col) && (colstat[j]==wcolstat[j]); j++ ); |
for ( j = 0; (j<col) && (colstat[j]==wcolstat[j]); j++ ); |
if ( j < col ) continue; |
if ( j < col ) continue; |
|
|
inv = invm(remgq((GQ)m1,md),md); |
inv = invm(remgq((GQ)m1,md),md); |
m2 = utogz(md); mulgz(m1,m2,&m3); |
m2 = utogz(md); mulgz(m1,m2,&m3); |
for ( i = 0; i < row; i++ ) |
for ( i = 0; i < row; i++ ) |
switch ( stat[i] ) { |
switch ( stat[i] ) { |
case 1: |
case 1: |
/* consistency check */ |
/* consistency check */ |
ri = (GZ *)BDY(r)[i]; wmi = wmat[i]; |
ri = (GZ *)BDY(r)[i]; wmi = wmat[i]; |
for ( j = 0; j < col; j++ ) if ( colstat[j] ) break; |
for ( j = 0; j < col; j++ ) if ( colstat[j] ) break; |
h = md-remgq((GQ)dn[i],md); |
h = md-remgq((GQ)dn[i],md); |
for ( j++, k = 0; j < col; j++ ) |
for ( j++, k = 0; j < col; j++ ) |
if ( !colstat[j] ) { |
if ( !colstat[j] ) { |
t = remgq((GQ)ri[k],md); |
t = remgq((GQ)ri[k],md); |
DMAR(wmi[i],h,t,md,t1); |
DMAR(wmi[i],h,t,md,t1); |
if ( t1 ) break; |
if ( t1 ) break; |
} |
} |
if ( j == col ) { stat[i]++; break; } |
if ( j == col ) { stat[i]++; break; } |
else { |
else { |
/* fall to the case 0 */ |
/* fall to the case 0 */ |
stat[i] = 0; |
stat[i] = 0; |
} |
} |
case 0: |
case 0: |
tmi = tmat[i]; wmi = wmat[i]; |
tmi = tmat[i]; wmi = wmat[i]; |
for ( j = k = 0; j < col; j++ ) |
for ( j = k = 0; j < col; j++ ) |
if ( !colstat[j] ) { |
if ( !colstat[j] ) { |
if ( tmi[k] ) { |
if ( tmi[k] ) { |
/* f3 = f1+m1*(m1 mod m2)^(-1)*(f2 - f1 mod m2) */ |
/* f3 = f1+m1*(m1 mod m2)^(-1)*(f2 - f1 mod m2) */ |
t = remgq((GQ)tmi[k],md); |
t = remgq((GQ)tmi[k],md); |
if ( wmi[j] >= t ) t = wmi[j]-t; |
if ( wmi[j] >= t ) t = wmi[j]-t; |
else t = md-(t-wmi[j]); |
else t = md-(t-wmi[j]); |
DMAR(t,inv,0,md,t1) |
DMAR(t,inv,0,md,t1) |
u = utogz(t1); mulgz(m1,u,&s); |
u = utogz(t1); mulgz(m1,u,&s); |
addgz(tmi[k],s,&u); tmi[k] = u; |
addgz(tmi[k],s,&u); tmi[k] = u; |
} else if ( wmi[j] ) { |
} else if ( wmi[j] ) { |
/* f3 = m1*(m1 mod m2)^(-1)*f2 */ |
/* f3 = m1*(m1 mod m2)^(-1)*f2 */ |
DMAR(wmi[j],inv,0,md,t) |
DMAR(wmi[j],inv,0,md,t) |
u = utogz(t); mulgz(m1,u,&s); tmi[k] = s; |
u = utogz(t); mulgz(m1,u,&s); tmi[k] = s; |
} |
} |
k++; |
k++; |
} |
} |
break; |
break; |
case 2: default: |
case 2: default: |
break; |
break; |
} |
} |
m1 = m3; |
m1 = m3; |
if ( ind % 4 ) |
if ( ind % 4 ) |
ret = 0; |
ret = 0; |
else |
else |
ret = gz_intmtoratm2(crmat,m1,*nm,dn,stat); |
ret = gz_intmtoratm2(crmat,m1,*nm,dn,stat); |
if ( ret ) { |
if ( ret ) { |
*rindp = rind = (int *)MALLOC_ATOMIC(row*sizeof(int)); |
*rindp = rind = (int *)MALLOC_ATOMIC(row*sizeof(int)); |
*cindp = cind = (int *)MALLOC_ATOMIC((col-row)*sizeof(int)); |
*cindp = cind = (int *)MALLOC_ATOMIC((col-row)*sizeof(int)); |
for ( j = k = l = 0; j < col; j++ ) |
for ( j = k = l = 0; j < col; j++ ) |
if ( colstat[j] ) rind[k++] = j; |
if ( colstat[j] ) rind[k++] = j; |
else cind[l++] = j; |
else cind[l++] = j; |
return gz_gensolve_check2(mat,*nm,dn,rind,cind); |
return gz_gensolve_check2(mat,*nm,dn,rind,cind); |
} |
} |
} |
} |
} |
} |
} |
} |
|
|
int gz_generic_gauss_elim_direct(MAT mat,MAT *nm,GZ *dn,int **rindp,int **cindp){ |
int gz_generic_gauss_elim_direct(MAT mat,MAT *nm,GZ *dn,int **rindp,int **cindp){ |
GZ **bmat,*s; |
GZ **bmat,*s; |
GZ u,v,w,x,d,t,y; |
GZ u,v,w,x,d,t,y; |
int row,col,i,j,k,l,m,rank; |
int row,col,i,j,k,l,m,rank; |
int *colstat,*colpos,*cind; |
int *colstat,*colpos,*cind; |
MAT r,in; |
MAT r,in; |
|
|
row = mat->row; col = mat->col; |
row = mat->row; col = mat->col; |
MKMAT(in,row,col); |
MKMAT(in,row,col); |
for ( i = 0; i < row; i++ ) |
for ( i = 0; i < row; i++ ) |
for ( j = 0; j < col; j++ ) in->body[i][j] = mat->body[i][j]; |
for ( j = 0; j < col; j++ ) in->body[i][j] = mat->body[i][j]; |
bmat = (GZ **)in->body; |
bmat = (GZ **)in->body; |
colstat = (int *)MALLOC_ATOMIC(col*sizeof(int)); |
colstat = (int *)MALLOC_ATOMIC(col*sizeof(int)); |
*rindp = colpos = (int *)MALLOC_ATOMIC(row*sizeof(int)); |
*rindp = colpos = (int *)MALLOC_ATOMIC(row*sizeof(int)); |
for ( j = 0, rank = 0, d = ONEGZ; j < col; j++ ) { |
for ( j = 0, rank = 0, d = ONEGZ; j < col; j++ ) { |
for ( i = rank; i < row && !bmat[i][j]; i++ ); |
for ( i = rank; i < row && !bmat[i][j]; i++ ); |
if ( i == row ) { colstat[j] = 0; continue; } |
if ( i == row ) { colstat[j] = 0; continue; } |
else { colstat[j] = 1; colpos[rank] = j; } |
else { colstat[j] = 1; colpos[rank] = j; } |
if ( i != rank ) |
if ( i != rank ) |
for ( k = j; k < col; k++ ) { |
for ( k = j; k < col; k++ ) { |
t = bmat[i][k]; bmat[i][k] = bmat[rank][k]; bmat[rank][k] = t; |
t = bmat[i][k]; bmat[i][k] = bmat[rank][k]; bmat[rank][k] = t; |
} |
} |
for ( i = rank+1, v = bmat[rank][j]; i < row; i++ ) |
for ( i = rank+1, v = bmat[rank][j]; i < row; i++ ) |
for ( k = j, u = bmat[i][j]; k < col; k++ ) { |
for ( k = j, u = bmat[i][j]; k < col; k++ ) { |
mulgz(bmat[i][k],v,&w); mulgz(bmat[rank][k],u,&x); |
mulgz(bmat[i][k],v,&w); mulgz(bmat[rank][k],u,&x); |
subgz(w,x,&y); divsgz(y,d,&bmat[i][k]); |
subgz(w,x,&y); divsgz(y,d,&bmat[i][k]); |
} |
} |
d = v; rank++; |
d = v; rank++; |
} |
} |
*dn = d; |
*dn = d; |
s = (GZ *)MALLOC(col*sizeof(GZ)); |
s = (GZ *)MALLOC(col*sizeof(GZ)); |
for ( i = rank-1; i >= 0; i-- ) { |
for ( i = rank-1; i >= 0; i-- ) { |
for ( k = colpos[i]; k < col; k++ ) mulgz(bmat[i][k],d,&s[k]); |
for ( k = colpos[i]; k < col; k++ ) mulgz(bmat[i][k],d,&s[k]); |
for ( m = rank-1; m > i; m-- ) { |
for ( m = rank-1; m > i; m-- ) { |
for ( k = colpos[m], u = bmat[i][k]; k < col; k++ ) { |
for ( k = colpos[m], u = bmat[i][k]; k < col; k++ ) { |
mulgz(bmat[m][k],u,&w); subgz(s[k],w,&x); s[k] = x; |
mulgz(bmat[m][k],u,&w); subgz(s[k],w,&x); s[k] = x; |
} |
} |
} |
} |
for ( k = colpos[i], u = bmat[i][k]; k < col; k++ ) |
for ( k = colpos[i], u = bmat[i][k]; k < col; k++ ) |
divgz(s[k],u,&bmat[i][k]); |
divgz(s[k],u,&bmat[i][k]); |
} |
} |
*cindp = cind = (int *)MALLOC_ATOMIC((col-rank)*sizeof(int)); |
*cindp = cind = (int *)MALLOC_ATOMIC((col-rank)*sizeof(int)); |
MKMAT(r,rank,col-rank); *nm = r; |
MKMAT(r,rank,col-rank); *nm = r; |
for ( j = 0, k = 0; j < col; j++ ) |
for ( j = 0, k = 0; j < col; j++ ) |
if ( !colstat[j] ) { |
if ( !colstat[j] ) { |
cind[k] = j; |
cind[k] = j; |
for ( i = 0; i < rank; i++ ) r->body[i][k] = bmat[i][j]; |
for ( i = 0; i < rank; i++ ) r->body[i][k] = bmat[i][j]; |
k++; |
k++; |
} |
} |
return rank; |
return rank; |
} |
} |
|
|
int gz_intmtoratm(MAT mat,GZ md,MAT nm,GZ *dn) |
int gz_intmtoratm(MAT mat,GZ md,MAT nm,GZ *dn) |
{ |
{ |
GZ t,s,b,dn0,dn1,nm1,q,u,unm,udn,dmy; |
GZ t,s,b,dn0,dn1,nm1,q,u,unm,udn,dmy; |
int i,j,k,l,row,col,sgn,ret; |
int i,j,k,l,row,col,sgn,ret; |
GZ **rmat,**tmat,*tmi,*nmk; |
GZ **rmat,**tmat,*tmi,*nmk; |
|
|
if ( UNIGZ(md) ) |
if ( UNIGZ(md) ) |
return 0; |
return 0; |
row = mat->row; col = mat->col; |
row = mat->row; col = mat->col; |
bshiftgz(md,1,&t); |
bshiftgz(md,1,&t); |
isqrtgz(t,&s); |
isqrtgz(t,&s); |
bshiftgz(s,64,&b); |
bshiftgz(s,64,&b); |
if ( !b ) b = ONEGZ; |
if ( !b ) b = ONEGZ; |
dn0 = ONEGZ; |
dn0 = ONEGZ; |
tmat = (GZ **)mat->body; |
tmat = (GZ **)mat->body; |
rmat = (GZ **)nm->body; |
rmat = (GZ **)nm->body; |
for ( i = 0; i < row; i++ ) |
for ( i = 0; i < row; i++ ) |
for ( j = 0, tmi = tmat[i]; j < col; j++ ) |
for ( j = 0, tmi = tmat[i]; j < col; j++ ) |
if ( tmi[j] ) { |
if ( tmi[j] ) { |
mulgz(tmi[j],dn0,&s); |
mulgz(tmi[j],dn0,&s); |
divqrgz(s,md,&dmy,&u); |
divqrgz(s,md,&dmy,&u); |
ret = gz_inttorat(u,md,b,&sgn,&unm,&udn); |
ret = gz_inttorat(u,md,b,&sgn,&unm,&udn); |
if ( !ret ) return 0; |
if ( !ret ) return 0; |
else { |
else { |
if ( sgn < 0 ) chsgngz(unm,&nm1); |
if ( sgn < 0 ) chsgngz(unm,&nm1); |
else nm1 = unm; |
else nm1 = unm; |
dn1 = udn; |
dn1 = udn; |
if ( !UNIGZ(dn1) ) { |
if ( !UNIGZ(dn1) ) { |
for ( k = 0; k < i; k++ ) |
for ( k = 0; k < i; k++ ) |
for ( l = 0, nmk = rmat[k]; l < col; l++ ) { |
for ( l = 0, nmk = rmat[k]; l < col; l++ ) { |
mulgz(nmk[l],dn1,&q); nmk[l] = q; |
mulgz(nmk[l],dn1,&q); nmk[l] = q; |
} |
} |
for ( l = 0, nmk = rmat[i]; l < j; l++ ) { |
for ( l = 0, nmk = rmat[i]; l < j; l++ ) { |
mulgz(nmk[l],dn1,&q); nmk[l] = q; |
mulgz(nmk[l],dn1,&q); nmk[l] = q; |
} |
} |
} |
} |
rmat[i][j] = nm1; |
rmat[i][j] = nm1; |
mulgz(dn0,dn1,&q); dn0 = q; |
mulgz(dn0,dn1,&q); dn0 = q; |
} |
} |
} |
} |
*dn = dn0; |
*dn = dn0; |
return 1; |
return 1; |
} |
} |
|
|
int gz_intmtoratm2(MAT mat,GZ md,MAT nm,GZ *dn,int *stat) |
int gz_intmtoratm2(MAT mat,GZ md,MAT nm,GZ *dn,int *stat) |
{ |
{ |
int row,col,i,j,ret; |
int row,col,i,j,ret; |
GZ dn0,dn1,t,s,b; |
GZ dn0,dn1,t,s,b; |
GZ *w,*tmi; |
GZ *w,*tmi; |
GZ **tmat; |
GZ **tmat; |
|
|
bshiftgz(md,1,&t); |
bshiftgz(md,1,&t); |
isqrtgz(t,&s); |
isqrtgz(t,&s); |
bshiftgz(s,64,&b); |
bshiftgz(s,64,&b); |
tmat = (GZ **)mat->body; |
tmat = (GZ **)mat->body; |
if ( UNIGZ(md) ) return 0; |
if ( UNIGZ(md) ) return 0; |
row = mat->row; col = mat->col; |
row = mat->row; col = mat->col; |
dn0 = ONEGZ; |
dn0 = ONEGZ; |
for ( i = 0; i < row; i++ ) |
for ( i = 0; i < row; i++ ) |
if ( cmpgz(dn[i],dn0) > 0 ) dn0 = dn[i]; |
if ( cmpgz(dn[i],dn0) > 0 ) dn0 = dn[i]; |
w = (GZ *)MALLOC(col*sizeof(GZ)); |
w = (GZ *)MALLOC(col*sizeof(GZ)); |
for ( i = 0; i < row; i++ ) |
for ( i = 0; i < row; i++ ) |
if ( stat[i] == 0 ) { |
if ( stat[i] == 0 ) { |
for ( j = 0, tmi = tmat[i]; j < col; j++ ) |
for ( j = 0, tmi = tmat[i]; j < col; j++ ) |
mulgz(tmi[j],dn0,&w[j]); |
mulgz(tmi[j],dn0,&w[j]); |
ret = gz_intvtoratv(w,col,md,b,BDY(nm)[i],&dn[i]); |
ret = gz_intvtoratv(w,col,md,b,BDY(nm)[i],&dn[i]); |
if ( ret ) { |
if ( ret ) { |
stat[i] = 1; |
stat[i] = 1; |
mulgz(dn0,dn[i],&t); dn[i] = t; dn0 = t; |
mulgz(dn0,dn[i],&t); dn[i] = t; dn0 = t; |
} |
} |
} |
} |
for ( i = 0; i < row; i++ ) if ( !stat[i] ) break; |
for ( i = 0; i < row; i++ ) if ( !stat[i] ) break; |
if ( i == row ) return 1; |
if ( i == row ) return 1; |
else return 0; |
else return 0; |
} |
} |
|
|
int gz_intvtoratv(GZ *v,int n,GZ md,GZ b,GZ *nm,GZ *dn) |
int gz_intvtoratv(GZ *v,int n,GZ md,GZ b,GZ *nm,GZ *dn) |
{ |
{ |
GZ dn0,dn1,q,s,u,nm1,unm,udn,dmy; |
GZ dn0,dn1,q,s,u,nm1,unm,udn,dmy; |
GZ *nmk; |
GZ *nmk; |
int j,l,col,ret,sgn; |
int j,l,col,ret,sgn; |
|
|
for ( j = 0; j < n; j++ ) nm[j] = 0; |
for ( j = 0; j < n; j++ ) nm[j] = 0; |
dn0 = ONEGZ; |
dn0 = ONEGZ; |
for ( j = 0; j < n; j++ ) { |
for ( j = 0; j < n; j++ ) { |
if ( !v[j] ) continue; |
if ( !v[j] ) continue; |
mulgz(v[j],dn0,&s); |
mulgz(v[j],dn0,&s); |
divqrgz(s,md,&dmy,&u); |
divqrgz(s,md,&dmy,&u); |
ret = gz_inttorat(u,md,b,&sgn,&unm,&udn); |
ret = gz_inttorat(u,md,b,&sgn,&unm,&udn); |
if ( !ret ) return 0; |
if ( !ret ) return 0; |
if ( sgn < 0 ) chsgngz(unm,&nm1); |
if ( sgn < 0 ) chsgngz(unm,&nm1); |
else nm1 = unm; |
else nm1 = unm; |
dn1 = udn; |
dn1 = udn; |
if ( !UNIGZ(dn1) ) |
if ( !UNIGZ(dn1) ) |
for ( l = 0; l < j; l++ ) { |
for ( l = 0; l < j; l++ ) { |
mulgz(nm[l],dn1,&q); nm[l] = q; |
mulgz(nm[l],dn1,&q); nm[l] = q; |
} |
} |
nm[j] = nm1; |
nm[j] = nm1; |
mulgz(dn0,dn1,&q); dn0 = q; |
mulgz(dn0,dn1,&q); dn0 = q; |
} |
} |
*dn = dn0; |
*dn = dn0; |
return 1; |
return 1; |
} |
} |
|
|
/* assuming 0 < c < m */ |
/* assuming 0 < c < m */ |
|
|
int gz_inttorat(GZ c,GZ m,GZ b,int *sgnp,GZ *nmp,GZ *dnp) |
int gz_inttorat(GZ c,GZ m,GZ b,int *sgnp,GZ *nmp,GZ *dnp) |
{ |
{ |
GZ qq,t,u1,v1,r1; |
GZ qq,t,u1,v1,r1; |
GZ q,u2,v2,r2; |
GZ q,u2,v2,r2; |
|
|
u1 = 0; v1 = ONEGZ; u2 = m; v2 = c; |
u1 = 0; v1 = ONEGZ; u2 = m; v2 = c; |
while ( cmpgz(v2,b) >= 0 ) { |
while ( cmpgz(v2,b) >= 0 ) { |
divqrgz(u2,v2,&q,&r2); u2 = v2; v2 = r2; |
divqrgz(u2,v2,&q,&r2); u2 = v2; v2 = r2; |
mulgz(q,v1,&t); subgz(u1,t,&r1); u1 = v1; v1 = r1; |
mulgz(q,v1,&t); subgz(u1,t,&r1); u1 = v1; v1 = r1; |
} |
} |
if ( cmpgz(v1,b) >= 0 ) return 0; |
if ( cmpgz(v1,b) >= 0 ) return 0; |
else { |
else { |
*nmp = v2; |
*nmp = v2; |
if ( mpz_sgn(BDY(v1))<0 ) { |
if ( mpz_sgn(BDY(v1))<0 ) { |
*sgnp = -1; chsgngz(v1,dnp); |
*sgnp = -1; chsgngz(v1,dnp); |
} else { |
} else { |
*sgnp = 1; *dnp = v1; |
*sgnp = 1; *dnp = v1; |
} |
} |
return 1; |
return 1; |
} |
} |
} |
} |
|
|
extern int f4_nocheck; |
extern int f4_nocheck; |
|
|
int gz_gensolve_check(MAT mat,MAT nm,GZ dn,int *rind,int *cind) |
int gz_gensolve_check(MAT mat,MAT nm,GZ dn,int *rind,int *cind) |
{ |
{ |
int row,col,rank,clen,i,j,k,l; |
int row,col,rank,clen,i,j,k,l; |
GZ s,t; |
GZ s,t; |
GZ *w; |
GZ *w; |
GZ *mati,*nmk; |
GZ *mati,*nmk; |
|
|
if ( f4_nocheck ) return 1; |
if ( f4_nocheck ) return 1; |
row = mat->row; col = mat->col; rank = nm->row; clen = nm->col; |
row = mat->row; col = mat->col; rank = nm->row; clen = nm->col; |
w = (GZ *)MALLOC(clen*sizeof(GZ)); |
w = (GZ *)MALLOC(clen*sizeof(GZ)); |
for ( i = 0; i < row; i++ ) { |
for ( i = 0; i < row; i++ ) { |
mati = (GZ *)mat->body[i]; |
mati = (GZ *)mat->body[i]; |
bzero(w,clen*sizeof(GZ)); |
bzero(w,clen*sizeof(GZ)); |
for ( k = 0; k < rank; k++ ) |
for ( k = 0; k < rank; k++ ) |
for ( l = 0, nmk = (GZ *)nm->body[k]; l < clen; l++ ) { |
for ( l = 0, nmk = (GZ *)nm->body[k]; l < clen; l++ ) { |
mulgz(mati[rind[k]],nmk[l],&t); addgz(w[l],t,&s); w[l] = s; |
mulgz(mati[rind[k]],nmk[l],&t); addgz(w[l],t,&s); w[l] = s; |
} |
} |
for ( j = 0; j < clen; j++ ) { |
for ( j = 0; j < clen; j++ ) { |
mulgz(dn,mati[cind[j]],&t); |
mulgz(dn,mati[cind[j]],&t); |
if ( cmpgz(w[j],t) ) break; |
if ( cmpgz(w[j],t) ) break; |
} |
} |
if ( j != clen ) break; |
if ( j != clen ) break; |
} |
} |
if ( i != row ) return 0; |
if ( i != row ) return 0; |
else return 1; |
else return 1; |
} |
} |
|
|
int gz_gensolve_check2(MAT mat,MAT nm,GZ *dn,int *rind,int *cind) |
int gz_gensolve_check2(MAT mat,MAT nm,GZ *dn,int *rind,int *cind) |
{ |
{ |
int row,col,rank,clen,i,j,k,l; |
int row,col,rank,clen,i,j,k,l; |
GZ s,t,u,d; |
GZ s,t,u,d; |
GZ *w,*m; |
GZ *w,*m; |
GZ *mati,*nmk; |
GZ *mati,*nmk; |
|
|
if ( f4_nocheck ) return 1; |
if ( f4_nocheck ) return 1; |
row = mat->row; col = mat->col; rank = nm->row; clen = nm->col; |
row = mat->row; col = mat->col; rank = nm->row; clen = nm->col; |
w = (GZ *)MALLOC(clen*sizeof(GZ)); |
w = (GZ *)MALLOC(clen*sizeof(GZ)); |
m = (GZ *)MALLOC(clen*sizeof(GZ)); |
m = (GZ *)MALLOC(clen*sizeof(GZ)); |
for ( d = dn[0], i = 1; i < rank; i++ ) { |
for ( d = dn[0], i = 1; i < rank; i++ ) { |
lcmgz(d,dn[i],&t); d = t; |
lcmgz(d,dn[i],&t); d = t; |
} |
} |
for ( i = 0; i < rank; i++ ) divsgz(d,dn[i],&m[i]); |
for ( i = 0; i < rank; i++ ) divsgz(d,dn[i],&m[i]); |
for ( i = 0; i < row; i++ ) { |
for ( i = 0; i < row; i++ ) { |
mati = (GZ *)mat->body[i]; |
mati = (GZ *)mat->body[i]; |
bzero(w,clen*sizeof(GZ)); |
bzero(w,clen*sizeof(GZ)); |
for ( k = 0; k < rank; k++ ) { |
for ( k = 0; k < rank; k++ ) { |
mulgz(mati[rind[k]],m[k],&u); |
mulgz(mati[rind[k]],m[k],&u); |
for ( l = 0, nmk = (GZ *)nm->body[k]; l < clen; l++ ) { |
for ( l = 0, nmk = (GZ *)nm->body[k]; l < clen; l++ ) { |
mulgz(u,nmk[l],&t); addgz(w[l],t,&s); w[l] = s; |
mulgz(u,nmk[l],&t); addgz(w[l],t,&s); w[l] = s; |
} |
} |
} |
} |
for ( j = 0; j < clen; j++ ) { |
for ( j = 0; j < clen; j++ ) { |
mulgz(d,mati[cind[j]],&t); |
mulgz(d,mati[cind[j]],&t); |
if ( cmpgz(w[j],t) ) break; |
if ( cmpgz(w[j],t) ) break; |
} |
} |
if ( j != clen ) break; |
if ( j != clen ) break; |
} |
} |
if ( i != row ) return 0; |
if ( i != row ) return 0; |
else return 1; |
else return 1; |
} |
} |
|
|
void isqrtgz(GZ a,GZ *r) |
void isqrtgz(GZ a,GZ *r) |
{ |
{ |
int k; |
int k; |
GZ x,t,x2,xh,quo,rem; |
GZ x,t,x2,xh,quo,rem; |
Q two; |
Q two; |
|
|
if ( !a ) *r = 0; |
if ( !a ) *r = 0; |
else if ( UNIGZ(a) ) *r = ONEGZ; |
else if ( UNIGZ(a) ) *r = ONEGZ; |
else { |
else { |
k = gz_bits((GQ)a); /* a <= 2^k-1 */ |
k = gz_bits((GQ)a); /* a <= 2^k-1 */ |
bshiftgz(ONEGZ,-((k>>1)+(k&1)),&x); /* a <= x^2 */ |
bshiftgz(ONEGZ,-((k>>1)+(k&1)),&x); /* a <= x^2 */ |
STOQ(2,two); |
STOQ(2,two); |
while ( 1 ) { |
while ( 1 ) { |
pwrgz(x,two,&t); |
pwrgz(x,two,&t); |
if ( cmpgz(t,a) <= 0 ) { |
if ( cmpgz(t,a) <= 0 ) { |
*r = x; return; |
*r = x; return; |
} else { |
} else { |
if ( mpz_tstbit(BDY(x),0) ) addgz(x,a,&t); |
if ( mpz_tstbit(BDY(x),0) ) addgz(x,a,&t); |
else t = a; |
else t = a; |
bshiftgz(x,-1,&x2); divqrgz(t,x2,&quo,&rem); |
bshiftgz(x,-1,&x2); divqrgz(t,x2,&quo,&rem); |
bshiftgz(x,1,&xh); addgz(quo,xh,&x); |
bshiftgz(x,1,&xh); addgz(quo,xh,&x); |
} |
} |
} |
} |
} |
} |
} |
} |
|
|
void bshiftgz(GZ a,int n,GZ *r) |
void bshiftgz(GZ a,int n,GZ *r) |
{ |
{ |
mpz_t t; |
mpz_t t; |
|
|
if ( !a ) *r = 0; |
if ( !a ) *r = 0; |
else if ( n == 0 ) *r = a; |
else if ( n == 0 ) *r = a; |
else if ( n < 0 ) { |
else if ( n < 0 ) { |
mpz_init(t); mpz_mul_2exp(t,BDY(a),-n); MPZTOGZ(t,*r); |
mpz_init(t); mpz_mul_2exp(t,BDY(a),-n); MPZTOGZ(t,*r); |
} else { |
} else { |
mpz_init(t); mpz_fdiv_q_2exp(t,BDY(a),n); |
mpz_init(t); mpz_fdiv_q_2exp(t,BDY(a),n); |
if ( !mpz_sgn(t) ) *r = 0; |
if ( !mpz_sgn(t) ) *r = 0; |
else MPZTOGZ(t,*r); |
else MPZTOGZ(t,*r); |
} |
} |
} |
} |
|
|
void addlf(GZ a,GZ b,GZ *c) |
void addlf(GZ a,GZ b,GZ *c) |
Line 1274 void addlf(GZ a,GZ b,GZ *c) |
|
Line 1275 void addlf(GZ a,GZ b,GZ *c) |
|
|
|
addgz(a,b,c); |
addgz(a,b,c); |
if ( !lf_lazy ) { |
if ( !lf_lazy ) { |
remgz(*c,current_mod_lf,&t); *c = t; |
// remgz(*c,current_mod_lf,&t); *c = t; |
|
if ( cmpgz(*c,current_mod_lf) >= 0 ) { |
|
subgz(*c,current_mod_lf,&t); *c = t; |
|
} |
} |
} |
} |
} |
|
|
Line 1323 void lmtolf(LM a,GZ *b) |
|
Line 1327 void lmtolf(LM a,GZ *b) |
|
{ |
{ |
Q q; |
Q q; |
|
|
NTOQ(BDY(a),1,q); *b = ztogz(q); |
if ( !a ) *b = 0; |
|
else { |
|
NTOQ(BDY(a),1,q); *b = ztogz(q); |
|
} |
} |
} |
|
|
void setmod_lf(N p) |
void setmod_lf(N p) |
Line 1331 void setmod_lf(N p) |
|
Line 1338 void setmod_lf(N p) |
|
Q q; |
Q q; |
|
|
NTOQ(p,1,q); current_mod_lf = ztogz(q); |
NTOQ(p,1,q); current_mod_lf = ztogz(q); |
|
current_mod_lf_size = mpz_size(BDY(current_mod_lf))+1; |
} |
} |
|
|
void simplf_force(GZ a,GZ *b) |
void simplf_force(GZ a,GZ *b) |