===================================================================
RCS file: /home/cvs/OpenXM_contrib2/asir2018/engine/Q.c,v
retrieving revision 1.9
retrieving revision 1.22
diff -u -p -r1.9 -r1.22
--- OpenXM_contrib2/asir2018/engine/Q.c	2018/10/19 23:27:38	1.9
+++ 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.8 2018/10/02 09:06:15 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)
 {
@@ -30,7 +35,7 @@ void gc_free(void *p,size_t size)
 
 void init_gmpq()
 {
-  mp_set_memory_functions(Risa_GC_malloc_atomic,gc_realloc,gc_free);
+  mp_set_memory_functions(Risa_GC_malloc,gc_realloc,gc_free);
 
   mpz_init(ONEMPZ); mpz_set_ui(ONEMPZ,1); MPZTOZ(ONEMPZ,ONE);
   gmp_randinit_default(GMP_RAND);
@@ -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)
@@ -515,7 +525,7 @@ void factorialz(unsigned int n,Z *nr)
 {
   mpz_t a;
   mpz_init(a);
-  mpz_fac_ui(a,n);
+  mpz_fac_ui(a,(unsigned long)n);
   MPZTOZ(a,*nr);
 }
 
@@ -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 {
@@ -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,27 +1384,37 @@ 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;
   }
 }
 
 extern int f4_nocheck;
 
-int mpz_gensolve_check(MAT mat,mpz_t **nm,mpz_t dn,int rank,int *rind,int *cind)
+int mpz_gensolve_check(MAT mat,mpz_t **nm,mpz_t dn,int rank,int clen,int *rind,int *cind)
 {
-  int row,col,clen,i,j,k,l;
+  int row,col,i,j,k,l;
   mpz_t t;
   mpz_t *w;
   Z *mati;
   mpz_t *nmk;
 
   if ( f4_nocheck ) return 1;
-  row = mat->row; col = mat->col; clen = col-rank;
+  row = mat->row; col = mat->col;
   w = (mpz_t *)MALLOC(clen*sizeof(mpz_t));
   mpz_init(t);
   for ( i = 0; i < clen; i++ ) mpz_init(w[i]);
@@ -1469,7 +1505,7 @@ void isqrtz(Z a,Z *r)
   Z two;
 
   if ( !a ) *r = 0;
-  else if ( UNIQ(a) ) *r = ONE;
+  else if ( UNIZ(a) ) *r = ONE;
   else {
     k = z_bits((Q)a); /* a <= 2^k-1 */
     bshiftz(ONE,-((k>>1)+(k&1)),&x); /* a <= x^2 */
@@ -1592,7 +1628,7 @@ int generic_gauss_elim_hensel(MAT mat,MAT *nmmat,Z *dn
   Z wq;
 
 #if SIZEOF_LONG == 8
-  return generic_gauss_elim_hensel64(mat,nmmat,dn,rindp,cindp);
+  return generic_gauss_elim_hensel64(mat,nmmat,dn,rindp,cindp,0);
 #endif
 init_eg(&eg_mul1); init_eg(&eg_mul2);
   a0 = (Z **)mat->body;
@@ -1706,8 +1742,10 @@ get_eg(&tmp2); add_eg(&eg_mul2,&tmp1,&tmp2);
                 if ( !cinfo[j] )
                   cind[k++] = j;
               return rank;
-            }
+            } else
+              goto reset;
           } else {
+reset:
             period = period*3/2;
             count = 0;
           }
@@ -1743,6 +1781,9 @@ int generic_gauss_elim_hensel_dalg(MAT mat,DP *mb,MAT 
   Z wq;
   DP m;
 
+#if SIZEOF_LONG == 8
+  return generic_gauss_elim_hensel64(mat,nmmat,dn,rindp,cindp,mb);
+#endif
   a0 = (Z **)mat->body;
   row = mat->row; col = mat->col;
   w = (int **)almat(row,col);
@@ -1815,8 +1856,8 @@ int generic_gauss_elim_hensel_dalg(MAT mat,DP *mb,MAT 
         for ( i = 0; i < rank; i++ )
           for ( j = 0, bi = b[i], wi = wc[i]; j < ri; j++ )
             wi[j] = remqi((Q)bi[j],md);
-        /* wc = A^(-1)wc; wc is normalized */
-        solve_by_lu_mod(w,rank,md,wc,ri,1);
+        /* wc = A^(-1)wc; wc is not normalized */
+        solve_by_lu_mod(w,rank,md,wc,ri,0);
         /* x += q*wc */
         for ( i = 0; i < rank; i++ )
           for ( j = 0, wi = wc[i]; j < ri; j++ ) mul1addtoz(q,wi[j],&x[i][j]);
@@ -1829,17 +1870,13 @@ int generic_gauss_elim_hensel_dalg(MAT mat,DP *mb,MAT 
               mpz_init_set(uz,BDY(b[i][j]));
             else
               mpz_init_set_ui(uz,0);
-            for ( k = 0; k < rank; k++ ) {
-              if ( a[i][k] && wc[k][j] ) {
-                if ( wc[k][j] < 0 )
-                  mpz_addmul_ui(uz,BDY(a[i][k]),-wc[k][j]);
-                else
+            for ( k = 0; k < rank; k++ )
+              if ( a[i][k] && wc[k][j] )
                   mpz_submul_ui(uz,BDY(a[i][k]),wc[k][j]);
-              }
-            }
             MPZTOZ(uz,u);
             divsz(u,mdq,&b[i][j]);
           }
+
         count++;
         /* q = q*md */
         mulz(q,mdq,&u); q = u;
@@ -1859,8 +1896,10 @@ int generic_gauss_elim_hensel_dalg(MAT mat,DP *mb,MAT 
                 if ( !cinfo[j] )
                   cind[k++] = j;
               return rank;
-            }
+            } else
+              goto reset;
           } else {
+reset:
             period = period*3/2;
             count = 0;
           }
