[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.1 and 1.7

version 1.1, 2018/09/19 05:45:05 version 1.7, 2020/01/09 01:47:40
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: OpenXM_contrib2/asir2018/builtin/array.c,v 1.6 2019/12/13 14:40:49 fujimoto 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 206  void Plusolve_main(NODE arg,VECT *rp)
Line 213  void Plusolve_main(NODE arg,VECT *rp)
   v = (VECT)ARG0(arg); len = v->len;    v = (VECT)ARG0(arg); len = v->len;
   d = (Z *)BDY(v);    d = (Z *)BDY(v);
   rhs = (int *)MALLOC_ATOMIC(len*sizeof(int));    rhs = (int *)MALLOC_ATOMIC(len*sizeof(int));
   for ( i = 0; i < len; i++ ) rhs[i] = QTOS(d[i]);    for ( i = 0; i < len; i++ ) rhs[i] = ZTOS(d[i]);
   solve_l(ll,l,len,rhs,modulus);    solve_l(ll,l,len,rhs,modulus);
   solve_u(ul,u,len,rhs,modulus);    solve_u(ul,u,len,rhs,modulus);
   NEWVECT(r); r->len = len;    NEWVECT(r); r->len = len;
   r->body = (pointer *)MALLOC(len*sizeof(pointer));    r->body = (pointer *)MALLOC(len*sizeof(pointer));
   p = (Z *)r->body;    p = (Z *)r->body;
   for ( i = 0; i < len; i++ )    for ( i = 0; i < len; i++ )
     STOQ(rhs[i],p[i]);      STOZ(rhs[i],p[i]);
   *rp = r;    *rp = r;
 }  }
   
Line 565  void Psepmat_destructive(NODE arg,LIST *rp)
Line 572  void Psepmat_destructive(NODE arg,LIST *rp)
   
 void Psepvect(NODE arg,VECT *rp)  void Psepvect(NODE arg,VECT *rp)
 {  {
   sepvect((VECT)ARG0(arg),QTOS((Q)ARG1(arg)),rp);    sepvect((VECT)ARG0(arg),ZTOS((Q)ARG1(arg)),rp);
 }  }
   
 void sepvect(VECT v,int d,VECT *rp)  void sepvect(VECT v,int d,VECT *rp)
Line 601  void Pnewvect(NODE arg,VECT *rp)
Line 608  void Pnewvect(NODE arg,VECT *rp)
   NODE tn;    NODE tn;
   
   asir_assert(ARG0(arg),O_N,"newvect");    asir_assert(ARG0(arg),O_N,"newvect");
   len = QTOS((Q)ARG0(arg));    len = ZTOS((Q)ARG0(arg));
   if ( len < 0 )    if ( len < 0 )
     error("newvect : invalid size");      error("newvect : invalid size");
   MKVECT(vect,len);    MKVECT(vect,len);
Line 691  void Pnewbytearray(NODE arg,BYTEARRAY *rp)
Line 698  void Pnewbytearray(NODE arg,BYTEARRAY *rp)
       case O_N:        case O_N:
         if ( !RATN(ARG0(arg)) )          if ( !RATN(ARG0(arg)) )
           error("newbytearray : invalid argument");            error("newbytearray : invalid argument");
         len = QTOS((Q)ARG0(arg));          len = ZTOS((Q)ARG0(arg));
         if ( len < 0 )          if ( len < 0 )
           error("newbytearray : invalid size");            error("newbytearray : invalid size");
         MKBYTEARRAY(array,len);          MKBYTEARRAY(array,len);
