=================================================================== RCS file: /home/cvs/OpenXM_contrib2/asir2018/builtin/array.c,v retrieving revision 1.2 retrieving revision 1.9 diff -u -p -r1.2 -r1.9 --- OpenXM_contrib2/asir2018/builtin/array.c 2018/09/28 08:20:27 1.2 +++ OpenXM_contrib2/asir2018/builtin/array.c 2021/03/26 09:05:41 1.9 @@ -45,7 +45,7 @@ * DEVELOPER SHALL HAVE NO LIABILITY IN CONNECTION WITH THE USE, * PERFORMANCE OR NON-PERFORMANCE OF THE SOFTWARE. * - * $OpenXM: OpenXM_contrib2/asir2018/builtin/array.c,v 1.1 2018/09/19 05:45:05 noro Exp $ + * $OpenXM: OpenXM_contrib2/asir2018/builtin/array.c,v 1.8 2020/10/06 06:31:19 noro 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) { @@ -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 && 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) { NODE n0; @@ -1239,11 +1305,20 @@ 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; 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 = ZTOS((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; wmat = (int **)almat(row,col); @@ -1284,6 +1359,7 @@ void Pgeneric_gauss_elim_mod(NODE arg,LIST *rp) MKLIST(*rp,n0); } + void Pleqm(NODE arg,VECT *rp) { MAT m; @@ -1354,7 +1430,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 +1543,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 +1559,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 +1585,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 +1620,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 +1634,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 +1658,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 +1711,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 +1727,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 +1944,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 +1990,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 +2001,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 +2050,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 +2061,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 +2106,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 +2140,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 +2150,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 +2186,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 +2242,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 +2277,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 +2306,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 */ @@ -2215,7 +2467,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 +2476,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); } @@ -2442,7 +2694,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 +2704,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 +2740,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,7 +2751,7 @@ 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); @@ -2536,7 +2788,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 +2800,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); @@ -2721,7 +2973,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<