===================================================================
RCS file: /home/cvs/OpenXM_contrib2/asir2018/engine/Q.c,v
retrieving revision 1.13
retrieving revision 1.22
diff -u -p -r1.13 -r1.22
--- OpenXM_contrib2/asir2018/engine/Q.c	2019/03/18 07:00:33	1.13
+++ OpenXM_contrib2/asir2018/engine/Q.c	2022/04/03 01:30:27	1.22
@@ -1,11 +1,12 @@
-/* $OpenXM: OpenXM_contrib2/asir2018/engine/Q.c,v 1.12 2019/03/17 02:48:13 noro Exp $ */
+/* $OpenXM: OpenXM_contrib2/asir2018/engine/Q.c,v 1.21 2021/10/08 07:27:05 noro Exp $ */
 #include "ca.h"
+#include "parse.h"
 #include "gmp.h"
 #include "base.h"
 #include "inline.h"
 
 mpz_t ONEMPZ;
-Z ONE;
+extern Z ONE;
 int lf_lazy;
 Z current_mod_lf;
 int current_mod_lf_size;
@@ -17,6 +18,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 +98,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 +396,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)
@@ -558,7 +568,7 @@ void subq(Q n1,Q n2,Q *nr)
 
   if ( !n1 ) {
     if ( !n2 ) *nr = 0;
-    else if ( n1->z ) chsgnz((Z)n1,(Z *)nr);
+    else if ( n2->z ) chsgnz((Z)n2,(Z *)nr);
     else {
         mpq_init(t); mpq_neg(t,BDY(n2)); MPQTOQ(t,*nr);
       }
@@ -665,6 +675,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];
@@ -762,12 +776,15 @@ void lgp(P p,Z *g,Z *l);
 
 void ptozp(P p,int sgn,Q *c,P *pr)
 {
-  Z nm,dn;
+  Z nm,dn,nm1;
 
   if ( !p ) {
     *c = 0; *pr = 0;
   } else {
     lgp(p,&nm,&dn);
+    if ( sgn < 0 ) {
+      chsgnz(nm,&nm1); nm = nm1;
+    }
     divz(nm,dn,(Z *)c);
     divsp(CO,p,(P)*c,pr);
   }
@@ -780,10 +797,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 {
@@ -1229,7 +1246,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;
@@ -1337,6 +1354,15 @@ int mpz_inttorat(mpz_t c,mpz_t m,mpz_t b,mpz_t nm,mpz_
   }
   if ( mpz_cmp(v1,b) >= 0 ) return 0;
   else {
+    mpz_gcd(t,v1,v2);
+    if ( UNIMPZ(t) ) 
+      mpz_set_ui(r1,0); 
+    else {
+      /* v1 /= t, v2 /= t, t=c*v1-v2, r1=t%m */
+      mpz_divexact(v1,v1,t); mpz_divexact(v2,v2,t);
+      mpz_mul(t,c,v1); mpz_sub(t,t,v2); mpz_mod(r1,t,m);
+    }
+    if ( mpz_sgn(r1) ) return 0;
     if ( mpz_sgn(v1)<0  ) {
       mpz_neg(dn,v1); mpz_neg(nm,v2);
     } else {
@@ -1348,7 +1374,7 @@ int mpz_inttorat(mpz_t c,mpz_t m,mpz_t b,mpz_t nm,mpz_
 
 int inttorat(Z c,Z m,Z b,Z *nmp,Z *dnp)
 {
-  Z qq,t,u1,v1,r1;
+  Z qq,t,s,r,u1,v1,r1;
   Z q,u2,v2,r2;
 
   u1 = 0; v1 = ONE; u2 = m; v2 = c;
@@ -1358,11 +1384,21 @@ int inttorat(Z c,Z m,Z b,Z *nmp,Z *dnp)
   }
   if ( cmpz(v1,b) >= 0 ) return 0;
   else {
-    if ( mpz_sgn(BDY(v1))<0  ) {
-      chsgnz(v1,dnp); chsgnz(v2,nmp);
+    /* reduction and check */
+    /* v2/v1 = u2/u1, c*u1-u2 = 0 mod m? */
+    gcdz(v1,v2,&t);
+    if ( UNIZ(t) ) {
+      u1 = v1; u2 = v2; r = 0;
     } else {
-      *dnp = v1; *nmp = v2;
+      divsz(v1,t,&u1); divsz(v2,t,&u2);
+      mulz(c,u1,&t); subz(t,u2,&s); remz(s,m,&r);
     }
+    if ( r ) return 0;
+    if ( mpz_sgn(BDY(u1))<0  ) {
+      chsgnz(u1,dnp); chsgnz(u2,nmp);
+    } else {
+      *dnp = u1; *nmp = u2;
+    }
     return 1;
   }
 }
@@ -2241,7 +2277,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++ )