Line 701  void Pnewbytearray(NODE arg,BYTEARRAY *rp)
Line 708  void Pnewbytearray(NODE arg,BYTEARRAY *rp)
     }      }
   } else if ( ac == 2 ) {    } else if ( ac == 2 ) {
     asir_assert(ARG0(arg),O_N,"newbytearray");      asir_assert(ARG0(arg),O_N,"newbytearray");
     len = QTOS((Q)ARG0(arg));      len = ZTOS((Q)ARG0(arg));
     if ( len < 0 )      if ( len < 0 )
       error("newbytearray : invalid size");        error("newbytearray : invalid size");
     MKBYTEARRAY(array,len);      MKBYTEARRAY(array,len);
Line 715  void Pnewbytearray(NODE arg,BYTEARRAY *rp)
Line 722  void Pnewbytearray(NODE arg,BYTEARRAY *rp)
         if ( r <= len ) {          if ( r <= len ) {
           for ( i = 0, tn = BDY(list), vb = BDY(array); tn;            for ( i = 0, tn = BDY(list), vb = BDY(array); tn;
             i++, tn = NEXT(tn) )              i++, tn = NEXT(tn) )
             vb[i] = (unsigned char)QTOS((Q)BDY(tn));              vb[i] = (unsigned char)ZTOS((Q)BDY(tn));
         }          }
         break;          break;
       case O_STR:        case O_STR:
Line 746  void Pmemoryplot_to_coord(NODE arg,LIST *rp)
Line 753  void Pmemoryplot_to_coord(NODE arg,LIST *rp)
   
   asir_assert(ARG0(arg),O_LIST,"memoryplot_to_coord");    asir_assert(ARG0(arg),O_LIST,"memoryplot_to_coord");
   arg = BDY((LIST)ARG0(arg));    arg = BDY((LIST)ARG0(arg));
   len = QTOS((Q)ARG0(arg));    len = ZTOS((Q)ARG0(arg));
   blen = (len+7)/8;    blen = (len+7)/8;
   y = QTOS((Q)ARG1(arg));    y = ZTOS((Q)ARG1(arg));
   ba = (BYTEARRAY)ARG2(arg); a = ba->body;    ba = (BYTEARRAY)ARG2(arg); a = ba->body;
   r0 = 0;    r0 = 0;
   for ( j = 0; j < y; j++ )    for ( j = 0; j < y; j++ )
     for ( i = 0; i < len; i++ )      for ( i = 0; i < len; i++ )
       if ( MEMORY_GETPOINT(a,blen,i,j) ) {        if ( MEMORY_GETPOINT(a,blen,i,j) ) {
         NEXTNODE(r0,r);          NEXTNODE(r0,r);
         STOQ(i,iq); STOQ(j,jq);          STOZ(i,iq); STOZ(j,jq);
         n = mknode(2,iq,jq);          n = mknode(2,iq,jq);
         MKLIST(l,n);          MKLIST(l,n);
         BDY(r) = l;          BDY(r) = l;
Line 775  void Pnewmat(NODE arg,MAT *rp)
Line 782  void Pnewmat(NODE arg,MAT *rp)
   
   asir_assert(ARG0(arg),O_N,"newmat");    asir_assert(ARG0(arg),O_N,"newmat");
   asir_assert(ARG1(arg),O_N,"newmat");    asir_assert(ARG1(arg),O_N,"newmat");
   row = QTOS((Q)ARG0(arg)); col = QTOS((Q)ARG1(arg));    row = ZTOS((Q)ARG0(arg)); col = ZTOS((Q)ARG1(arg));
   if ( row < 0 || col < 0 )    if ( row < 0 || col < 0 )
     error("newmat : invalid size");      error("newmat : invalid size");
   MKMAT(m,row,col);    MKMAT(m,row,col);
Line 967  void Premainder(NODE arg,Obj *rp)
Line 974  void Premainder(NODE arg,Obj *rp)
         cmp(md,(P)a,(P *)rp); break;          cmp(md,(P)a,(P *)rp); break;
       case O_VECT:        case O_VECT:
         /* strage spec */          /* strage spec */
         smd = QTOS(md);          smd = ZTOS(md);
         v = (VECT)a; n = v->len; vb = v->body;          v = (VECT)a; n = v->len; vb = v->body;
         MKVECT(w,n); wb = w->body;          MKVECT(w,n); wb = w->body;
         for ( i = 0; i < n; i++ ) {          for ( i = 0; i < n; i++ ) {
Line 1015  void Psremainder(NODE arg,Obj *rp)
Line 1022  void Psremainder(NODE arg,Obj *rp)
       case O_P:        case O_P:
         cmp(md,(P)a,(P *)rp); break;          cmp(md,(P)a,(P *)rp); break;
       case O_VECT:        case O_VECT:
         smd = QTOS(md);          smd = ZTOS(md);
         v = (VECT)a; n = v->len; vb = v->body;          v = (VECT)a; n = v->len; vb = v->body;
         MKVECT(w,n); wb = w->body;          MKVECT(w,n); wb = w->body;
         for ( i = 0; i < n; i++ )          for ( i = 0; i < n; i++ )
Line 1049  void Psize(NODE arg,LIST *rp)
Line 1056  void Psize(NODE arg,LIST *rp)
     switch (OID(ARG0(arg))) {      switch (OID(ARG0(arg))) {
       case O_VECT:        case O_VECT:
         n = ((VECT)ARG0(arg))->len;          n = ((VECT)ARG0(arg))->len;
         STOQ(n,q); MKNODE(t,q,0);          STOZ(n,q); MKNODE(t,q,0);
         break;          break;
       case O_MAT:        case O_MAT:
         n = ((MAT)ARG0(arg))->row; m = ((MAT)ARG0(arg))->col;          n = ((MAT)ARG0(arg))->row; m = ((MAT)ARG0(arg))->col;
         STOQ(m,q); MKNODE(s,q,0); STOQ(n,q); MKNODE(t,q,s);          STOZ(m,q); MKNODE(s,q,0); STOZ(n,q); MKNODE(t,q,s);
         break;          break;
       case O_IMAT:        case O_IMAT:
         n = ((IMAT)ARG0(arg))->row; m = ((IMAT)ARG0(arg))->col;          n = ((IMAT)ARG0(arg))->row; m = ((IMAT)ARG0(arg))->col;
         STOQ(m,q); MKNODE(s,q,0); STOQ(n,q); MKNODE(t,q,s);          STOZ(m,q); MKNODE(s,q,0); STOZ(n,q); MKNODE(t,q,s);
         break;          break;
       default:        default:
         error("size : invalid argument"); break;          error("size : invalid argument"); break;
Line 1080  void Pdet(NODE arg,P *rp)
Line 1087  void Pdet(NODE arg,P *rp)
   else if ( argc(arg) == 1 )    else if ( argc(arg) == 1 )
     detp(CO,(P **)BDY(m),m->row,rp);      detp(CO,(P **)BDY(m),m->row,rp);
   else {    else {
     n = m->row; mod = QTOS((Q)ARG1(arg)); mat = (P **)BDY(m);      n = m->row; mod = ZTOS((Q)ARG1(arg)); mat = (P **)BDY(m);
     w = (P **)almat_pointer(n,n);      w = (P **)almat_pointer(n,n);
     for ( i = 0; i < n; i++ )      for ( i = 0; i < n; i++ )
       for ( j = 0; j < n; j++ )        for ( j = 0; j < n; j++ )
Line 1109  void Pinvmat(NODE arg,LIST *rp)
Line 1116  void Pinvmat(NODE arg,LIST *rp)
     nd = mknode(2,r,dn);      nd = mknode(2,r,dn);
     MKLIST(*rp,nd);      MKLIST(*rp,nd);
   } else {    } else {
     n = m->row; mod = QTOS((Q)ARG1(arg)); mat = (P **)BDY(m);      n = m->row; mod = ZTOS((Q)ARG1(arg)); mat = (P **)BDY(m);
     w = (P **)almat_pointer(n,n);      w = (P **)almat_pointer(n,n);
     for ( i = 0; i < n; i++ )      for ( i = 0; i < n; i++ )
       for ( j = 0; j < n; j++ )        for ( j = 0; j < n; j++ )
Line 1169  void Pgeneric_gauss_elim(NODE arg,LIST *rp)
Line 1176  void Pgeneric_gauss_elim(NODE arg,LIST *rp)
   MKVECT(rind,rank);    MKVECT(rind,rank);
   MKVECT(cind,t);    MKVECT(cind,t);
   for ( i = 0; i < rank; i++ ) {    for ( i = 0; i < rank; i++ ) {
     STOQ(ri[i],q);      STOZ(ri[i],q);
     BDY(rind)[i] = (pointer)q;      BDY(rind)[i] = (pointer)q;
   }    }
   for ( i = 0; i < t; i++ ) {    for ( i = 0; i < t; i++ ) {
     STOQ(ci[i],q);      STOZ(ci[i],q);
     BDY(cind)[i] = (pointer)q;      BDY(cind)[i] = (pointer)q;
   }    }
   n0 = mknode(4,nm,dn,rind,cind);    n0 = mknode(4,nm,dn,rind,cind);
Line 1195  void Pindep_rows_mod(NODE arg,VECT *rp)
Line 1202  void Pindep_rows_mod(NODE arg,VECT *rp)
   
   asir_assert(ARG0(arg),O_MAT,"indep_rows_mod");    asir_assert(ARG0(arg),O_MAT,"indep_rows_mod");
   asir_assert(ARG1(arg),O_N,"indep_rows_mod");    asir_assert(ARG1(arg),O_N,"indep_rows_mod");
   m = (MAT)ARG0(arg); md = QTOS((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 = (int **)almat(row,col);    wmat = (int **)almat(row,col);
   
Line 1211  void Pindep_rows_mod(NODE arg,VECT *rp)
Line 1218  void Pindep_rows_mod(NODE arg,VECT *rp)
   MKVECT(rind,rank);    MKVECT(rind,rank);
   rib = (Z *)rind->body;    rib = (Z *)rind->body;
   for ( j = 0; j < rank; j++ ) {    for ( j = 0; j < rank; j++ ) {
     STOQ(rowstat[j],rib[j]);      STOZ(rowstat[j],rib[j]);
   }    }
     *rp = rind;      *rp = rind;
 }  }
Line 1229  void Pindep_rows_mod(NODE arg,VECT *rp)
Line 1236  void Pindep_rows_mod(NODE arg,VECT *rp)
     B[I] <-> x_{R[I]}+B[I][0]x_{C[0]}+B[I][1]x_{C[1]}+...      B[I] <-> x_{R[I]}+B[I][0]x_{C[0]}+B[I][1]x_{C[1]}+...
 */  */
   
   #if SIZEOF_LONG==8
   void Pgeneric_gauss_elim_mod64(NODE arg,LIST *rp)
   {
     NODE n0;
     MAT m,mat;
     VECT rind,cind,rnum;
     Z **tmat;
     mp_limb_t **wmat,**row0;
     Z *rib,*cib,*rnb;
     int *colstat;
     mp_limb_t *p;
     Z q;
     mp_limb_t md;
     int i,j,k,l,row,col,t,rank;
   
     asir_assert(ARG0(arg),O_MAT,"generic_gauss_elim_mod64");
     asir_assert(ARG1(arg),O_N,"generic_gauss_elim_mod64");
     m = (MAT)ARG0(arg); md = ZTOS((Z)ARG1(arg));
     row = m->row; col = m->col; tmat = (Z **)m->body;
     wmat = (mp_limb_t **)almat64(row,col);
   
     row0 = (mp_limb_t **)ALLOCA(row*sizeof(mp_limb_t *));
     for ( i = 0; i < row; i++ ) row0[i] = wmat[i];
   
     colstat = (int *)MALLOC_ATOMIC(col*sizeof(int));
     for ( i = 0; i < row; i++ )
       for ( j = 0; j < col; j++ )
         wmat[i][j] = remqi64((Q)tmat[i][j],md);
     rank = generic_gauss_elim_mod64(wmat,row,col,md,colstat);
   
     MKVECT(rnum,rank);
     rnb = (Z *)rnum->body;
     for ( i = 0; i < rank; i++ )
       for ( j = 0, p = wmat[i]; j < row; j++ )
         if ( p == row0[j] )
           STOZ(j,rnb[i]);
   
     MKMAT(mat,rank,col-rank);
     tmat = (Z **)mat->body;
     for ( i = 0; i < rank; i++ )
       for ( j = k = 0; j < col; j++ )
         if ( !colstat[j] ) {
           UTOZ(wmat[i][j],tmat[i][k]); k++;
         }
   
     MKVECT(rind,rank);
     MKVECT(cind,col-rank);
     rib = (Z *)rind->body; cib = (Z *)cind->body;
     for ( j = k = l = 0; j < col; j++ )
       if ( colstat[j] ) {
         STOZ(j,rib[k]); k++;
       } else {
         STOZ(j,cib[l]); l++;
       }
     n0 = mknode(4,mat,rind,cind,rnum);
     MKLIST(*rp,n0);
   }
   #endif
   
 void Pgeneric_gauss_elim_mod(NODE arg,LIST *rp)  void Pgeneric_gauss_elim_mod(NODE arg,LIST *rp)
 {  {
   NODE n0;    NODE n0;
Line 1239  void Pgeneric_gauss_elim_mod(NODE arg,LIST *rp)
Line 1305  void Pgeneric_gauss_elim_mod(NODE arg,LIST *rp)
   Z *rib,*cib,*rnb;    Z *rib,*cib,*rnb;
   int *colstat,*p;    int *colstat,*p;
   Z q;    Z q;
     long mdl;
   int md,i,j,k,l,row,col,t,rank;    int md,i,j,k,l,row,col,t,rank;
   
   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");
   m = (MAT)ARG0(arg); md = QTOS((Z)ARG1(arg));  #if SIZEOF_LONG==8
     mdl = ZTOS((Z)ARG1(arg));
     if ( mdl >= ((mp_limb_t)1)<<32 ) {
       Pgeneric_gauss_elim_mod64(arg,rp);
       return;
     }
   #endif
     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 = (int **)almat(row,col);    wmat = (int **)almat(row,col);
   
Line 1261  void Pgeneric_gauss_elim_mod(NODE arg,LIST *rp)
Line 1336  void Pgeneric_gauss_elim_mod(NODE arg,LIST *rp)
   for ( i = 0; i < rank; i++ )    for ( i = 0; i < rank; i++ )
     for ( j = 0, p = wmat[i]; j < row; j++ )      for ( j = 0, p = wmat[i]; j < row; j++ )
       if ( p == row0[j] )        if ( p == row0[j] )
         STOQ(j,rnb[i]);          STOZ(j,rnb[i]);
   
   MKMAT(mat,rank,col-rank);    MKMAT(mat,rank,col-rank);
   tmat = (Z **)mat->body;    tmat = (Z **)mat->body;
   for ( i = 0; i < rank; i++ )    for ( i = 0; i < rank; i++ )
     for ( j = k = 0; j < col; j++ )      for ( j = k = 0; j < col; j++ )
       if ( !colstat[j] ) {        if ( !colstat[j] ) {
         UTOQ(wmat[i][j],tmat[i][k]); k++;          UTOZ(wmat[i][j],tmat[i][k]); k++;
       }        }
   
   MKVECT(rind,rank);    MKVECT(rind,rank);
Line 1276  void Pgeneric_gauss_elim_mod(NODE arg,LIST *rp)
Line 1351  void Pgeneric_gauss_elim_mod(NODE arg,LIST *rp)
   rib = (Z *)rind->body; cib = (Z *)cind->body;    rib = (Z *)rind->body; cib = (Z *)cind->body;
   for ( j = k = l = 0; j < col; j++ )    for ( j = k = l = 0; j < col; j++ )
     if ( colstat[j] ) {      if ( colstat[j] ) {
       STOQ(j,rib[k]); k++;        STOZ(j,rib[k]); k++;
     } else {      } else {
       STOQ(j,cib[l]); l++;        STOZ(j,cib[l]); l++;
     }      }
   n0 = mknode(4,mat,rind,cind,rnum);    n0 = mknode(4,mat,rind,cind,rnum);
   MKLIST(*rp,n0);    MKLIST(*rp,n0);
 }  }
   
   
 void Pleqm(NODE arg,VECT *rp)  void Pleqm(NODE arg,VECT *rp)
 {  {
   MAT m;    MAT m;
Line 1296  void Pleqm(NODE arg,VECT *rp)
Line 1372  void Pleqm(NODE arg,VECT *rp)
   
   asir_assert(ARG0(arg),O_MAT,"leqm");    asir_assert(ARG0(arg),O_MAT,"leqm");
   asir_assert(ARG1(arg),O_N,"leqm");    asir_assert(ARG1(arg),O_N,"leqm");
   m = (MAT)ARG0(arg); md = QTOS((Z)ARG1(arg));    m = (MAT)ARG0(arg); md = ZTOS((Z)ARG1(arg));
   row = m->row; col = m->col; mat = m->body;    row = m->row; col = m->col; mat = m->body;
   wmat = (int **)almat(row,col);    wmat = (int **)almat(row,col);
   for ( i = 0; i < row; i++ )    for ( i = 0; i < row; i++ )
Line 1311  void Pleqm(NODE arg,VECT *rp)
Line 1387  void Pleqm(NODE arg,VECT *rp)
     n = col - 1;      n = col - 1;
     MKVECT(vect,n);      MKVECT(vect,n);
     for ( i = 0, v = (Z *)vect->body; i < n; i++ ) {      for ( i = 0, v = (Z *)vect->body; i < n; i++ ) {
       t = (md-wmat[i][n])%md; STOQ(t,v[i]);        t = (md-wmat[i][n])%md; STOZ(t,v[i]);
     }      }
     *rp = vect;      *rp = vect;
   }    }
Line 1354  int gauss_elim_mod(int **mat,int row,int col,int md)
Line 1430  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 1635  void red_by_vect(int m,unsigned int *p,unsigned int *r
Line 1711  void red_by_vect(int m,unsigned int *p,unsigned int *r
     }      }
 }  }
   
 #if defined(__GNUC__) && SIZEOF_LONG==8  #if SIZEOF_LONG==8
 /* 64bit vector += UNIT vector(normalized) */  /* mp_limb_t vector += U32 vector(normalized)*U32 */
   
 void red_by_vect64(int m, U64 *p,unsigned int *c,U64 *r,unsigned int hc,int len)  void red_by_vect64(int m, mp_limb_t *p,unsigned int *c,mp_limb_t *r,unsigned int hc,int len)
 {  {
   U64 t;    mp_limb_t t;
   
   /* (p[0],c[0]) is normalized */    /* (p[0],c[0]) is normalized */
   *p++ = 0; *c++ = 0; r++; len--;    *p++ = 0; *c++ = 0; r++; len--;
Line 1651  void red_by_vect64(int m, U64 *p,unsigned int *c,U64 *
Line 1727  void red_by_vect64(int m, U64 *p,unsigned int *c,U64 *
       *p = t;        *p = t;
     }      }
 }  }
   
   /* mp_limb_t vector = (mp_limb_t vector+mp_limb_t vector*mp_limb_t)%mp_limb_t */
   
   void red_by_vect64mod(mp_limb_t m, mp_limb_t *p,mp_limb_t *r,mp_limb_t hc,int len)
   {
     *p++ = 0; r++; len--;
     for ( ; len; len--, r++, p++ )
       if ( *r )
         *p = muladdmod64(*r,hc,*p,m);
   }
   
   int generic_gauss_elim_mod64(mp_limb_t **mat,int row,int col,mp_limb_t md,int *colstat)
   {
     int i,j,k,l,rank;
     mp_limb_t inv,a;
     mp_limb_t *t,*pivot,*pk;
   
     for ( rank = 0, j = 0; j < col; j++ ) {
       for ( i = rank; i < row; i++ )
         if ( mat[i][j] )
           break;
       if ( i == row ) {
         colstat[j] = 0;
         continue;
       } else
         colstat[j] = 1;
       if ( i != rank ) {
         t = mat[i]; mat[i] = mat[rank]; mat[rank] = t;
       }
       pivot = mat[rank];
       inv = invmod64(pivot[j],md);
       for ( k = j, pk = pivot+k; k < col; k++, pk++ )
         if ( *pk )
           *pk = mulmod64(*pk,inv,md);
       for ( i = rank+1; i < row; i++ ) {
         t = mat[i];
         if ( a = t[j] )
           red_by_vect64mod(md,t+j,pivot+j,md-a,col-j);
       }
       rank++;
     }
     for ( j = col-1, l = rank-1; j >= 0; j-- )
       if ( colstat[j] ) {
         pivot = mat[l];
         for ( i = 0; i < l; i++ ) {
           t = mat[i];
           if ( a = t[j] )
             red_by_vect64mod(md,t+j,pivot+j,md-a,col-j);
         }
         l--;
       }
     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 2102  void Pleqm1(NODE arg,VECT *rp)
Line 2354  void Pleqm1(NODE arg,VECT *rp)
   
   asir_assert(ARG0(arg),O_MAT,"leqm1");    asir_assert(ARG0(arg),O_MAT,"leqm1");
   asir_assert(ARG1(arg),O_N,"leqm1");    asir_assert(ARG1(arg),O_N,"leqm1");
   m = (MAT)ARG0(arg); md = QTOS((Z)ARG1(arg));    m = (MAT)ARG0(arg); md = ZTOS((Z)ARG1(arg));
   row = m->row; col = m->col; mat = m->body;    row = m->row; col = m->col; mat = m->body;
   wmat = (int **)almat(row,col);    wmat = (int **)almat(row,col);
   for ( i = 0; i < row; i++ )    for ( i = 0; i < row; i++ )
Line 2117  void Pleqm1(NODE arg,VECT *rp)
Line 2369  void Pleqm1(NODE arg,VECT *rp)
     n = col - 1;      n = col - 1;
     MKVECT(vect,n);      MKVECT(vect,n);
     for ( i = 0, v = (Z *)vect->body; i < n; i++ ) {      for ( i = 0, v = (Z *)vect->body; i < n; i++ ) {
       t = (md-wmat[i][n])%md; STOQ(t,v[i]);        t = (md-wmat[i][n])%md; STOZ(t,v[i]);
     }      }
     *rp = vect;      *rp = vect;
   }    }
Line 2173  void Pgeninvm(NODE arg,LIST *rp)
Line 2425  void Pgeninvm(NODE arg,LIST *rp)
   
   asir_assert(ARG0(arg),O_MAT,"leqm1");    asir_assert(ARG0(arg),O_MAT,"leqm1");
   asir_assert(ARG1(arg),O_N,"leqm1");    asir_assert(ARG1(arg),O_N,"leqm1");
   m = (MAT)ARG0(arg); md = QTOS((Q)ARG1(arg));    m = (MAT)ARG0(arg); md = ZTOS((Q)ARG1(arg));
   row = m->row; col = m->col; mat = m->body;    row = m->row; col = m->col; mat = m->body;
   wmat = (unsigned int **)almat(row,col+row);    wmat = (unsigned int **)almat(row,col+row);
   for ( i = 0; i < row; i++ ) {    for ( i = 0; i < row; i++ ) {
Line 2188  void Pgeninvm(NODE arg,LIST *rp)
Line 2440  void Pgeninvm(NODE arg,LIST *rp)
     MKMAT(mat1,col,row); MKMAT(mat2,row-col,row);      MKMAT(mat1,col,row); MKMAT(mat2,row-col,row);
     for ( i = 0, tmat = (Z **)mat1->body; i < col; i++ )      for ( i = 0, tmat = (Z **)mat1->body; i < col; i++ )
       for ( j = 0; j < row; j++ )        for ( j = 0; j < row; j++ )
         UTOQ(wmat[i][j+col],tmat[i][j]);          UTOZ(wmat[i][j+col],tmat[i][j]);
     for ( tmat = (Z **)mat2->body; i < row; i++ )      for ( tmat = (Z **)mat2->body; i < row; i++ )
       for ( j = 0; j < row; j++ )        for ( j = 0; j < row; j++ )
         UTOQ(wmat[i][j+col],tmat[i-col][j]);          UTOZ(wmat[i][j+col],tmat[i-col][j]);
      MKNODE(node2,mat2,0); MKNODE(node1,mat1,node2); MKLIST(*rp,node1);       MKNODE(node2,mat2,0); MKNODE(node1,mat1,node2); MKLIST(*rp,node1);
   }    }
 }  }
Line 2244  void Psolve_by_lu_gfmmat(NODE arg,VECT *rp)
Line 2496  void Psolve_by_lu_gfmmat(NODE arg,VECT *rp)
   lu = (GFMMAT)ARG0(arg);    lu = (GFMMAT)ARG0(arg);
   perm = (Z *)BDY((VECT)ARG1(arg));    perm = (Z *)BDY((VECT)ARG1(arg));
   rhs = (Z *)BDY((VECT)ARG2(arg));    rhs = (Z *)BDY((VECT)ARG2(arg));
   md = (unsigned int)QTOS((Z)ARG3(arg));    md = (unsigned int)ZTOS((Z)ARG3(arg));
   n = lu->col;    n = lu->col;
   b = (unsigned int *)MALLOC_ATOMIC(n*sizeof(int));    b = (unsigned int *)MALLOC_ATOMIC(n*sizeof(int));
   sol = (unsigned int *)MALLOC_ATOMIC(n*sizeof(int));    sol = (unsigned int *)MALLOC_ATOMIC(n*sizeof(int));
   for ( i = 0; i < n; i++ )    for ( i = 0; i < n; i++ )
     b[i] = QTOS(rhs[QTOS(perm[i])]);      b[i] = ZTOS(rhs[ZTOS(perm[i])]);
   solve_by_lu_gfmmat(lu,md,b,sol);    solve_by_lu_gfmmat(lu,md,b,sol);
   MKVECT(r,n);    MKVECT(r,n);
   for ( i = 0, v = (Z *)r->body; i < n; i++ )    for ( i = 0, v = (Z *)r->body; i < n; i++ )
       UTOQ(sol[i],v[i]);        UTOZ(sol[i],v[i]);
   *rp = r;    *rp = r;
 }  }
   
Line 2308  void Plu_mat(NODE arg,LIST *rp)
Line 2560  void Plu_mat(NODE arg,LIST *rp)
   lu_dec_cr(m,lu,&dn,&iperm);    lu_dec_cr(m,lu,&dn,&iperm);
   MKVECT(perm,n);    MKVECT(perm,n);
   for ( i = 0, v = (Q *)perm->body; i < n; i++ )    for ( i = 0, v = (Q *)perm->body; i < n; i++ )
     STOQ(iperm[i],v[i]);      STOZ(iperm[i],v[i]);
   n0 = mknode(3,lu,dn,perm);    n0 = mknode(3,lu,dn,perm);
   MKLIST(*rp,n0);    MKLIST(*rp,n0);
 }  }
