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

Diff for /OpenXM_contrib2/asir2018/builtin/array.c between version 1.3 and 1.11

version 1.3, 2018/10/01 05:49:06 version 1.11, 2022/01/13 08:29:32
Line 45 
Line 45 
  * DEVELOPER SHALL HAVE NO LIABILITY IN CONNECTION WITH THE USE,   * DEVELOPER SHALL HAVE NO LIABILITY IN CONNECTION WITH THE USE,
  * PERFORMANCE OR NON-PERFORMANCE OF THE SOFTWARE.   * PERFORMANCE OR NON-PERFORMANCE OF THE SOFTWARE.
  *   *
  * $OpenXM: OpenXM_contrib2/asir2018/builtin/array.c,v 1.2 2018/09/28 08:20:27 noro Exp $   * $OpenXM: OpenXM_contrib2/asir2018/builtin/array.c,v 1.10 2022/01/13 08:15:02 noro Exp $
 */  */
 #include "ca.h"  #include "ca.h"
 #include "base.h"  #include "base.h"
Line 170  void solve_u(int *,ent **,int,int *,int);
Line 170  void solve_u(int *,ent **,int,int *,int);
 static int *ul,*ll;  static int *ul,*ll;
 static ent **u,**l;  static ent **u,**l;
 static int modulus;  static int modulus;
   #if defined(ANDROID)
   int getw(FILE *fp)
   {
           int x;
           return (fread((void *)&x, sizeof(x), 1, fp) == 1 ? x : EOF);
   }
   #endif
   
 void Plusolve_prep(NODE arg,Q *rp)  void Plusolve_prep(NODE arg,Q *rp)
 {  {
Line 615  void Pnewvect(NODE arg,VECT *rp)
Line 622  void Pnewvect(NODE arg,VECT *rp)
       return;        return;
     }      }
 #endif  #endif
     for ( i = 0, tn = BDY(list), vb = BDY(vect); tn; i++, tn = NEXT(tn) )      for ( i = 0, tn = BDY(list), vb = BDY(vect); tn && i<len; i++, tn = NEXT(tn) )
       vb[i] = (pointer)BDY(tn);        vb[i] = (pointer)BDY(tn);
   }    }
   *rp = vect;    *rp = vect;
Line 1164  void Pgeneric_gauss_elim(NODE arg,LIST *rp)
Line 1171  void Pgeneric_gauss_elim(NODE arg,LIST *rp)
   if ( is_hensel )    if ( is_hensel )
     rank = generic_gauss_elim_hensel(m,&nm,&dn,&ri,&ci);      rank = generic_gauss_elim_hensel(m,&nm,&dn,&ri,&ci);
   else    else
     rank = generic_gauss_elim64(m,&nm,&dn,&ri,&ci);      rank = generic_gauss_elim(m,&nm,&dn,&ri,&ci);
   t = col-rank;    t = col-rank;
   MKVECT(rind,rank);    MKVECT(rind,rank);
   MKVECT(cind,t);    MKVECT(cind,t);
Line 1243  void Pgeneric_gauss_elim_mod64(NODE arg,LIST *rp)
Line 1250  void Pgeneric_gauss_elim_mod64(NODE arg,LIST *rp)
   Z q;    Z q;
   mp_limb_t md;    mp_limb_t md;
   int i,j,k,l,row,col,t,rank;    int i,j,k,l,row,col,t,rank;
     Obj val;
     int asis = 0;
   
   asir_assert(ARG0(arg),O_MAT,"generic_gauss_elim_mod64");    asir_assert(ARG0(arg),O_MAT,"generic_gauss_elim_mod64");
   asir_assert(ARG1(arg),O_N,"generic_gauss_elim_mod64");    asir_assert(ARG1(arg),O_N,"generic_gauss_elim_mod64");
     if ( get_opt("asis",&val) && val ) asis = 1;
   m = (MAT)ARG0(arg); md = ZTOS((Z)ARG1(arg));    m = (MAT)ARG0(arg); md = ZTOS((Z)ARG1(arg));
   row = m->row; col = m->col; tmat = (Z **)m->body;    row = m->row; col = m->col; tmat = (Z **)m->body;
   wmat = (mp_limb_t **)almat64(row,col);    wmat = (mp_limb_t **)almat64(row,col);
