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