version 1.2, 2018/09/28 08:20:27 |
version 1.5, 2019/03/28 06:44:04 |
|
|
* 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.1 2018/09/19 05:45:05 noro Exp $ |
* $OpenXM: OpenXM_contrib2/asir2018/builtin/array.c,v 1.4 2018/10/19 23:27:38 noro Exp $ |
*/ |
*/ |
#include "ca.h" |
#include "ca.h" |
#include "base.h" |
#include "base.h" |
Line 1229 void Pindep_rows_mod(NODE arg,VECT *rp) |
|
Line 1229 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 1298 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 = 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; |
row = m->row; col = m->col; tmat = (Z **)m->body; |
wmat = (int **)almat(row,col); |
wmat = (int **)almat(row,col); |
|
|
Line 1284 void Pgeneric_gauss_elim_mod(NODE arg,LIST *rp) |
|
Line 1352 void Pgeneric_gauss_elim_mod(NODE arg,LIST *rp) |
|
MKLIST(*rp,n0); |
MKLIST(*rp,n0); |
} |
} |
|
|
|
|
void Pleqm(NODE arg,VECT *rp) |
void Pleqm(NODE arg,VECT *rp) |
{ |
{ |
MAT m; |
MAT m; |
Line 1635 void red_by_vect(int m,unsigned int *p,unsigned int *r |
|
Line 1704 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 1650 void red_by_vect64(int m, U64 *p,unsigned int *c,U64 * |
|
Line 1719 void red_by_vect64(int m, U64 *p,unsigned int *c,U64 * |
|
if ( t < *p ) (*c)++; |
if ( t < *p ) (*c)++; |
*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 |
|
|