Line 2327  void Plu_gfmmat(NODE arg,LIST *rp)
Line 2579  void Plu_gfmmat(NODE arg,LIST *rp)
   
   asir_assert(ARG0(arg),O_MAT,"lu_gfmmat");    asir_assert(ARG0(arg),O_MAT,"lu_gfmmat");
   asir_assert(ARG1(arg),O_N,"lu_gfmmat");    asir_assert(ARG1(arg),O_N,"lu_gfmmat");
   m = (MAT)ARG0(arg); md = (unsigned int)QTOS((Q)ARG1(arg));    m = (MAT)ARG0(arg); md = (unsigned int)ZTOS((Q)ARG1(arg));
   mat_to_gfmmat(m,md,&mm);    mat_to_gfmmat(m,md,&mm);
   row = m->row;    row = m->row;
   col = m->col;    col = m->col;
Line 2338  void Plu_gfmmat(NODE arg,LIST *rp)
Line 2590  void Plu_gfmmat(NODE arg,LIST *rp)
   else {    else {
     MKVECT(perm,row);      MKVECT(perm,row);
     for ( i = 0, v = (Z *)perm->body; i < row; i++ )      for ( i = 0, v = (Z *)perm->body; i < row; i++ )
       STOQ(iperm[i],v[i]);        STOZ(iperm[i],v[i]);
     n0 = mknode(2,mm,perm);      n0 = mknode(2,mm,perm);
   }    }
   MKLIST(*rp,n0);    MKLIST(*rp,n0);