Line 1258  void Pgeneric_gauss_elim_mod64(NODE arg,LIST *rp)
Line 1268  void Pgeneric_gauss_elim_mod64(NODE arg,LIST *rp)
     for ( j = 0; j < col; j++ )      for ( j = 0; j < col; j++ )
       wmat[i][j] = remqi64((Q)tmat[i][j],md);        wmat[i][j] = remqi64((Q)tmat[i][j],md);
   rank = generic_gauss_elim_mod64(wmat,row,col,md,colstat);    rank = generic_gauss_elim_mod64(wmat,row,col,md,colstat);
     if ( asis ) {
       MKMAT(mat,rank,col);
       tmat = (Z **)mat->body;
       for ( i = 0; i < rank; i++ )
         for ( j = 0; j < col; j++ ) {
           UTOZ(wmat[i][j],tmat[i][j]);
         }
       *rp = (LIST)mat;
       return;
     }
   
   MKVECT(rnum,rank);    MKVECT(rnum,rank);
   rnb = (Z *)rnum->body;    rnb = (Z *)rnum->body;
Line 1300  void Pgeneric_gauss_elim_mod(NODE arg,LIST *rp)
Line 1320  void Pgeneric_gauss_elim_mod(NODE arg,LIST *rp)
   Z q;    Z q;
   long mdl;    long mdl;
   int md,i,j,k,l,row,col,t,rank;    int md,i,j,k,l,row,col,t,rank;
     Obj val;
     int asis = 0;
   
   asir_assert(ARG0(arg),O_MAT,"generic_gauss_elim_mod");    asir_assert(ARG0(arg),O_MAT,"generic_gauss_elim_mod");
   asir_assert(ARG1(arg),O_N,"generic_gauss_elim_mod");    asir_assert(ARG1(arg),O_N,"generic_gauss_elim_mod");
     if ( get_opt("asis",&val) && val ) asis = 1;
 #if SIZEOF_LONG==8  #if SIZEOF_LONG==8
   mdl = ZTOS((Z)ARG1(arg));    mdl = ZTOS((Z)ARG1(arg));
   if ( mdl >= ((mp_limb_t)1)<<32 ) {    if ( mdl >= ((mp_limb_t)1)<<32 ) {
Line 1323  void Pgeneric_gauss_elim_mod(NODE arg,LIST *rp)
Line 1346  void Pgeneric_gauss_elim_mod(NODE arg,LIST *rp)
     for ( j = 0; j < col; j++ )      for ( j = 0; j < col; j++ )
       wmat[i][j] = remqi((Q)tmat[i][j],md);        wmat[i][j] = remqi((Q)tmat[i][j],md);
   rank = generic_gauss_elim_mod(wmat,row,col,md,colstat);    rank = generic_gauss_elim_mod(wmat,row,col,md,colstat);
     if ( asis ) {
       MKMAT(mat,rank,col);
       tmat = (Z **)mat->body;
       for ( i = 0; i < rank; i++ )
         for ( j = 0; j < col; j++ ) {
           UTOZ(wmat[i][j],tmat[i][j]);
         }
       *rp = (LIST)mat;
       return;
     }
   
   MKVECT(rnum,rank);    MKVECT(rnum,rank);
   rnb = (Z *)rnum->body;    rnb = (Z *)rnum->body;
Line 1423  int gauss_elim_mod(int **mat,int row,int col,int md)
Line 1456  int gauss_elim_mod(int **mat,int row,int col,int md)
     return -1;      return -1;
 }  }
   
 struct oEGT eg_mod,eg_elim,eg_elim1,eg_elim2,eg_chrem,eg_gschk,eg_intrat,eg_symb;  struct oEGT eg_mod,eg_elim,eg_elim1,eg_elim2,eg_chrem,eg_gschk,eg_intrat,eg_symb,eg_back,eg_fore;
 struct oEGT eg_conv;  struct oEGT eg_conv;
   
 #if 0  #if 0
