version 1.14, 2001/09/04 05:14:03 |
version 1.23, 2001/10/01 01:58:01 |
|
|
* 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/asir2000/builtin/array.c,v 1.13 2001/06/07 05:14:48 noro Exp $ |
* $OpenXM: OpenXM_contrib2/asir2000/builtin/array.c,v 1.22 2001/09/17 08:37:30 noro Exp $ |
*/ |
*/ |
#include "ca.h" |
#include "ca.h" |
#include "base.h" |
#include "base.h" |
Line 67 void mat_to_gfmmat(MAT,unsigned int,GFMMAT *); |
|
Line 67 void mat_to_gfmmat(MAT,unsigned int,GFMMAT *); |
|
|
|
int generic_gauss_elim_mod(int **,int,int,int,int *); |
int generic_gauss_elim_mod(int **,int,int,int,int *); |
int generic_gauss_elim(MAT ,MAT *,Q *,int **,int **); |
int generic_gauss_elim(MAT ,MAT *,Q *,int **,int **); |
|
void reduce_sp_by_red_mod_compress (int *,CDP *,int *,int,int,int); |
|
|
int gauss_elim_mod(int **,int,int,int); |
int gauss_elim_mod(int **,int,int,int); |
int gauss_elim_mod1(int **,int,int,int); |
int gauss_elim_mod1(int **,int,int,int); |
int gauss_elim_geninv_mod(unsigned int **,int,int,int); |
int gauss_elim_geninv_mod(unsigned int **,int,int,int); |
int gauss_elim_geninv_mod_swap(unsigned int **,int,int,unsigned int,unsigned int ***,int **); |
int gauss_elim_geninv_mod_swap(unsigned int **,int,int,unsigned int,unsigned int ***,int **); |
void Pnewvect(), Pnewmat(), Psepvect(), Psize(), Pdet(), Pleqm(), Pleqm1(), Pgeninvm(); |
void Pnewvect(), Pnewmat(), Psepvect(), Psize(), Pdet(), Pleqm(), Pleqm1(), Pgeninvm(); |
|
void Pinvmat(); |
void Pnewbytearray(); |
void Pnewbytearray(); |
|
|
void Pgeneric_gauss_elim_mod(); |
void Pgeneric_gauss_elim_mod(); |
Line 109 struct ftab array_tab[] = { |
|
Line 111 struct ftab array_tab[] = { |
|
{"vtol",Pvtol,1}, |
{"vtol",Pvtol,1}, |
{"size",Psize,1}, |
{"size",Psize,1}, |
{"det",Pdet,-2}, |
{"det",Pdet,-2}, |
|
{"invmat",Pinvmat,-2}, |
{"leqm",Pleqm,2}, |
{"leqm",Pleqm,2}, |
{"leqm1",Pleqm1,2}, |
{"leqm1",Pleqm1,2}, |
{"geninvm",Pgeninvm,2}, |
{"geninvm",Pgeninvm,2}, |
|
|
} |
} |
} |
} |
|
|
|
void Pinvmat(arg,rp) |
|
NODE arg; |
|
LIST *rp; |
|
{ |
|
MAT m,r; |
|
int n,i,j,mod; |
|
P dn; |
|
P **mat,**imat,**w; |
|
NODE nd; |
|
|
|
m = (MAT)ARG0(arg); |
|
asir_assert(m,O_MAT,"invmat"); |
|
if ( m->row != m->col ) |
|
error("invmat : non-square matrix"); |
|
else if ( argc(arg) == 1 ) { |
|
n = m->row; |
|
invmatp(CO,(P **)BDY(m),n,&imat,&dn); |
|
NEWMAT(r); r->row = n; r->col = n; r->body = (pointer **)imat; |
|
nd = mknode(2,r,dn); |
|
MKLIST(*rp,nd); |
|
} else { |
|
n = m->row; mod = QTOS((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++ ) |
|
ptomp(mod,mat[i][j],&w[i][j]); |
|
#if 0 |
|
detmp(CO,mod,w,n,&d); |
|
mptop(d,rp); |
|
#else |
|
error("not implemented yet"); |
|
#endif |
|
} |
|
} |
|
|
/* |
/* |
input : a row x col matrix A |
input : a row x col matrix A |
A[I] <-> A[I][0]*x_0+A[I][1]*x_1+... |
A[I] <-> A[I][0]*x_0+A[I][1]*x_1+... |
|
|
ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1 |
ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1 |
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 ) { DMAR(zzz,hc,*tj,md,*tj) } tj++; s++; |
} |
} |
} |
} |
|
|
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 ) { DMAR(zzz,hc,*tj,md,*tj) } tj++; s++; |
} |
} |
} |
} |
} |
} |
} |
} |
|
|
|
/* |
|
mat[i] : compressed reducers (i=0,...,nred-1) |
|
mat[0] < mat[1] < ... < mat[nred-1] w.r.t the term order |
|
*/ |
|
|
|
int red_by_compress(m,p,r,ri,hc,len) |
|
int m; |
|
unsigned int *p; |
|
register unsigned int *r; |
|
register unsigned int *ri; |
|
unsigned int hc; |
|
register int len; |
|
{ |
|
unsigned int up,lo; |
|
unsigned int dmy; |
|
unsigned int *pj; |
|
|
|
p[*ri] = 0; r++; ri++; |
|
for ( len--; len; len--, r++, ri++ ) { |
|
pj = p+ *ri; |
|
DMA(*r,hc,*pj,up,lo); |
|
if ( up ) { |
|
DSAB(m,up,lo,dmy,*pj); |
|
} else |
|
*pj = lo; |
|
} |
|
} |
|
|
|
/* p -= hc*r */ |
|
|
|
int red_by_vect(m,p,r,hc,len) |
|
int m; |
|
unsigned int *p,*r; |
|
unsigned int hc; |
|
int len; |
|
{ |
|
register unsigned int up,lo; |
|
unsigned int dmy; |
|
|
|
*p++ = 0; r++; len--; |
|
for ( ; len; len--, r++, p++ ) |
|
if ( *r ) { |
|
DMA(*r,hc,*p,up,lo); |
|
if ( up ) { |
|
DSAB(m,up,lo,dmy,*p); |
|
} else |
|
*p = lo; |
|
} |
|
} |
|
|
|
extern unsigned int **psca; |
|
|
|
void reduce_sp_by_red_mod_compress (sp,redmat,ind,nred,col,md) |
|
int *sp; |
|
CDP *redmat; |
|
int *ind; |
|
int nred,col; |
|
int md; |
|
{ |
|
int i,j,k,len; |
|
unsigned int *tj; |
|
CDP ri; |
|
unsigned int hc,up,lo,up1,lo1,c; |
|
unsigned int *usp; |
|
|
|
usp = (unsigned int *)sp; |
|
/* reduce the spolys by redmat */ |
|
for ( i = nred-1; i >= 0; i-- ) { |
|
/* reduce sp by redmat[i] */ |
|
usp[ind[i]] %= md; |
|
if ( hc = usp[ind[i]] ) { |
|
/* sp = sp-hc*redmat[i] */ |
|
hc = md-hc; |
|
ri = redmat[i]; |
|
len = ri->len; |
|
red_by_compress(md,usp,psca[ri->psindex],ri->body,hc,len); |
|
} |
|
} |
|
for ( i = 0; i < col; i++ ) |
|
if ( usp[i] >= md ) |
|
usp[i] %= md; |
|
} |
|
|
#define ONE_STEP2 if ( zzz = *pk ) { DMAR(zzz,a,*tk,md,*tk) } pk++; tk++; |
#define ONE_STEP2 if ( zzz = *pk ) { DMAR(zzz,a,*tk,md,*tk) } pk++; tk++; |
|
|
int generic_gauss_elim_mod(mat,row,col,md,colstat) |
int generic_gauss_elim_mod(mat0,row,col,md,colstat) |
int **mat; |
int **mat0; |
int row,col,md; |
int row,col,md; |
int *colstat; |
int *colstat; |
{ |
{ |
int i,j,k,l,inv,a,rank,zzz; |
int i,j,k,l,inv,a,rank,zzz; |
int *t,*pivot,*pk,*tk; |
unsigned int *t,*pivot,*pk,*tk; |
|
unsigned int **mat; |
|
|
|
mat = (unsigned int **)mat0; |
for ( rank = 0, j = 0; j < col; j++ ) { |
for ( rank = 0, j = 0; j < col; j++ ) { |
for ( i = rank; i < row && !mat[i][j]; i++ ); |
for ( i = rank; i < row; i++ ) |
|
mat[i][j] %= md; |
|
for ( i = rank; i < row; i++ ) |
|
if ( mat[i][j] ) |
|
break; |
if ( i == row ) { |
if ( i == row ) { |
colstat[j] = 0; |
colstat[j] = 0; |
continue; |
continue; |
|
|
inv = invm(pivot[j],md); |
inv = invm(pivot[j],md); |
for ( k = j, pk = pivot+k; k < col; k++, pk++ ) |
for ( k = j, pk = pivot+k; k < col; k++, pk++ ) |
if ( *pk ) { |
if ( *pk ) { |
|
if ( *pk >= md ) |
|
*pk %= md; |
DMAR(*pk,inv,0,md,*pk) |
DMAR(*pk,inv,0,md,*pk) |
} |
} |
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] ) |
a = md - a; pk = pivot+j; tk = t+j; |
red_by_vect(md,t+j,pivot+j,md-a,col-j); |
k = col-j; |
|
for ( ; k >= 64; k -= 64 ) { |
|
ONE_STEP2 ONE_STEP2 ONE_STEP2 ONE_STEP2 |
|
ONE_STEP2 ONE_STEP2 ONE_STEP2 ONE_STEP2 |
|
ONE_STEP2 ONE_STEP2 ONE_STEP2 ONE_STEP2 |
|
ONE_STEP2 ONE_STEP2 ONE_STEP2 ONE_STEP2 |
|
ONE_STEP2 ONE_STEP2 ONE_STEP2 ONE_STEP2 |
|
ONE_STEP2 ONE_STEP2 ONE_STEP2 ONE_STEP2 |
|
ONE_STEP2 ONE_STEP2 ONE_STEP2 ONE_STEP2 |
|
ONE_STEP2 ONE_STEP2 ONE_STEP2 ONE_STEP2 |
|
ONE_STEP2 ONE_STEP2 ONE_STEP2 ONE_STEP2 |
|
ONE_STEP2 ONE_STEP2 ONE_STEP2 ONE_STEP2 |
|
ONE_STEP2 ONE_STEP2 ONE_STEP2 ONE_STEP2 |
|
ONE_STEP2 ONE_STEP2 ONE_STEP2 ONE_STEP2 |
|
ONE_STEP2 ONE_STEP2 ONE_STEP2 ONE_STEP2 |
|
ONE_STEP2 ONE_STEP2 ONE_STEP2 ONE_STEP2 |
|
ONE_STEP2 ONE_STEP2 ONE_STEP2 ONE_STEP2 |
|
ONE_STEP2 ONE_STEP2 ONE_STEP2 ONE_STEP2 |
|
} |
|
for ( ; k >= 0; k -- ) { |
|
if ( zzz = *pk ) { DMAR(zzz,a,*tk,md,*tk) } pk++; tk++; |
|
} |
|
} |
|
} |
} |
rank++; |
rank++; |
} |
} |
|
|
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] ) { |
t[j] %= md; |
a = md-a; pk = pivot+j; tk = t+j; |
if ( a = t[j] ) |
k = col-j; |
red_by_vect(md,t+j,pivot+j,md-a,col-j); |
for ( ; k >= 64; k -= 64 ) { |
|
ONE_STEP2 ONE_STEP2 ONE_STEP2 ONE_STEP2 |
|
ONE_STEP2 ONE_STEP2 ONE_STEP2 ONE_STEP2 |
|
ONE_STEP2 ONE_STEP2 ONE_STEP2 ONE_STEP2 |
|
ONE_STEP2 ONE_STEP2 ONE_STEP2 ONE_STEP2 |
|
ONE_STEP2 ONE_STEP2 ONE_STEP2 ONE_STEP2 |
|
ONE_STEP2 ONE_STEP2 ONE_STEP2 ONE_STEP2 |
|
ONE_STEP2 ONE_STEP2 ONE_STEP2 ONE_STEP2 |
|
ONE_STEP2 ONE_STEP2 ONE_STEP2 ONE_STEP2 |
|
ONE_STEP2 ONE_STEP2 ONE_STEP2 ONE_STEP2 |
|
ONE_STEP2 ONE_STEP2 ONE_STEP2 ONE_STEP2 |
|
ONE_STEP2 ONE_STEP2 ONE_STEP2 ONE_STEP2 |
|
ONE_STEP2 ONE_STEP2 ONE_STEP2 ONE_STEP2 |
|
ONE_STEP2 ONE_STEP2 ONE_STEP2 ONE_STEP2 |
|
ONE_STEP2 ONE_STEP2 ONE_STEP2 ONE_STEP2 |
|
ONE_STEP2 ONE_STEP2 ONE_STEP2 ONE_STEP2 |
|
ONE_STEP2 ONE_STEP2 ONE_STEP2 ONE_STEP2 |
|
} |
|
for ( ; k >= 0; k -- ) { |
|
if ( zzz = *pk ) { DMAR(zzz,a,*tk,md,*tk) } pk++; tk++; |
|
} |
|
} |
|
} |
} |
l--; |
l--; |
|
} |
|
for ( j = 0, l = 0; l < rank; j++ ) |
|
if ( colstat[j] ) { |
|
t = mat[l]; |
|
for ( k = j; k < col; k++ ) |
|
if ( t[k] >= md ) |
|
t[k] %= md; |
|
l++; |
} |
} |
return rank; |
return rank; |
} |
} |