version 1.1, 2018/09/19 05:45:05 |
version 1.9, 2021/03/26 09:05:41 |
|
|
* 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.8 2020/10/06 06:31:19 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 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 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 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 1467 void lu_dec_cr(MAT mat,MAT lu,Q *dn,int **perm) |
|
Line 1543 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 1483 void reduce_reducers_mod(int **mat,int row,int col,int |
|
Line 1559 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 1509 void reduce_reducers_mod(int **mat,int row,int col,int |
|
Line 1585 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 1544 void pre_reduce_mod(int **mat,int row,int col,int nred |
|
Line 1620 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 1558 void pre_reduce_mod(int **mat,int row,int col,int nred |
|
Line 1634 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 1582 void reduce_sp_by_red_mod(int *sp,int **redmat,int *in |
|
Line 1658 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 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] ) != 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 |
#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 1692 void reduce_sp_by_red_mod_compress (int *sp,CDP *redma |
|
Line 1944 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 1738 int generic_gauss_elim_mod(int **mat0,int row,int col, |
|
Line 1990 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 1749 int generic_gauss_elim_mod(int **mat0,int row,int col, |
|
Line 2001 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 1798 int generic_gauss_elim_mod2(int **mat0,int row,int col |
|
Line 2050 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 1809 int generic_gauss_elim_mod2(int **mat0,int row,int col |
|
Line 2061 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 1854 int indep_rows_mod(int **mat0,int row,int col,int md,i |
|
Line 2106 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 1888 int generic_gauss_elim_sf(int **mat0,int row,int col,i |
|
Line 2140 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 1898 int generic_gauss_elim_sf(int **mat0,int row,int col,i |
|
Line 2150 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 1934 int lu_gfmmat(GFMMAT mat,unsigned int md,int *perm) |
|
Line 2186 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 1990 int find_lhs_and_lu_mod(unsigned int **a,int row,int c |
|
Line 2242 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 2025 int lu_mod(unsigned int **a,int n,unsigned int md,int |
|
Line 2277 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 2054 void solve_by_lu_mod(int **a,int n,int md,int **b,int |
|
Line 2306 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 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 2215 int gauss_elim_geninv_mod(unsigned int **mat,int row,i |
|
Line 2467 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 2224 int gauss_elim_geninv_mod(unsigned int **mat,int row,i |
|
Line 2476 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 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 2442 int gauss_elim_geninv_mod_swap(unsigned int **mat,int |
|
Line 2694 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 2452 int gauss_elim_geninv_mod_swap(unsigned int **mat,int |
|
Line 2704 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 2488 void Pgeninv_sf_swap(NODE arg,LIST *rp) |
|
Line 2740 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 2499 void Pgeninv_sf_swap(NODE arg,LIST *rp) |
|
Line 2751 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); |
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 2536 int gauss_elim_geninv_sf_swap(int **mat,int row,int co |
|
Line 2788 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 2548 int gauss_elim_geninv_sf_swap(int **mat,int row,int co |
|
Line 2800 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 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 2721 int generate_ONB_polynomial(UP2 *rp,int m,int type) |
|
Line 2973 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; |
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"); |
} |
} |