Line 2351  void Pmat_to_gfmmat(NODE arg,GFMMAT *rp)
Line 2603  void Pmat_to_gfmmat(NODE arg,GFMMAT *rp)
   
   asir_assert(ARG0(arg),O_MAT,"mat_to_gfmmat");    asir_assert(ARG0(arg),O_MAT,"mat_to_gfmmat");
   asir_assert(ARG1(arg),O_N,"mat_to_gfmmat");    asir_assert(ARG1(arg),O_N,"mat_to_gfmmat");
   m = (MAT)ARG0(arg); md = (unsigned int)QTOS((Q)ARG1(arg));    m = (MAT)ARG0(arg); md = (unsigned int)ZTOS((Q)ARG1(arg));
   mat_to_gfmmat(m,md,rp);    mat_to_gfmmat(m,md,rp);
 }  }
   
Line 2390  void Pgeninvm_swap(NODE arg,LIST *rp)
Line 2642  void Pgeninvm_swap(NODE arg,LIST *rp)
   
   asir_assert(ARG0(arg),O_MAT,"geninvm_swap");    asir_assert(ARG0(arg),O_MAT,"geninvm_swap");
   asir_assert(ARG1(arg),O_N,"geninvm_swap");    asir_assert(ARG1(arg),O_N,"geninvm_swap");
   m = (MAT)ARG0(arg); md = QTOS((Z)ARG1(arg));    m = (MAT)ARG0(arg); md = ZTOS((Z)ARG1(arg));
   row = m->row; col = m->col; mat = m->body;    row = m->row; col = m->col; mat = m->body;
   wmat = (unsigned int **)almat(row,col+row);    wmat = (unsigned int **)almat(row,col+row);
   for ( i = 0; i < row; i++ ) {    for ( i = 0; i < row; i++ ) {
Line 2406  void Pgeninvm_swap(NODE arg,LIST *rp)
Line 2658  void Pgeninvm_swap(NODE arg,LIST *rp)
     MKMAT(mat1,col,col);      MKMAT(mat1,col,col);
     for ( i = 0, tmat = (Z **)mat1->body; i < col; i++ )      for ( i = 0, tmat = (Z **)mat1->body; i < col; i++ )
       for ( j = 0; j < col; j++ )        for ( j = 0; j < col; j++ )
         UTOQ(invmat[i][j],tmat[i][j]);          UTOZ(invmat[i][j],tmat[i][j]);
     MKVECT(vect1,row);      MKVECT(vect1,row);
     for ( i = 0, tvect = (Z *)vect1->body; i < row; i++ )      for ( i = 0, tvect = (Z *)vect1->body; i < row; i++ )
       STOQ(index[i],tvect[i]);        STOZ(index[i],tvect[i]);
      MKNODE(node2,vect1,0); MKNODE(node1,mat1,node2); MKLIST(*rp,node1);       MKNODE(node2,vect1,0); MKNODE(node1,mat1,node2); MKLIST(*rp,node1);
   }    }
 }  }
