[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.3 and 1.5

version 1.3, 2018/09/24 22:26:43 version 1.5, 2018/09/28 08:20:28
Line 353  void pwrz(Z n1,Z n,Z *nr)
Line 353  void pwrz(Z n1,Z n,Z *nr)
   } else if ( !smallz(n) ) {    } else if ( !smallz(n) ) {
     error("exponent too big."); *nr = 0;      error("exponent too big."); *nr = 0;
   } else if ( n1->z && mpz_sgn(BDY((Z)n))>0 ) {    } else if ( n1->z && mpz_sgn(BDY((Z)n))>0 ) {
     mpz_init(z); mpz_pow_ui(z,BDY(n1),QTOS(n)); MPZTOZ(z,*nr);      mpz_init(z); mpz_pow_ui(z,BDY(n1),ZTOS(n)); MPZTOZ(z,*nr);
   } else {    } else {
     MPZTOMPQ(BDY(n1),q); MPQTOQ(q,r);      MPZTOMPQ(BDY(n1),q); MPQTOQ(q,r);
     pwrq(r,(Q)n,&p); *nr = (Z)p;      pwrq(r,(Q)n,&p); *nr = (Z)p;
Line 453  void gcdvz_estimate(VECT v,Z *q)
Line 453  void gcdvz_estimate(VECT v,Z *q)
     else addz(s,b[i],&u);      else addz(s,b[i],&u);
     s = u;      s = u;
   }    }
   for ( i = 0, t = 0; i < n; i++ ) {    for ( t = 0; i < n; i++ ) {
     if ( b[i] && mpz_sgn(BDY(b[i]))<0 ) subz(t,b[i],&u);      if ( b[i] && mpz_sgn(BDY(b[i]))<0 ) subz(t,b[i],&u);
     else addz(t,b[i],&u);      else addz(t,b[i],&u);
     t = u;      t = u;
Line 461  void gcdvz_estimate(VECT v,Z *q)
Line 461  void gcdvz_estimate(VECT v,Z *q)
   gcdz(s,t,q);    gcdz(s,t,q);
 }  }
   
   void gcdv_mpz_estimate(mpz_t g,mpz_t *b,int n)
   {
     int m,m2,i,j;
     mpz_t s,t;
   
     mpz_init(g);
     for ( i = 0, m = 0; i < n; i++ )
       if ( mpz_sgn(b[i]) ) m++;
     if ( !m ) {
       mpz_set_ui(g,0);
       return;
     }
     if ( m == 1 ) {
       for ( i = 0, m = 0; i < n; i++ )
         if ( mpz_sgn(b[i]) ) break;
       if ( mpz_sgn(b[i])<0 ) mpz_neg(g,b[i]);
       else mpz_set(g,b[i]);
       return ;
     }
     m2 = m/2;
     mpz_init_set_ui(s,0);
     for ( i = j = 0; j < m2; i++ ) {
       if ( mpz_sgn(b[i]) ) {
         if ( mpz_sgn(b[i])<0 )
           mpz_sub(s,s,b[i]);
         else
           mpz_add(s,s,b[i]);
         j++;
       }
     }
     mpz_init_set_ui(t,0);
     for ( ; i < n; i++ ) {
       if ( mpz_sgn(b[i]) ) {
         if ( mpz_sgn(b[i])<0 )
           mpz_sub(t,t,b[i]);
         else
           mpz_add(t,t,b[i]);
       }
     }
     mpz_gcd(g,s,t);
   }
   
   
 void factorialz(unsigned int n,Z *nr)  void factorialz(unsigned int n,Z *nr)
 {  {
   mpz_t a;    mpz_t a;
Line 605  void pwrq(Q n1,Q n,Q *nr)
Line 648  void pwrq(Q n1,Q n,Q *nr)
   } else if ( !smallz((Z)n) ) {    } else if ( !smallz((Z)n) ) {
     error("exponent too big."); *nr = 0;      error("exponent too big."); *nr = 0;
   } else {    } else {
     e = QTOS(n);      e = ZTOS(n);
     if ( e < 0 ) {      if ( e < 0 ) {
       e = -e;        e = -e;
       if ( n1->z ) {        if ( n1->z ) {
Line 662  void mkbc(int n,Z *t)
Line 705  void mkbc(int n,Z *t)
   Z c,d,iq;    Z c,d,iq;
   
   for ( t[0] = ONE, i = 1; i <= n/2; i++ ) {    for ( t[0] = ONE, i = 1; i <= n/2; i++ ) {
     STOQ(n-i+1,c); mulz(t[i-1],c,&d);      STOZ(n-i+1,c); mulz(t[i-1],c,&d);
     STOQ(i,iq); divsz(d,iq,&t[i]);      STOZ(i,iq); divsz(d,iq,&t[i]);
   }    }
   for ( ; i <= n; i++ )    for ( ; i <= n; i++ )
     t[i] = t[n-i];      t[i] = t[n-i];
Line 1326  void isqrtz(Z a,Z *r)
Line 1369  void isqrtz(Z a,Z *r)
   else {    else {
     k = z_bits((Q)a); /* a <= 2^k-1 */      k = z_bits((Q)a); /* a <= 2^k-1 */
     bshiftz(ONE,-((k>>1)+(k&1)),&x); /* a <= x^2 */      bshiftz(ONE,-((k>>1)+(k&1)),&x); /* a <= x^2 */
     STOQ(2,two);      STOZ(2,two);
     while ( 1 ) {      while ( 1 ) {
       pwrz(x,two,&t);        pwrz(x,two,&t);
       if ( cmpz(t,a) <= 0 ) {        if ( cmpz(t,a) <= 0 ) {
Line 1450  init_eg(&eg_mul1); init_eg(&eg_mul2);
Line 1493  init_eg(&eg_mul1); init_eg(&eg_mul2);
   w = (int **)almat(row,col);    w = (int **)almat(row,col);
   for ( ind = 0; ; ind++ ) {    for ( ind = 0; ; ind++ ) {
     md = get_lprime(ind);      md = get_lprime(ind);
     STOQ(md,mdq);      STOZ(md,mdq);
     for ( i = 0; i < row; i++ )      for ( i = 0; i < row; i++ )
       for ( j = 0, ai = a0[i], wi = w[i]; j < col; j++ )        for ( j = 0, ai = a0[i], wi = w[i]; j < col; j++ )
         wi[j] = remqi((Q)ai[j],md);          wi[j] = remqi((Q)ai[j],md);
Line 1598  int generic_gauss_elim_hensel_dalg(MAT mat,DP *mb,MAT 
Line 1641  int generic_gauss_elim_hensel_dalg(MAT mat,DP *mb,MAT 
   w = (int **)almat(row,col);    w = (int **)almat(row,col);
   for ( ind = 0; ; ind++ ) {    for ( ind = 0; ; ind++ ) {
     md = get_lprime(ind);      md = get_lprime(ind);
     STOQ(md,mdq);      STOZ(md,mdq);
     for ( i = 0; i < row; i++ )      for ( i = 0; i < row; i++ )
       for ( j = 0, ai = a0[i], wi = w[i]; j < col; j++ )        for ( j = 0, ai = a0[i], wi = w[i]; j < col; j++ )
         wi[j] = remqi((Q)ai[j],md);          wi[j] = remqi((Q)ai[j],md);
Line 1624  int generic_gauss_elim_hensel_dalg(MAT mat,DP *mb,MAT 
Line 1667  int generic_gauss_elim_hensel_dalg(MAT mat,DP *mb,MAT 
     a = (Z **)almat_pointer(rank,rank); /* lhs mat */      a = (Z **)almat_pointer(rank,rank); /* lhs mat */
     MKMAT(bmat,rank,col-rank); b = (Z **)bmat->body; /* lhs mat */      MKMAT(bmat,rank,col-rank); b = (Z **)bmat->body; /* lhs mat */
     for ( j = li = ri = 0; j < col; j++ )      for ( j = li = ri = 0; j < col; j++ )
       if ( cinfo[j] ) {        if ( cinfo[j] > 0 ) {
         /* the column is in lhs */          /* the column is in lhs */
         for ( i = 0; i < rank; i++ ) {          for ( i = 0; i < rank; i++ ) {
           w[i][li] = w[i][j];            w[i][li] = w[i][j];
           a[i][li] = a0[rinfo[i]][j];            a[i][li] = a0[rinfo[i]][j];
         }          }
         li++;          li++;
       } else {        } else if ( !cinfo[j] ) {
         /* the column is in rhs */          /* the column is in rhs */
         for ( i = 0; i < rank; i++ )          for ( i = 0; i < rank; i++ )
           b[i][ri] = a0[rinfo[i]][j];            b[i][ri] = a0[rinfo[i]][j];
Legend:
Removed from v.1.3  
changed lines
  Added in v.1.5
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>