Line 1536  void lu_dec_cr(MAT mat,MAT lu,Q *dn,int **perm)
Line 1569  void lu_dec_cr(MAT mat,MAT lu,Q *dn,int **perm)
   
 int f4_nocheck;  int f4_nocheck;
   
 #define ONE_STEP1  if ( zzz = *s ) { DMAR(zzz,hc,*tj,md,*tj) } tj++; s++;  #define ONE_STEP1  if ( ( zzz = *s ) != 0 ) { DMAR(zzz,hc,*tj,md,*tj) } tj++; s++;
   
 void reduce_reducers_mod(int **mat,int row,int col,int md)  void reduce_reducers_mod(int **mat,int row,int col,int md)
 {  {
Line 1552  void reduce_reducers_mod(int **mat,int row,int col,int
Line 1585  void reduce_reducers_mod(int **mat,int row,int col,int
     ind[i] = j;      ind[i] = j;
     for ( l = i-1; l >= 0; l-- ) {      for ( l = i-1; l >= 0; l-- ) {
       /* reduce mat[i] by mat[l] */        /* reduce mat[i] by mat[l] */
       if ( hc = t[ind[l]] ) {        if ( ( hc = t[ind[l]] ) != 0 ) {
         /* mat[i] = mat[i]-hc*mat[l] */          /* mat[i] = mat[i]-hc*mat[l] */
         j = ind[l];          j = ind[l];
         s = mat[l]+j;          s = mat[l]+j;
Line 1578  void reduce_reducers_mod(int **mat,int row,int col,int
Line 1611  void reduce_reducers_mod(int **mat,int row,int col,int
           ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1            ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
         }          }
         for ( ; k > 0; k-- ) {          for ( ; k > 0; k-- ) {
           if ( zzz = *s ) { DMAR(zzz,hc,*tj,md,*tj) } tj++; s++;            if ( ( zzz = *s ) != 0 ) { DMAR(zzz,hc,*tj,md,*tj) } tj++; s++;
         }          }
       }        }
     }      }
Line 1613  void pre_reduce_mod(int **mat,int row,int col,int nred
Line 1646  void pre_reduce_mod(int **mat,int row,int col,int nred
         DMAR(t[k],inv,0,md,t[k])          DMAR(t[k],inv,0,md,t[k])
     for ( l = i-1; l >= 0; l-- ) {      for ( l = i-1; l >= 0; l-- ) {
       /* reduce mat[i] by mat[l] */        /* reduce mat[i] by mat[l] */
       if ( hc = t[ind[l]] ) {        if ( ( hc = t[ind[l]] ) != 0 ) {
         /* mat[i] = mat[i]-hc*mat[l] */          /* mat[i] = mat[i]-hc*mat[l] */
         for ( k = ind[l], hc = md-hc, s = mat[l]+k, tk = t+k;          for ( k = ind[l], hc = md-hc, s = mat[l]+k, tk = t+k;
           k < col; k++, tk++, s++ )            k < col; k++, tk++, s++ )
Line 1627  void pre_reduce_mod(int **mat,int row,int col,int nred
Line 1660  void pre_reduce_mod(int **mat,int row,int col,int nred
     t = mat[i];      t = mat[i];
     for ( l = nred-1; l >= 0; l-- ) {      for ( l = nred-1; l >= 0; l-- ) {
       /* reduce mat[i] by mat[l] */        /* reduce mat[i] by mat[l] */
       if ( hc = t[ind[l]] ) {        if ( ( hc = t[ind[l]] ) != 0 ) {
         /* mat[i] = mat[i]-hc*mat[l] */          /* mat[i] = mat[i]-hc*mat[l] */
         for ( k = ind[l], hc = md-hc, s = mat[l]+k, tk = t+k;          for ( k = ind[l], hc = md-hc, s = mat[l]+k, tk = t+k;
           k < col; k++, tk++, s++ )            k < col; k++, tk++, s++ )
Line 1651  void reduce_sp_by_red_mod(int *sp,int **redmat,int *in
Line 1684  void reduce_sp_by_red_mod(int *sp,int **redmat,int *in
   /* reduce the spolys by redmat */    /* reduce the spolys by redmat */
   for ( i = nred-1; i >= 0; i-- ) {    for ( i = nred-1; i >= 0; i-- ) {
     /* reduce sp by redmat[i] */      /* reduce sp by redmat[i] */
     if ( hc = sp[ind[i]] ) {      if ( ( hc = sp[ind[i]] ) != 0 ) {
       /* sp = sp-hc*redmat[i] */        /* sp = sp-hc*redmat[i] */
       j = ind[i];        j = ind[i];
       hc = md-hc;        hc = md-hc;
       s = redmat[i]+j;        s = redmat[i]+j;
       tj = sp+j;        tj = sp+j;
       for ( k = col-j; k > 0; k-- ) {        for ( k = col-j; k > 0; k-- ) {
         if ( zzz = *s ) { DMAR(zzz,hc,*tj,md,*tj) } tj++; s++;          if ( ( zzz = *s ) != 0 ) { DMAR(zzz,hc,*tj,md,*tj) } tj++; s++;
       }        }
     }      }
   }    }
Line 1756  int generic_gauss_elim_mod64(mp_limb_t **mat,int row,i
Line 1789  int generic_gauss_elim_mod64(mp_limb_t **mat,int row,i
         *pk = mulmod64(*pk,inv,md);          *pk = mulmod64(*pk,inv,md);
     for ( i = rank+1; i < row; i++ ) {      for ( i = rank+1; i < row; i++ ) {
       t = mat[i];        t = mat[i];
       if ( a = t[j] )        if ( ( a = t[j] ) != 0 )
         red_by_vect64mod(md,t+j,pivot+j,md-a,col-j);          red_by_vect64mod(md,t+j,pivot+j,md-a,col-j);
     }      }
     rank++;      rank++;
Line 1766  int generic_gauss_elim_mod64(mp_limb_t **mat,int row,i
Line 1799  int generic_gauss_elim_mod64(mp_limb_t **mat,int row,i
       pivot = mat[l];        pivot = mat[l];
       for ( i = 0; i < l; i++ ) {        for ( i = 0; i < l; i++ ) {
         t = mat[i];          t = mat[i];
         if ( a = t[j] )          if ( ( a = t[j] ) != 0 )
           red_by_vect64mod(md,t+j,pivot+j,md-a,col-j);            red_by_vect64mod(md,t+j,pivot+j,md-a,col-j);
       }        }
       l--;        l--;
Line 1774  int generic_gauss_elim_mod64(mp_limb_t **mat,int row,i
Line 1807  int generic_gauss_elim_mod64(mp_limb_t **mat,int row,i
   return rank;    return rank;
 }  }
   
   int find_lhs_and_lu_mod64(mp_limb_t **a,int row,int col,
     mp_limb_t md,int **rinfo,int **cinfo)
   {
     int i,j,k,d;
     int *rp,*cp;
     mp_limb_t *t,*pivot;
     mp_limb_t inv,m;
   
     *rinfo = rp = (int *)MALLOC_ATOMIC(row*sizeof(int));
     *cinfo = cp = (int *)MALLOC_ATOMIC(col*sizeof(int));
     for ( i = 0; i < row; i++ )
       rp[i] = i;
     for ( k = 0, d = 0; k < col; k++ ) {
       for ( i = d; i < row && !a[i][k]; i++ );
       if ( i == row ) {
         cp[k] = 0;
         continue;
       } else
         cp[k] = 1;
       if ( i != d ) {
         j = rp[i]; rp[i] = rp[d]; rp[d] = j;
         t = a[i]; a[i] = a[d]; a[d] = t;
       }
       pivot = a[d];
       pivot[k] = inv = invmod64(pivot[k],md);
       for ( i = d+1; i < row; i++ ) {
         t = a[i];
         if ( (m = t[k]) != 0 ) {
           t[k] = mulmod64(inv,m,md);
           for ( j = k+1, m = md - t[k]; j < col; j++ )
             if ( pivot[j] ) {
               t[j] = muladdmod64(m,pivot[j],t[j],md);
             }
         }
       }
       d++;
     }
     return d;
   }
   
   int lu_mod64(mp_limb_t **a,int n,mp_limb_t md,int **rinfo)
   {
     int i,j,k;
     int *rp;
     mp_limb_t *t,*pivot;
     mp_limb_t inv,m;
   
     *rinfo = rp = (int *)MALLOC_ATOMIC(n*sizeof(int));
     for ( i = 0; i < n; i++ ) rp[i] = i;
     for ( k = 0; k < n; k++ ) {
       for ( i = k; i < n && !a[i][k]; i++ );
       if ( i == n ) return 0;
       if ( i != k ) {
         j = rp[i]; rp[i] = rp[k]; rp[k] = j;
         t = a[i]; a[i] = a[k]; a[k] = t;
       }
       pivot = a[k];
       inv = invmod64(pivot[k],md);
       for ( i = k+1; i < n; i++ ) {
         t = a[i];
         if ( (m = t[k]) != 0 ) {
           t[k] = mulmod64(inv,m,md);
           for ( j = k+1, m = md - t[k]; j < n; j++ )
             if ( pivot[j] ) {
               t[j] = muladdmod64(m,pivot[j],t[j],md);
             }
         }
       }
     }
     return 1;
   }
   
   /*
     Input
     a : n x n matrix; a result of LU-decomposition
     md : modulus
     b : n x l matrix
    Output
     b = a^(-1)b
    */
   
   void solve_by_lu_mod64(mp_limb_t **a,int n,mp_limb_t md,mp_limb_signed_t **b,int l,int normalize)
   {
     mp_limb_t *y,*c;
     int i,j,k;
     mp_limb_t t,m,m2;
   
     y = (mp_limb_t *)MALLOC_ATOMIC(n*sizeof(mp_limb_t));
     c = (mp_limb_t *)MALLOC_ATOMIC(n*sizeof(mp_limb_t));
     m2 = md/2;
     for ( k = 0; k < l; k++ ) {
       /* copy b[.][k] to c */
       for ( i = 0; i < n; i++ )
         c[i] = b[i][k];
       /* solve Ly=c */
       for ( i = 0; i < n; i++ ) {
         for ( t = c[i], j = 0; j < i; j++ )
           if ( a[i][j] ) {
             m = md - a[i][j];
             t = muladdmod64(m,y[j],t,md);
           }
         y[i] = t;
       }
       /* solve Uc=y */
       for ( i = n-1; i >= 0; i-- ) {
         for ( t = y[i], j =i+1; j < n; j++ )
           if ( a[i][j] ) {
             m = md - a[i][j];
             t = muladdmod64(m,c[j],t,md);
           }
         /* a[i][i] = 1/U[i][i] */
         c[i] = mulmod64(t,a[i][i],md);
       }
       /* copy c to b[.][k] with normalization */
       if ( normalize )
         for ( i = 0; i < n; i++ )
           b[i][k] = (mp_limb_signed_t)(c[i]>m2 ? c[i]-md : c[i]);
       else
         for ( i = 0; i < n; i++ )
           b[i][k] = (mp_limb_signed_t)c[i];
     }
   }
 #endif  #endif
   
 void red_by_vect_sf(int m,unsigned int *p,unsigned int *r,unsigned int hc,int len)  void red_by_vect_sf(int m,unsigned int *p,unsigned int *r,unsigned int hc,int len)
Line 1815  void reduce_sp_by_red_mod_compress (int *sp,CDP *redma
Line 1970  void reduce_sp_by_red_mod_compress (int *sp,CDP *redma
   for ( i = nred-1; i >= 0; i-- ) {    for ( i = nred-1; i >= 0; i-- ) {
     /* reduce sp by redmat[i] */      /* reduce sp by redmat[i] */
     usp[ind[i]] %= md;      usp[ind[i]] %= md;
     if ( hc = usp[ind[i]] ) {      if ( ( hc = usp[ind[i]] ) != 0 ) {
       /* sp = sp-hc*redmat[i] */        /* sp = sp-hc*redmat[i] */
       hc = md-hc;        hc = md-hc;
       ri = redmat[i];        ri = redmat[i];
Line 1861  int generic_gauss_elim_mod(int **mat0,int row,int col,
Line 2016  int generic_gauss_elim_mod(int **mat0,int row,int col,
       }        }
     for ( i = rank+1; i < row; i++ ) {      for ( i = rank+1; i < row; i++ ) {
       t = mat[i];        t = mat[i];
       if ( a = t[j] )        if ( ( a = t[j] ) != 0 )
         red_by_vect(md,t+j,pivot+j,md-a,col-j);          red_by_vect(md,t+j,pivot+j,md-a,col-j);
     }      }
     rank++;      rank++;
Line 1872  int generic_gauss_elim_mod(int **mat0,int row,int col,
Line 2027  int generic_gauss_elim_mod(int **mat0,int row,int col,
       for ( i = 0; i < l; i++ ) {        for ( i = 0; i < l; i++ ) {
         t = mat[i];          t = mat[i];
         t[j] %= md;          t[j] %= md;
         if ( a = t[j] )          if ( ( a = t[j] ) != 0 )
           red_by_vect(md,t+j,pivot+j,md-a,col-j);            red_by_vect(md,t+j,pivot+j,md-a,col-j);
       }        }
       l--;        l--;
Line 1921  int generic_gauss_elim_mod2(int **mat0,int row,int col
Line 2076  int generic_gauss_elim_mod2(int **mat0,int row,int col
       }        }
     for ( i = rank+1; i < row; i++ ) {      for ( i = rank+1; i < row; i++ ) {
       t = mat[i];        t = mat[i];
       if ( a = t[j] )        if ( ( a = t[j] ) != 0 )
         red_by_vect(md,t+j,pivot+j,md-a,col-j);          red_by_vect(md,t+j,pivot+j,md-a,col-j);
     }      }
     rank++;      rank++;
Line 1932  int generic_gauss_elim_mod2(int **mat0,int row,int col
Line 2087  int generic_gauss_elim_mod2(int **mat0,int row,int col
       for ( i = 0; i < l; i++ ) {        for ( i = 0; i < l; i++ ) {
         t = mat[i];          t = mat[i];
         t[j] %= md;          t[j] %= md;
         if ( a = t[j] )          if ( ( a = t[j] ) != 0 )
           red_by_vect(md,t+j,pivot+j,md-a,col-j);            red_by_vect(md,t+j,pivot+j,md-a,col-j);
       }        }
       l--;        l--;
Line 1977  int indep_rows_mod(int **mat0,int row,int col,int md,i
Line 2132  int indep_rows_mod(int **mat0,int row,int col,int md,i
       }        }
     for ( i = rank+1; i < row; i++ ) {      for ( i = rank+1; i < row; i++ ) {
       t = mat[i];        t = mat[i];
       if ( a = t[j] )        if ( ( a = t[j] ) != 0 )
         red_by_vect(md,t+j,pivot+j,md-a,col-j);          red_by_vect(md,t+j,pivot+j,md-a,col-j);
     }      }
     rank++;      rank++;
Line 2011  int generic_gauss_elim_sf(int **mat0,int row,int col,i
Line 2166  int generic_gauss_elim_sf(int **mat0,int row,int col,i
         *pk = _mulsf(*pk,inv);          *pk = _mulsf(*pk,inv);
     for ( i = rank+1; i < row; i++ ) {      for ( i = rank+1; i < row; i++ ) {
       t = mat[i];        t = mat[i];
       if ( a = t[j] )        if ( ( a = t[j] ) != 0 )
         red_by_vect_sf(md,t+j,pivot+j,_chsgnsf(a),col-j);          red_by_vect_sf(md,t+j,pivot+j,_chsgnsf(a),col-j);
     }      }
     rank++;      rank++;
Line 2021  int generic_gauss_elim_sf(int **mat0,int row,int col,i
Line 2176  int generic_gauss_elim_sf(int **mat0,int row,int col,i
       pivot = mat[l];        pivot = mat[l];
       for ( i = 0; i < l; i++ ) {        for ( i = 0; i < l; i++ ) {
         t = mat[i];          t = mat[i];
         if ( a = t[j] )          if ( ( a = t[j] ) != 0 )
           red_by_vect_sf(md,t+j,pivot+j,_chsgnsf(a),col-j);            red_by_vect_sf(md,t+j,pivot+j,_chsgnsf(a),col-j);
       }        }
       l--;        l--;
Line 2057  int lu_gfmmat(GFMMAT mat,unsigned int md,int *perm)
Line 2212  int lu_gfmmat(GFMMAT mat,unsigned int md,int *perm)
     pivot[k] = inv = invm(pivot[k],md);      pivot[k] = inv = invm(pivot[k],md);
     for ( i = k+1; i < row; i++ ) {      for ( i = k+1; i < row; i++ ) {
       t = a[i];        t = a[i];
       if ( m = t[k] ) {        if ( ( m = t[k] ) != 0 ) {
         DMAR(inv,m,0,md,t[k])          DMAR(inv,m,0,md,t[k])
         for ( j = k+1, m = md - t[k]; j < col; j++ )          for ( j = k+1, m = md - t[k]; j < col; j++ )
           if ( pivot[j] ) {            if ( pivot[j] ) {
Line 2113  int find_lhs_and_lu_mod(unsigned int **a,int row,int c
Line 2268  int find_lhs_and_lu_mod(unsigned int **a,int row,int c
     pivot[k] = inv = invm(pivot[k],md);      pivot[k] = inv = invm(pivot[k],md);
     for ( i = d+1; i < row; i++ ) {      for ( i = d+1; i < row; i++ ) {
       t = a[i];        t = a[i];
       if ( m = t[k] ) {        if ( ( m = t[k] ) != 0 ) {
         DMAR(inv,m,0,md,t[k])          DMAR(inv,m,0,md,t[k])
         for ( j = k+1, m = md - t[k]; j < col; j++ )          for ( j = k+1, m = md - t[k]; j < col; j++ )
           if ( pivot[j] ) {            if ( pivot[j] ) {
Line 2148  int lu_mod(unsigned int **a,int n,unsigned int md,int 
Line 2303  int lu_mod(unsigned int **a,int n,unsigned int md,int 
     inv = invm(pivot[k],md);      inv = invm(pivot[k],md);
     for ( i = k+1; i < n; i++ ) {      for ( i = k+1; i < n; i++ ) {
       t = a[i];        t = a[i];
       if ( m = t[k] ) {        if ( ( m = t[k] ) != 0 ) {
         DMAR(inv,m,0,md,t[k])          DMAR(inv,m,0,md,t[k])
         for ( j = k+1, m = md - t[k]; j < n; j++ )          for ( j = k+1, m = md - t[k]; j < n; j++ )
           if ( pivot[j] ) {            if ( pivot[j] ) {
Line 2177  void solve_by_lu_mod(int **a,int n,int md,int **b,int 
Line 2332  void solve_by_lu_mod(int **a,int n,int md,int **b,int 
   int i,j,k;    int i,j,k;
   unsigned int t,m,m2;    unsigned int t,m,m2;
   
   y = (int *)MALLOC_ATOMIC(n*sizeof(int));    y = (unsigned int *)MALLOC_ATOMIC(n*sizeof(int));
   c = (int *)MALLOC_ATOMIC(n*sizeof(int));    c = (unsigned int *)MALLOC_ATOMIC(n*sizeof(int));
   m2 = md>>1;    m2 = md>>1;
   for ( k = 0; k < l; k++ ) {    for ( k = 0; k < l; k++ ) {
     /* copy b[.][k] to c */      /* copy b[.][k] to c */
Line 2338  int gauss_elim_geninv_mod(unsigned int **mat,int row,i
Line 2493  int gauss_elim_geninv_mod(unsigned int **mat,int row,i
       pivot[k] = dmar(pivot[k],inv,0,md);        pivot[k] = dmar(pivot[k],inv,0,md);
     for ( i = j+1; i < row; i++ ) {      for ( i = j+1; i < row; i++ ) {
       t = mat[i];        t = mat[i];
       if ( a = t[j] )        if ( ( a = t[j] ) != 0  )
         for ( k = j, a = md - a; k < m; k++ )          for ( k = j, a = md - a; k < m; k++ )
           t[k] = dmar(pivot[k],a,t[k],md);            t[k] = dmar(pivot[k],a,t[k],md);
     }      }
Line 2347  int gauss_elim_geninv_mod(unsigned int **mat,int row,i
Line 2502  int gauss_elim_geninv_mod(unsigned int **mat,int row,i
     pivot = mat[j];      pivot = mat[j];
     for ( i = j-1; i >= 0; i-- ) {      for ( i = j-1; i >= 0; i-- ) {
       t = mat[i];        t = mat[i];
       if ( a = t[j] )        if ( ( a = t[j] ) != 0 )
         for ( k = j, a = md - a; k < m; k++ )          for ( k = j, a = md - a; k < m; k++ )
           t[k] = dmar(pivot[k],a,t[k],md);            t[k] = dmar(pivot[k],a,t[k],md);
     }      }
Line 2565  int gauss_elim_geninv_mod_swap(unsigned int **mat,int 
Line 2720  int gauss_elim_geninv_mod_swap(unsigned int **mat,int 
         pivot[k] = (unsigned int)dmar(pivot[k],inv,0,md);          pivot[k] = (unsigned int)dmar(pivot[k],inv,0,md);
     for ( i = j+1; i < row; i++ ) {      for ( i = j+1; i < row; i++ ) {
       t = mat[i];        t = mat[i];
       if ( a = t[j] )        if ( ( a = t[j] ) != 0 )
         for ( k = j, a = md - a; k < m; k++ )          for ( k = j, a = md - a; k < m; k++ )
           if ( pivot[k] )            if ( pivot[k] )
             t[k] = dmar(pivot[k],a,t[k],md);              t[k] = dmar(pivot[k],a,t[k],md);
Line 2575  int gauss_elim_geninv_mod_swap(unsigned int **mat,int 
Line 2730  int gauss_elim_geninv_mod_swap(unsigned int **mat,int 
     pivot = mat[j];      pivot = mat[j];
     for ( i = j-1; i >= 0; i-- ) {      for ( i = j-1; i >= 0; i-- ) {
       t = mat[i];        t = mat[i];
       if ( a = t[j] )        if ( ( a = t[j] ) != 0 )
         for ( k = j, a = md - a; k < m; k++ )          for ( k = j, a = md - a; k < m; k++ )
           if ( pivot[k] )            if ( pivot[k] )
             t[k] = dmar(pivot[k],a,t[k],md);              t[k] = dmar(pivot[k],a,t[k],md);
Line 2611  void Pgeninv_sf_swap(NODE arg,LIST *rp)
Line 2766  void Pgeninv_sf_swap(NODE arg,LIST *rp)
   for ( i = 0; i < row; i++ ) {    for ( i = 0; i < row; i++ ) {
     bzero((char *)wmat[i],(col+row)*sizeof(int));      bzero((char *)wmat[i],(col+row)*sizeof(int));
     for ( j = 0; j < col; j++ )      for ( j = 0; j < col; j++ )
       if ( q = (GFS)mat[i][j] )        if ( ( q = (GFS)mat[i][j] ) != 0 )
         wmat[i][j] = FTOIF(CONT(q));          wmat[i][j] = FTOIF(CONT(q));
     wmat[i][col+i] = _onesf();      wmat[i][col+i] = _onesf();
   }    }
Line 2622  void Pgeninv_sf_swap(NODE arg,LIST *rp)
Line 2777  void Pgeninv_sf_swap(NODE arg,LIST *rp)
     MKMAT(mat1,col,col);      MKMAT(mat1,col,col);
     for ( i = 0, tmat = (GFS **)mat1->body; i < col; i++ )      for ( i = 0, tmat = (GFS **)mat1->body; i < col; i++ )
       for ( j = 0; j < col; j++ )        for ( j = 0; j < col; j++ )
         if ( t = invmat[i][j] ) {          if ( ( t = invmat[i][j] ) != 0 ) {
           MKGFS(IFTOF(t),tmat[i][j]);            MKGFS(IFTOF(t),tmat[i][j]);
         }          }
     MKVECT(vect1,row);      MKVECT(vect1,row);
Line 2659  int gauss_elim_geninv_sf_swap(int **mat,int row,int co
Line 2814  int gauss_elim_geninv_sf_swap(int **mat,int row,int co
         pivot[k] = _mulsf(pivot[k],inv);          pivot[k] = _mulsf(pivot[k],inv);
     for ( i = j+1; i < row; i++ ) {      for ( i = j+1; i < row; i++ ) {
       t = mat[i];        t = mat[i];
       if ( a = t[j] )        if ( ( a = t[j] ) != 0 )
         for ( k = j, a = _chsgnsf(a); k < m; k++ )          for ( k = j, a = _chsgnsf(a); k < m; k++ )
           if ( pivot[k] ) {            if ( pivot[k] ) {
             u = _mulsf(pivot[k],a);              u = _mulsf(pivot[k],a);
Line 2671  int gauss_elim_geninv_sf_swap(int **mat,int row,int co
Line 2826  int gauss_elim_geninv_sf_swap(int **mat,int row,int co
     pivot = mat[j];      pivot = mat[j];
     for ( i = j-1; i >= 0; i-- ) {      for ( i = j-1; i >= 0; i-- ) {
       t = mat[i];        t = mat[i];
       if ( a = t[j] )        if ( ( a = t[j] ) != 0 )
         for ( k = j, a = _chsgnsf(a); k < m; k++ )          for ( k = j, a = _chsgnsf(a); k < m; k++ )
           if ( pivot[k] ) {            if ( pivot[k] ) {
             u = _mulsf(pivot[k],a);              u = _mulsf(pivot[k],a);
Line 2844  int generate_ONB_polynomial(UP2 *rp,int m,int type)
Line 2999  int generate_ONB_polynomial(UP2 *rp,int m,int type)
       for ( i = 0; i < w; i++ )        for ( i = 0; i < w; i++ )
         f->b[i] = 0xffffffff;          f->b[i] = 0xffffffff;
       /* mask the top word if necessary */        /* mask the top word if necessary */
       if ( r = (m+1)&31 )        if ( ( r = (m+1)&31 ) != 0 )
         f->b[w-1] &= (1<<r)-1;          f->b[w-1] &= (1<<r)-1;
       return 0;        return 0;
       break;        break;
Legend:
Removed from v.1.3  
changed lines
  Added in v.1.11
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>