Line 2504  void Pgeninv_sf_swap(NODE arg,LIST *rp)
Line 2756  void Pgeninv_sf_swap(NODE arg,LIST *rp)
         }          }
     MKVECT(vect1,row);      MKVECT(vect1,row);
     for ( i = 0, tvect = (Z *)vect1->body; i < row; i++ )      for ( i = 0, tvect = (Z *)vect1->body; i < row; i++ )
       STOQ(index[i],tvect[i]);        STOZ(index[i],tvect[i]);
     MKNODE(node2,vect1,0); MKNODE(node1,mat1,node2); MKLIST(*rp,node1);      MKNODE(node2,vect1,0); MKNODE(node1,mat1,node2); MKLIST(*rp,node1);
   }    }
 }  }
Line 2596  void Pnbpoly_up2(NODE arg,GF2N *rp)
Line 2848  void Pnbpoly_up2(NODE arg,GF2N *rp)
   int m,type,ret;    int m,type,ret;
   UP2 r;    UP2 r;
   
   m = QTOS((Z)ARG0(arg));    m = ZTOS((Z)ARG0(arg));
   type = QTOS((Z)ARG1(arg));    type = ZTOS((Z)ARG1(arg));
   ret = generate_ONB_polynomial(&r,m,type);    ret = generate_ONB_polynomial(&r,m,type);
   if ( ret == 0 )    if ( ret == 0 )
     MKGF2N(r,*rp);      MKGF2N(r,*rp);