@@ -2017,7 +2056,7 @@ RESET:
         for ( j = k = l = 0; j < col; j++ )
           if ( colstat[j] ) rind[k++] = j;  
           else cind[l++] = j;
-        if ( mpz_gensolve_check(mat,num,den,rank,rind,cind) ) {
+        if ( mpz_gensolve_check(mat,num,den,rank,col-rank,rind,cind) ) {
           MKMAT(r,rank,col-rank); *nm = r;
           for ( i = 0; i < rank; i++ )
             for ( j = 0; j < col-rank; j++ ) {
@@ -2131,7 +2170,7 @@ RESET:
 }
 #endif
 
-int generic_gauss_elim_hensel64(MAT mat,MAT *nmmat,Z *dn,int **rindp,int **cindp)
+int generic_gauss_elim_hensel64(MAT mat,MAT *nmmat,Z *dn,int **rindp,int **cindp,DP *mb)
 {
   MAT r;
   Z z;
@@ -2151,6 +2190,7 @@ int generic_gauss_elim_hensel64(MAT mat,MAT *nmmat,Z *
   int count;
   int ret;
   int period;
+  DP m;
 
   a0 = (Z **)mat->body;
   row = mat->row; col = mat->col;
@@ -2168,10 +2208,23 @@ int generic_gauss_elim_hensel64(MAT mat,MAT *nmmat,Z *
     if ( DP_Print > 3 ) {
       fprintf(asir_out,"done.\n"); fflush(asir_out);
     }
+
+    if ( mb ) {
+      /* this part is added for inv_or_split_dalg */
+      for ( i = 0; i < col-1; i++ ) {
+        if ( !cinfo[i] ) {
+          m = mb[i];
+          for ( j = i+1; j < col-1; j++ )
+            if ( dp_redble(mb[j],m) )
+              cinfo[j] = -1;
+        }
+      }
+    }
+  
     a = (mpz_t **)mpz_allocmat(rank,rank); /* lhs mat */
     b = (mpz_t **)mpz_allocmat(rank,col-rank);
     for ( j = li = ri = 0; j < col; j++ )
-      if ( cinfo[j] ) {
+      if ( cinfo[j] > 0 ) {
         /* the column is in lhs */
         for ( i = 0; i < rank; i++ ) {
           w[i][li] = w[i][j];
@@ -2181,7 +2234,7 @@ int generic_gauss_elim_hensel64(MAT mat,MAT *nmmat,Z *
             mpz_set_ui(a[i][li],0);
         }
         li++;
-      } else {
+      } else if ( !cinfo[j] ) {
         /* the column is in rhs */
         for ( i = 0; i < rank; i++ ) {
           if ( a0[rinfo[i]][j] )
@@ -2224,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++ )
@@ -2249,16 +2302,15 @@ int generic_gauss_elim_hensel64(MAT mat,MAT *nmmat,Z *
         count++;
         /* q = q*md */
         mpz_mul_ui(q,q,md);
-        fprintf(stderr,".");
         if ( count == period ) {
           ret = mpz_intmtoratm(x,rank,ri,q,nm,den);
           if ( ret ) {
             for ( j = k = l = 0; j < col; j++ )
-              if ( cinfo[j] )
+              if ( cinfo[j] > 0 )
                 rind[k++] = j;  
-              else
+              else if ( !cinfo[j] )
                 cind[l++] = j;
-            ret = mpz_gensolve_check(mat,nm,den,rank,rind,cind);
+            ret = mpz_gensolve_check(mat,nm,den,rank,ri,rind,cind);
             if ( ret ) {
               *rindp = rind;
               *cindp = cind;
@@ -2272,8 +2324,10 @@ int generic_gauss_elim_hensel64(MAT mat,MAT *nmmat,Z *
                 }
               MPZTOZ(den,*dn);
               return rank;
-            }
+            } else
+              goto reset;
           } else {
+reset:
             fprintf(stderr,"F");
             period = period*3/2;
             count = 0;