=================================================================== RCS file: /home/cvs/OpenXM_contrib2/asir2018/builtin/array.c,v retrieving revision 1.1 retrieving revision 1.10 diff -u -p -r1.1 -r1.10 --- OpenXM_contrib2/asir2018/builtin/array.c 2018/09/19 05:45:05 1.1 +++ OpenXM_contrib2/asir2018/builtin/array.c 2022/01/13 08:15:02 1.10 @@ -45,7 +45,7 @@ * DEVELOPER SHALL HAVE NO LIABILITY IN CONNECTION WITH THE USE, * PERFORMANCE OR NON-PERFORMANCE OF THE SOFTWARE. * - * $OpenXM$ + * $OpenXM: OpenXM_contrib2/asir2018/builtin/array.c,v 1.9 2021/03/26 09:05:41 ohara Exp $ */ #include "ca.h" #include "base.h" @@ -170,6 +170,13 @@ void solve_u(int *,ent **,int,int *,int); static int *ul,*ll; static ent **u,**l; 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) { @@ -206,14 +213,14 @@ void Plusolve_main(NODE arg,VECT *rp) v = (VECT)ARG0(arg); len = v->len; d = (Z *)BDY(v); 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_u(ul,u,len,rhs,modulus); NEWVECT(r); r->len = len; r->body = (pointer *)MALLOC(len*sizeof(pointer)); p = (Z *)r->body; for ( i = 0; i < len; i++ ) - STOQ(rhs[i],p[i]); + STOZ(rhs[i],p[i]); *rp = r; } @@ -565,7 +572,7 @@ void Psepmat_destructive(NODE arg,LIST *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) @@ -601,7 +608,7 @@ void Pnewvect(NODE arg,VECT *rp) NODE tn; asir_assert(ARG0(arg),O_N,"newvect"); - len = QTOS((Q)ARG0(arg)); + len = ZTOS((Q)ARG0(arg)); if ( len < 0 ) error("newvect : invalid size"); MKVECT(vect,len); @@ -615,7 +622,7 @@ void Pnewvect(NODE arg,VECT *rp) return; } #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 && ibody; r0 = 0; for ( j = 0; j < y; j++ ) for ( i = 0; i < len; i++ ) if ( MEMORY_GETPOINT(a,blen,i,j) ) { NEXTNODE(r0,r); - STOQ(i,iq); STOQ(j,jq); + STOZ(i,iq); STOZ(j,jq); n = mknode(2,iq,jq); MKLIST(l,n); BDY(r) = l; @@ -775,7 +782,7 @@ void Pnewmat(NODE arg,MAT *rp) asir_assert(ARG0(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 ) error("newmat : invalid size"); MKMAT(m,row,col); @@ -967,7 +974,7 @@ void Premainder(NODE arg,Obj *rp) cmp(md,(P)a,(P *)rp); break; case O_VECT: /* strage spec */ - smd = QTOS(md); + smd = ZTOS(md); v = (VECT)a; n = v->len; vb = v->body; MKVECT(w,n); wb = w->body; for ( i = 0; i < n; i++ ) { @@ -1015,7 +1022,7 @@ void Psremainder(NODE arg,Obj *rp) case O_P: cmp(md,(P)a,(P *)rp); break; case O_VECT: - smd = QTOS(md); + smd = ZTOS(md); v = (VECT)a; n = v->len; vb = v->body; MKVECT(w,n); wb = w->body; for ( i = 0; i < n; i++ ) @@ -1049,15 +1056,15 @@ void Psize(NODE arg,LIST *rp) switch (OID(ARG0(arg))) { case O_VECT: n = ((VECT)ARG0(arg))->len; - STOQ(n,q); MKNODE(t,q,0); + STOZ(n,q); MKNODE(t,q,0); break; case O_MAT: 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; case O_IMAT: 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; default: error("size : invalid argument"); break; @@ -1080,7 +1087,7 @@ void Pdet(NODE arg,P *rp) else if ( argc(arg) == 1 ) detp(CO,(P **)BDY(m),m->row,rp); 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); for ( i = 0; i < n; i++ ) for ( j = 0; j < n; j++ ) @@ -1109,7 +1116,7 @@ void Pinvmat(NODE arg,LIST *rp) nd = mknode(2,r,dn); MKLIST(*rp,nd); } 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); for ( i = 0; i < n; i++ ) for ( j = 0; j < n; j++ ) @@ -1169,11 +1176,11 @@ void Pgeneric_gauss_elim(NODE arg,LIST *rp) MKVECT(rind,rank); MKVECT(cind,t); for ( i = 0; i < rank; i++ ) { - STOQ(ri[i],q); + STOZ(ri[i],q); BDY(rind)[i] = (pointer)q; } for ( i = 0; i < t; i++ ) { - STOQ(ci[i],q); + STOZ(ci[i],q); BDY(cind)[i] = (pointer)q; } n0 = mknode(4,nm,dn,rind,cind); @@ -1195,7 +1202,7 @@ void Pindep_rows_mod(NODE arg,VECT *rp) asir_assert(ARG0(arg),O_MAT,"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; wmat = (int **)almat(row,col); @@ -1211,7 +1218,7 @@ void Pindep_rows_mod(NODE arg,VECT *rp) MKVECT(rind,rank); rib = (Z *)rind->body; for ( j = 0; j < rank; j++ ) { - STOQ(rowstat[j],rib[j]); + STOZ(rowstat[j],rib[j]); } *rp = rind; } @@ -1229,6 +1236,78 @@ 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]}+... */ +#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; + Obj val; + int asis = 0; + + asir_assert(ARG0(arg),O_MAT,"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)); + 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); + if ( asis ) { + MKMAT(mat,row,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); + 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) { NODE n0; @@ -1239,11 +1318,23 @@ void Pgeneric_gauss_elim_mod(NODE arg,LIST *rp) Z *rib,*cib,*rnb; int *colstat,*p; Z q; + long mdl; 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(ARG1(arg),O_N,"generic_gauss_elim_mod"); - m = (MAT)ARG0(arg); md = QTOS((Z)ARG1(arg)); + if ( get_opt("asis",&val) && val ) asis = 1; +#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; wmat = (int **)almat(row,col); @@ -1255,20 +1346,30 @@ void Pgeneric_gauss_elim_mod(NODE arg,LIST *rp) for ( j = 0; j < col; j++ ) wmat[i][j] = remqi((Q)tmat[i][j],md); rank = generic_gauss_elim_mod(wmat,row,col,md,colstat); + if ( asis ) { + MKMAT(mat,row,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); rnb = (Z *)rnum->body; for ( i = 0; i < rank; i++ ) for ( j = 0, p = wmat[i]; j < row; j++ ) if ( p == row0[j] ) - STOQ(j,rnb[i]); + 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] ) { - UTOQ(wmat[i][j],tmat[i][k]); k++; + UTOZ(wmat[i][j],tmat[i][k]); k++; } MKVECT(rind,rank); @@ -1276,14 +1377,15 @@ void Pgeneric_gauss_elim_mod(NODE arg,LIST *rp) rib = (Z *)rind->body; cib = (Z *)cind->body; for ( j = k = l = 0; j < col; j++ ) if ( colstat[j] ) { - STOQ(j,rib[k]); k++; + STOZ(j,rib[k]); k++; } else { - STOQ(j,cib[l]); l++; + STOZ(j,cib[l]); l++; } n0 = mknode(4,mat,rind,cind,rnum); MKLIST(*rp,n0); } + void Pleqm(NODE arg,VECT *rp) { MAT m; @@ -1296,7 +1398,7 @@ void Pleqm(NODE arg,VECT *rp) asir_assert(ARG0(arg),O_MAT,"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; wmat = (int **)almat(row,col); for ( i = 0; i < row; i++ ) @@ -1311,7 +1413,7 @@ void Pleqm(NODE arg,VECT *rp) n = col - 1; MKVECT(vect,n); 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; } @@ -1354,7 +1456,7 @@ int gauss_elim_mod(int **mat,int row,int col,int md) 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; #if 0 @@ -1467,7 +1569,7 @@ void lu_dec_cr(MAT mat,MAT lu,Q *dn,int **perm) 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) { @@ -1483,7 +1585,7 @@ void reduce_reducers_mod(int **mat,int row,int col,int ind[i] = j; for ( l = i-1; l >= 0; 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] */ j = ind[l]; s = mat[l]+j; @@ -1509,7 +1611,7 @@ void reduce_reducers_mod(int **mat,int row,int col,int ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1 } 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++; } } } @@ -1544,7 +1646,7 @@ void pre_reduce_mod(int **mat,int row,int col,int nred DMAR(t[k],inv,0,md,t[k]) for ( l = i-1; l >= 0; 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] */ for ( k = ind[l], hc = md-hc, s = mat[l]+k, tk = t+k; k < col; k++, tk++, s++ ) @@ -1558,7 +1660,7 @@ void pre_reduce_mod(int **mat,int row,int col,int nred t = mat[i]; for ( l = nred-1; l >= 0; 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] */ for ( k = ind[l], hc = md-hc, s = mat[l]+k, tk = t+k; k < col; k++, tk++, s++ ) @@ -1582,14 +1684,14 @@ void reduce_sp_by_red_mod(int *sp,int **redmat,int *in /* reduce the spolys by redmat */ for ( i = nred-1; i >= 0; i-- ) { /* reduce sp by redmat[i] */ - if ( hc = sp[ind[i]] ) { + if ( ( hc = sp[ind[i]] ) != 0 ) { /* sp = sp-hc*redmat[i] */ j = ind[i]; hc = md-hc; s = redmat[i]+j; tj = sp+j; 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++; } } } @@ -1635,12 +1737,12 @@ void red_by_vect(int m,unsigned int *p,unsigned int *r } } -#if defined(__GNUC__) && SIZEOF_LONG==8 -/* 64bit vector += UNIT vector(normalized) */ +#if SIZEOF_LONG==8 +/* 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; r++; len--; @@ -1651,6 +1753,182 @@ void red_by_vect64(int m, U64 *p,unsigned int *c,U64 * *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] ) != 0 ) + 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] ) != 0 ) + 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 void red_by_vect_sf(int m,unsigned int *p,unsigned int *r,unsigned int hc,int len) @@ -1692,7 +1970,7 @@ void reduce_sp_by_red_mod_compress (int *sp,CDP *redma for ( i = nred-1; i >= 0; i-- ) { /* reduce sp by redmat[i] */ usp[ind[i]] %= md; - if ( hc = usp[ind[i]] ) { + if ( ( hc = usp[ind[i]] ) != 0 ) { /* sp = sp-hc*redmat[i] */ hc = md-hc; ri = redmat[i]; @@ -1738,7 +2016,7 @@ int generic_gauss_elim_mod(int **mat0,int row,int col, } for ( i = rank+1; i < row; 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); } rank++; @@ -1749,7 +2027,7 @@ int generic_gauss_elim_mod(int **mat0,int row,int col, for ( i = 0; i < l; i++ ) { t = mat[i]; t[j] %= md; - if ( a = t[j] ) + if ( ( a = t[j] ) != 0 ) red_by_vect(md,t+j,pivot+j,md-a,col-j); } l--; @@ -1798,7 +2076,7 @@ int generic_gauss_elim_mod2(int **mat0,int row,int col } for ( i = rank+1; i < row; 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); } rank++; @@ -1809,7 +2087,7 @@ int generic_gauss_elim_mod2(int **mat0,int row,int col for ( i = 0; i < l; i++ ) { t = mat[i]; t[j] %= md; - if ( a = t[j] ) + if ( ( a = t[j] ) != 0 ) red_by_vect(md,t+j,pivot+j,md-a,col-j); } l--; @@ -1854,7 +2132,7 @@ int indep_rows_mod(int **mat0,int row,int col,int md,i } for ( i = rank+1; i < row; 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); } rank++; @@ -1888,7 +2166,7 @@ int generic_gauss_elim_sf(int **mat0,int row,int col,i *pk = _mulsf(*pk,inv); for ( i = rank+1; i < row; 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); } rank++; @@ -1898,7 +2176,7 @@ int generic_gauss_elim_sf(int **mat0,int row,int col,i pivot = mat[l]; for ( i = 0; i < l; 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); } l--; @@ -1934,7 +2212,7 @@ int lu_gfmmat(GFMMAT mat,unsigned int md,int *perm) pivot[k] = inv = invm(pivot[k],md); for ( i = k+1; i < row; i++ ) { t = a[i]; - if ( m = t[k] ) { + if ( ( m = t[k] ) != 0 ) { DMAR(inv,m,0,md,t[k]) for ( j = k+1, m = md - t[k]; j < col; j++ ) if ( pivot[j] ) { @@ -1990,7 +2268,7 @@ int find_lhs_and_lu_mod(unsigned int **a,int row,int c pivot[k] = inv = invm(pivot[k],md); for ( i = d+1; i < row; i++ ) { t = a[i]; - if ( m = t[k] ) { + if ( ( m = t[k] ) != 0 ) { DMAR(inv,m,0,md,t[k]) for ( j = k+1, m = md - t[k]; j < col; j++ ) if ( pivot[j] ) { @@ -2025,7 +2303,7 @@ int lu_mod(unsigned int **a,int n,unsigned int md,int inv = invm(pivot[k],md); for ( i = k+1; i < n; i++ ) { t = a[i]; - if ( m = t[k] ) { + if ( ( m = t[k] ) != 0 ) { DMAR(inv,m,0,md,t[k]) for ( j = k+1, m = md - t[k]; j < n; j++ ) if ( pivot[j] ) { @@ -2054,8 +2332,8 @@ void solve_by_lu_mod(int **a,int n,int md,int **b,int int i,j,k; unsigned int t,m,m2; - y = (int *)MALLOC_ATOMIC(n*sizeof(int)); - c = (int *)MALLOC_ATOMIC(n*sizeof(int)); + y = (unsigned int *)MALLOC_ATOMIC(n*sizeof(int)); + c = (unsigned int *)MALLOC_ATOMIC(n*sizeof(int)); m2 = md>>1; for ( k = 0; k < l; k++ ) { /* copy b[.][k] to c */ @@ -2102,7 +2380,7 @@ void Pleqm1(NODE arg,VECT *rp) asir_assert(ARG0(arg),O_MAT,"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; wmat = (int **)almat(row,col); for ( i = 0; i < row; i++ ) @@ -2117,7 +2395,7 @@ void Pleqm1(NODE arg,VECT *rp) n = col - 1; MKVECT(vect,n); 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; } @@ -2173,7 +2451,7 @@ void Pgeninvm(NODE arg,LIST *rp) asir_assert(ARG0(arg),O_MAT,"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; wmat = (unsigned int **)almat(row,col+row); for ( i = 0; i < row; i++ ) { @@ -2188,10 +2466,10 @@ void Pgeninvm(NODE arg,LIST *rp) MKMAT(mat1,col,row); MKMAT(mat2,row-col,row); for ( i = 0, tmat = (Z **)mat1->body; i < col; i++ ) 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 ( 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); } } @@ -2215,7 +2493,7 @@ int gauss_elim_geninv_mod(unsigned int **mat,int row,i pivot[k] = dmar(pivot[k],inv,0,md); for ( i = j+1; i < row; i++ ) { t = mat[i]; - if ( a = t[j] ) + if ( ( a = t[j] ) != 0 ) for ( k = j, a = md - a; k < m; k++ ) t[k] = dmar(pivot[k],a,t[k],md); } @@ -2224,7 +2502,7 @@ int gauss_elim_geninv_mod(unsigned int **mat,int row,i pivot = mat[j]; for ( i = j-1; i >= 0; i-- ) { t = mat[i]; - if ( a = t[j] ) + if ( ( a = t[j] ) != 0 ) for ( k = j, a = md - a; k < m; k++ ) t[k] = dmar(pivot[k],a,t[k],md); } @@ -2244,16 +2522,16 @@ void Psolve_by_lu_gfmmat(NODE arg,VECT *rp) lu = (GFMMAT)ARG0(arg); perm = (Z *)BDY((VECT)ARG1(arg)); rhs = (Z *)BDY((VECT)ARG2(arg)); - md = (unsigned int)QTOS((Z)ARG3(arg)); + md = (unsigned int)ZTOS((Z)ARG3(arg)); n = lu->col; b = (unsigned int *)MALLOC_ATOMIC(n*sizeof(int)); sol = (unsigned int *)MALLOC_ATOMIC(n*sizeof(int)); 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); MKVECT(r,n); for ( i = 0, v = (Z *)r->body; i < n; i++ ) - UTOQ(sol[i],v[i]); + UTOZ(sol[i],v[i]); *rp = r; } @@ -2308,7 +2586,7 @@ void Plu_mat(NODE arg,LIST *rp) lu_dec_cr(m,lu,&dn,&iperm); MKVECT(perm,n); 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); MKLIST(*rp,n0); } @@ -2327,7 +2605,7 @@ void Plu_gfmmat(NODE arg,LIST *rp) asir_assert(ARG0(arg),O_MAT,"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); row = m->row; col = m->col; @@ -2338,7 +2616,7 @@ void Plu_gfmmat(NODE arg,LIST *rp) else { MKVECT(perm,row); 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); } MKLIST(*rp,n0); @@ -2351,7 +2629,7 @@ void Pmat_to_gfmmat(NODE arg,GFMMAT *rp) asir_assert(ARG0(arg),O_MAT,"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); } @@ -2390,7 +2668,7 @@ void Pgeninvm_swap(NODE arg,LIST *rp) asir_assert(ARG0(arg),O_MAT,"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; wmat = (unsigned int **)almat(row,col+row); for ( i = 0; i < row; i++ ) { @@ -2406,10 +2684,10 @@ void Pgeninvm_swap(NODE arg,LIST *rp) MKMAT(mat1,col,col); for ( i = 0, tmat = (Z **)mat1->body; i < col; i++ ) for ( j = 0; j < col; j++ ) - UTOQ(invmat[i][j],tmat[i][j]); + UTOZ(invmat[i][j],tmat[i][j]); MKVECT(vect1,row); 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); } } @@ -2442,7 +2720,7 @@ int gauss_elim_geninv_mod_swap(unsigned int **mat,int pivot[k] = (unsigned int)dmar(pivot[k],inv,0,md); for ( i = j+1; i < row; i++ ) { t = mat[i]; - if ( a = t[j] ) + if ( ( a = t[j] ) != 0 ) for ( k = j, a = md - a; k < m; k++ ) if ( pivot[k] ) t[k] = dmar(pivot[k],a,t[k],md); @@ -2452,7 +2730,7 @@ int gauss_elim_geninv_mod_swap(unsigned int **mat,int pivot = mat[j]; for ( i = j-1; i >= 0; i-- ) { t = mat[i]; - if ( a = t[j] ) + if ( ( a = t[j] ) != 0 ) for ( k = j, a = md - a; k < m; k++ ) if ( pivot[k] ) t[k] = dmar(pivot[k],a,t[k],md); @@ -2488,7 +2766,7 @@ void Pgeninv_sf_swap(NODE arg,LIST *rp) for ( i = 0; i < row; i++ ) { bzero((char *)wmat[i],(col+row)*sizeof(int)); 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][col+i] = _onesf(); } @@ -2499,12 +2777,12 @@ void Pgeninv_sf_swap(NODE arg,LIST *rp) MKMAT(mat1,col,col); for ( i = 0, tmat = (GFS **)mat1->body; i < col; i++ ) for ( j = 0; j < col; j++ ) - if ( t = invmat[i][j] ) { + if ( ( t = invmat[i][j] ) != 0 ) { MKGFS(IFTOF(t),tmat[i][j]); } MKVECT(vect1,row); 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); } } @@ -2536,7 +2814,7 @@ int gauss_elim_geninv_sf_swap(int **mat,int row,int co pivot[k] = _mulsf(pivot[k],inv); for ( i = j+1; i < row; i++ ) { t = mat[i]; - if ( a = t[j] ) + if ( ( a = t[j] ) != 0 ) for ( k = j, a = _chsgnsf(a); k < m; k++ ) if ( pivot[k] ) { u = _mulsf(pivot[k],a); @@ -2548,7 +2826,7 @@ int gauss_elim_geninv_sf_swap(int **mat,int row,int co pivot = mat[j]; for ( i = j-1; i >= 0; i-- ) { t = mat[i]; - if ( a = t[j] ) + if ( ( a = t[j] ) != 0 ) for ( k = j, a = _chsgnsf(a); k < m; k++ ) if ( pivot[k] ) { u = _mulsf(pivot[k],a); @@ -2596,8 +2874,8 @@ void Pnbpoly_up2(NODE arg,GF2N *rp) int m,type,ret; UP2 r; - m = QTOS((Z)ARG0(arg)); - type = QTOS((Z)ARG1(arg)); + m = ZTOS((Z)ARG0(arg)); + type = ZTOS((Z)ARG1(arg)); ret = generate_ONB_polynomial(&r,m,type); if ( ret == 0 ) MKGF2N(r,*rp); @@ -2611,7 +2889,7 @@ void Px962_irredpoly_up2(NODE arg,GF2N *rp) GF2N prev; UP2 r; - m = QTOS((Q)ARG0(arg)); + m = ZTOS((Q)ARG0(arg)); prev = (GF2N)ARG1(arg); if ( !prev ) { w = (m>>5)+1; NEWUP2(r,w); r->w = 0; @@ -2636,7 +2914,7 @@ void Pirredpoly_up2(NODE arg,GF2N *rp) GF2N prev; UP2 r; - m = QTOS((Q)ARG0(arg)); + m = ZTOS((Q)ARG0(arg)); prev = (GF2N)ARG1(arg); if ( !prev ) { w = (m>>5)+1; NEWUP2(r,w); r->w = 0; @@ -2665,8 +2943,8 @@ void Pmat_swap_row_destructive(NODE arg, MAT *m) asir_assert(ARG1(arg),O_N,"mat_swap_row_destructive"); asir_assert(ARG2(arg),O_N,"mat_swap_row_destructive"); mat = (MAT)ARG0(arg); - i1 = QTOS((Q)ARG1(arg)); - i2 = QTOS((Q)ARG2(arg)); + i1 = ZTOS((Q)ARG1(arg)); + i2 = ZTOS((Q)ARG2(arg)); if ( i1 < 0 || i2 < 0 || i1 >= mat->row || i2 >= mat->row ) error("mat_swap_row_destructive : Out of range"); t = mat->body[i1]; @@ -2686,8 +2964,8 @@ void Pmat_swap_col_destructive(NODE arg, MAT *m) asir_assert(ARG1(arg),O_N,"mat_swap_col_destructive"); asir_assert(ARG2(arg),O_N,"mat_swap_col_destructive"); mat = (MAT)ARG0(arg); - j1 = QTOS((Q)ARG1(arg)); - j2 = QTOS((Q)ARG2(arg)); + j1 = ZTOS((Q)ARG1(arg)); + j2 = ZTOS((Q)ARG2(arg)); if ( j1 < 0 || j2 < 0 || j1 >= mat->col || j2 >= mat->col ) error("mat_swap_col_destructive : Out of range"); n = mat->row; @@ -2721,7 +2999,7 @@ int generate_ONB_polynomial(UP2 *rp,int m,int type) for ( i = 0; i < w; i++ ) f->b[i] = 0xffffffff; /* mask the top word if necessary */ - if ( r = (m+1)&31 ) + if ( ( r = (m+1)&31 ) != 0 ) f->b[w-1] &= (1<= mat->col) { error("mat_col : Out of range"); }