[BACK]Return to Q.c CVS log [TXT][DIR] Up to [local] / OpenXM_contrib2 / asir2018 / engine

Diff for /OpenXM_contrib2/asir2018/engine/Q.c between version 1.11 and 1.17

version 1.11, 2019/03/03 05:21:17 version 1.17, 2020/01/09 01:47:41
Line 1 
Line 1 
 /* $OpenXM: OpenXM_contrib2/asir2018/engine/Q.c,v 1.10 2018/12/10 22:24:42 noro Exp $ */  /* $OpenXM: OpenXM_contrib2/asir2018/engine/Q.c,v 1.16 2019/05/22 02:41:16 noro Exp $ */
 #include "ca.h"  #include "ca.h"
 #include "gmp.h"  #include "gmp.h"
 #include "base.h"  #include "base.h"
 #include "inline.h"  #include "inline.h"
   
 mpz_t ONEMPZ;  mpz_t ONEMPZ;
 Z ONE;  extern Z ONE;
 int lf_lazy;  int lf_lazy;
 Z current_mod_lf;  Z current_mod_lf;
 int current_mod_lf_size;  int current_mod_lf_size;
Line 515  void factorialz(unsigned int n,Z *nr)
Line 515  void factorialz(unsigned int n,Z *nr)
 {  {
   mpz_t a;    mpz_t a;
   mpz_init(a);    mpz_init(a);
   mpz_fac_ui(a,n);    mpz_fac_ui(a,(unsigned long)n);
   MPZTOZ(a,*nr);    MPZTOZ(a,*nr);
 }  }
   
Line 558  void subq(Q n1,Q n2,Q *nr)
Line 558  void subq(Q n1,Q n2,Q *nr)
   
   if ( !n1 ) {    if ( !n1 ) {
     if ( !n2 ) *nr = 0;      if ( !n2 ) *nr = 0;
     else if ( n1->z ) chsgnz((Z)n1,(Z *)nr);      else if ( n2->z ) chsgnz((Z)n2,(Z *)nr);
     else {      else {
         mpq_init(t); mpq_neg(t,BDY(n2)); MPQTOQ(t,*nr);          mpq_init(t); mpq_neg(t,BDY(n2)); MPQTOQ(t,*nr);
       }        }
Line 762  void lgp(P p,Z *g,Z *l);
Line 762  void lgp(P p,Z *g,Z *l);
   
 void ptozp(P p,int sgn,Q *c,P *pr)  void ptozp(P p,int sgn,Q *c,P *pr)
 {  {
   Z nm,dn;    Z nm,dn,nm1;
   
   if ( !p ) {    if ( !p ) {
     *c = 0; *pr = 0;      *c = 0; *pr = 0;
   } else {    } else {
     lgp(p,&nm,&dn);      lgp(p,&nm,&dn);
       if ( sgn < 0 ) {
         chsgnz(nm,&nm1); nm = nm1;
       }
     divz(nm,dn,(Z *)c);      divz(nm,dn,(Z *)c);
     divsp(CO,p,(P)*c,pr);      divsp(CO,p,(P)*c,pr);
   }    }
Line 1337  int mpz_inttorat(mpz_t c,mpz_t m,mpz_t b,mpz_t nm,mpz_
Line 1340  int mpz_inttorat(mpz_t c,mpz_t m,mpz_t b,mpz_t nm,mpz_
   }    }
   if ( mpz_cmp(v1,b) >= 0 ) return 0;    if ( mpz_cmp(v1,b) >= 0 ) return 0;
   else {    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  ) {      if ( mpz_sgn(v1)<0  ) {
       mpz_neg(dn,v1); mpz_neg(nm,v2);        mpz_neg(dn,v1); mpz_neg(nm,v2);
     } else {      } else {
Line 1348  int mpz_inttorat(mpz_t c,mpz_t m,mpz_t b,mpz_t nm,mpz_
Line 1360  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)  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;    Z q,u2,v2,r2;
   
   u1 = 0; v1 = ONE; u2 = m; v2 = c;    u1 = 0; v1 = ONE; u2 = m; v2 = c;
Line 1358  int inttorat(Z c,Z m,Z b,Z *nmp,Z *dnp)
Line 1370  int inttorat(Z c,Z m,Z b,Z *nmp,Z *dnp)
   }    }
   if ( cmpz(v1,b) >= 0 ) return 0;    if ( cmpz(v1,b) >= 0 ) return 0;
   else {    else {
     if ( mpz_sgn(BDY(v1))<0  ) {      /* reduction and check */
       chsgnz(v1,dnp); chsgnz(v2,nmp);      /* v2/v1 = u2/u1, c*u1-u2 = 0 mod m? */
       gcdz(v1,v2,&t);
       if ( UNIZ(t) ) {
         u1 = v1; u2 = v2; r = 0;
     } else {      } 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;      return 1;
   }    }
 }  }
   
 extern int f4_nocheck;  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 t;
   mpz_t *w;    mpz_t *w;
   Z *mati;    Z *mati;
   mpz_t *nmk;    mpz_t *nmk;
   
   if ( f4_nocheck ) return 1;    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));    w = (mpz_t *)MALLOC(clen*sizeof(mpz_t));
   mpz_init(t);    mpz_init(t);
   for ( i = 0; i < clen; i++ ) mpz_init(w[i]);    for ( i = 0; i < clen; i++ ) mpz_init(w[i]);
Line 2020  RESET:
Line 2042  RESET:
         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 ( 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;            MKMAT(r,rank,col-rank); *nm = r;
           for ( i = 0; i < rank; i++ )            for ( i = 0; i < rank; i++ )
             for ( j = 0; j < col-rank; j++ ) {              for ( j = 0; j < col-rank; j++ ) {
Line 2274  int generic_gauss_elim_hensel64(MAT mat,MAT *nmmat,Z *
Line 2296  int generic_gauss_elim_hensel64(MAT mat,MAT *nmmat,Z *
                 rind[k++] = j;                  rind[k++] = j;
               else if ( !cinfo[j] )                else if ( !cinfo[j] )
                 cind[l++] = 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 ) {              if ( ret ) {
               *rindp = rind;                *rindp = rind;
               *cindp = cind;                *cindp = cind;
Legend:
Removed from v.1.11  
changed lines
  Added in v.1.17
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>