=================================================================== RCS file: /home/cvs/OpenXM_contrib2/asir2018/engine/Q.c,v retrieving revision 1.17 retrieving revision 1.21 diff -u -p -r1.17 -r1.21 --- OpenXM_contrib2/asir2018/engine/Q.c 2020/01/09 01:47:41 1.17 +++ OpenXM_contrib2/asir2018/engine/Q.c 2021/10/08 07:27:05 1.21 @@ -1,4 +1,4 @@ -/* $OpenXM: OpenXM_contrib2/asir2018/engine/Q.c,v 1.16 2019/05/22 02:41:16 noro Exp $ */ +/* $OpenXM: OpenXM_contrib2/asir2018/engine/Q.c,v 1.20 2020/12/04 08:09:33 noro Exp $ */ #include "ca.h" #include "gmp.h" #include "base.h" @@ -17,6 +17,10 @@ extern int DP_Print; void isqrtz(Z a,Z *r); void bshiftz(Z a,int n,Z *r); +int mpz_inttorat(mpz_t c,mpz_t m,mpz_t b,mpz_t nm,mpz_t dn); +int generic_gauss_elim_hensel64(MAT mat,MAT *nmmat,Z *dn,int **rindp,int **cindp,DP *mb); +int find_lhs_and_lu_mod64(mp_limb_t **a,int row,int col,mp_limb_t md,int **rinfo,int **cinfo); +void solve_by_lu_mod64(mp_limb_t **a,int n,mp_limb_t md,mp_limb_signed_t **b,int l,int normalize); void *gc_realloc(void *p,size_t osize,size_t nsize) { @@ -93,8 +97,11 @@ void dnq(Q q,Z *r) int sgnq(Q q) { if ( !q ) return 0; - else if ( q->z ) return mpz_sgn(BDY((Z)q)); - else return mpz_sgn(mpq_numref(BDY(q))); + else if ( q->z ) { + int sgn; + sgn = mpz_sgn(BDY((Z)q)); + return sgn; + } else return mpz_sgn(mpq_numref(BDY(q))); } Q mpqtozq(mpq_t a) @@ -388,6 +395,8 @@ int cmpz(Z q1,Z q2) else if ( sgn < 0 ) return -1; else return 0; } + /* XXX */ + return 0; } void gcdz(Z n1,Z n2,Z *nq) @@ -665,6 +674,10 @@ void pwrq(Q n1,Q n,Q *nr) nm[0] = mpq_denref(BDY(n1))[0]; dn[0] = mpq_numref(BDY(n1))[0]; } + if ( mpz_sgn(dn)<0 ) { + mpz_neg(nm,nm); + mpz_neg(dn,dn); + } } else { if ( n1->z ) { nm[0] = BDY((Z)n1)[0]; @@ -783,10 +796,10 @@ void lgp(P p,Z *g,Z *l) if ( NUM(p) ) { if ( ((Q)p)->z ) { - MPZTOZ(BDY((Z)p),*g); + absz((Z)p,g); *l = ONE; } else { - MPZTOZ(mpq_numref(BDY((Q)p)),*g); + MPZTOZ(mpq_numref(BDY((Q)p)),g1); absz(g1,g); MPZTOZ(mpq_denref(BDY((Q)p)),*l); } } else { @@ -1232,7 +1245,7 @@ int intmtoratm(MAT mat,Z md,MAT nm,Z *dn) return 0; row = mat->row; col = mat->col; bshiftz(md,1,&t); - isqrt(t,&s); + isqrtz(t,&s); bshiftz(s,64,&b); if ( !b ) b = ONE; dn0 = ONE; @@ -2263,7 +2276,7 @@ int generic_gauss_elim_hensel64(MAT mat,MAT *nmmat,Z * for ( j = 0, bi = b[i], wi = wc[i]; j < ri; j++ ) wi[j] = mpz_fdiv_ui(bi[j],md); /* wc = A^(-1)wc; wc is not normalized */ - solve_by_lu_mod64(w,rank,md,wc,ri,0); + solve_by_lu_mod64(w,rank,md,(mp_limb_signed_t **)wc,ri,0); /* x += q*wc */ for ( i = 0; i < rank; i++ ) for ( j = 0, wi = wc[i]; j < ri; j++ )