Line 2611  void Px962_irredpoly_up2(NODE arg,GF2N *rp)
Line 2863  void Px962_irredpoly_up2(NODE arg,GF2N *rp)
   GF2N prev;    GF2N prev;
   UP2 r;    UP2 r;
   
   m = QTOS((Q)ARG0(arg));    m = ZTOS((Q)ARG0(arg));
   prev = (GF2N)ARG1(arg);    prev = (GF2N)ARG1(arg);
   if ( !prev ) {    if ( !prev ) {
     w = (m>>5)+1; NEWUP2(r,w); r->w = 0;      w = (m>>5)+1; NEWUP2(r,w); r->w = 0;
Line 2636  void Pirredpoly_up2(NODE arg,GF2N *rp)
Line 2888  void Pirredpoly_up2(NODE arg,GF2N *rp)
   GF2N prev;    GF2N prev;
   UP2 r;    UP2 r;
   
   m = QTOS((Q)ARG0(arg));    m = ZTOS((Q)ARG0(arg));
   prev = (GF2N)ARG1(arg);    prev = (GF2N)ARG1(arg);
   if ( !prev ) {    if ( !prev ) {
     w = (m>>5)+1; NEWUP2(r,w); r->w = 0;      w = (m>>5)+1; NEWUP2(r,w); r->w = 0;
Line 2665  void Pmat_swap_row_destructive(NODE arg, MAT *m)
Line 2917  void Pmat_swap_row_destructive(NODE arg, MAT *m)
   asir_assert(ARG1(arg),O_N,"mat_swap_row_destructive");    asir_assert(ARG1(arg),O_N,"mat_swap_row_destructive");
   asir_assert(ARG2(arg),O_N,"mat_swap_row_destructive");    asir_assert(ARG2(arg),O_N,"mat_swap_row_destructive");
   mat = (MAT)ARG0(arg);    mat = (MAT)ARG0(arg);
   i1 = QTOS((Q)ARG1(arg));    i1 = ZTOS((Q)ARG1(arg));
   i2 = QTOS((Q)ARG2(arg));    i2 = ZTOS((Q)ARG2(arg));
   if ( i1 < 0 || i2 < 0 || i1 >= mat->row || i2 >= mat->row )    if ( i1 < 0 || i2 < 0 || i1 >= mat->row || i2 >= mat->row )
     error("mat_swap_row_destructive : Out of range");      error("mat_swap_row_destructive : Out of range");
   t = mat->body[i1];    t = mat->body[i1];
Line 2686  void Pmat_swap_col_destructive(NODE arg, MAT *m)
Line 2938  void Pmat_swap_col_destructive(NODE arg, MAT *m)
   asir_assert(ARG1(arg),O_N,"mat_swap_col_destructive");    asir_assert(ARG1(arg),O_N,"mat_swap_col_destructive");
   asir_assert(ARG2(arg),O_N,"mat_swap_col_destructive");    asir_assert(ARG2(arg),O_N,"mat_swap_col_destructive");
   mat = (MAT)ARG0(arg);    mat = (MAT)ARG0(arg);
   j1 = QTOS((Q)ARG1(arg));    j1 = ZTOS((Q)ARG1(arg));
   j2 = QTOS((Q)ARG2(arg));    j2 = ZTOS((Q)ARG2(arg));
   if ( j1 < 0 || j2 < 0 || j1 >= mat->col || j2 >= mat->col )    if ( j1 < 0 || j2 < 0 || j1 >= mat->col || j2 >= mat->col )
     error("mat_swap_col_destructive : Out of range");      error("mat_swap_col_destructive : Out of range");
   n = mat->row;    n = mat->row;
Line 2999  void Pnd_det(NODE arg,P *rp)
Line 3251  void Pnd_det(NODE arg,P *rp)
   if ( argc(arg) == 1 )    if ( argc(arg) == 1 )
     nd_det(0,ARG0(arg),rp);      nd_det(0,ARG0(arg),rp);
   else    else
     nd_det(QTOS((Q)ARG1(arg)),ARG0(arg),rp);      nd_det(ZTOS((Q)ARG1(arg)),ARG0(arg),rp);
 }  }
   
 void Pmat_col(NODE arg,VECT *rp)  void Pmat_col(NODE arg,VECT *rp)
Line 3011  void Pmat_col(NODE arg,VECT *rp)
Line 3263  void Pmat_col(NODE arg,VECT *rp)
   asir_assert(ARG0(arg),O_MAT,"mat_col");    asir_assert(ARG0(arg),O_MAT,"mat_col");
   asir_assert(ARG1(arg),O_N,"mat_col");    asir_assert(ARG1(arg),O_N,"mat_col");
   mat = (MAT)ARG0(arg);    mat = (MAT)ARG0(arg);
   j = QTOS((Q)ARG1(arg));    j = ZTOS((Q)ARG1(arg));
   if ( j < 0 || j >= mat->col) {    if ( j < 0 || j >= mat->col) {
     error("mat_col : Out of range");      error("mat_col : Out of range");
   }    }

Legend:
Removed from v.1.1  
changed lines
  Added in v.1.7

FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>