OpenXM_contrib2/asir2000/builtin/array.c - diff

Return to array.c CVS log

Up to [local] / OpenXM_contrib2 / asir2000 / builtin

Diff for /OpenXM_contrib2/asir2000/builtin/array.c between version 1.17 and 1.29

version 1.17, 2001/09/10 05:55:13

version 1.29, 2003/06/09 16:18:09

Line 45

* DEVELOPER SHALL HAVE NO LIABILITY IN CONNECTION WITH THE USE,

* PERFORMANCE OR NON-PERFORMANCE OF THE SOFTWARE.

* $OpenXM$

* $OpenXM: OpenXM_contrib2/asir2000/builtin/array.c,v 1.28 2003/05/29 16:44:59 saito Exp $

#include "ca.h"

#include "base.h"

Line 59

extern int DP_Print; /* XXX */

void inner_product_mat_int_mod(Q **,int **,int,int,int,Q *);

void solve_by_lu_mod(int **,int,int,int **,int);

void solve_by_lu_gfmmat(GFMMAT,unsigned int,unsigned int *,unsigned int *);

int lu_gfmmat(GFMMAT,unsigned int,int *);

void mat_to_gfmmat(MAT,unsigned int,GFMMAT *);

int generic_gauss_elim_mod(int **,int,int,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_mod1(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 **);

void Pnewvect(), Pnewmat(), Psepvect(), Psize(), Pdet(), Pleqm(), Pleqm1(), Pgeninvm();

void Pinvmat();

void Pnewbytearray();

void Pgeneric_gauss_elim();

void Pgeneric_gauss_elim_mod();

void Pmat_to_gfmmat(),Plu_gfmmat(),Psolve_by_lu_gfmmat();

void Pgeninvm_swap(), Premainder(), Psremainder(), Pvtol();

void Pgeninv_sf_swap();

void sepvect();

void Pmulmat_gf2n();

void Pbconvmat_gf2n();

Line 92 void Pirredpoly_up2();

Line 82 void Pirredpoly_up2();

void Pnbpoly_up2();

void Pqsort();

void Pexponent_vector();

void Pmat_swap_row_destructive();

void Pmat_swap_col_destructive();

void Pvect();

void Pmat();

void Pmatc();

struct ftab array_tab[] = {

{"solve_by_lu_gfmmat",Psolve_by_lu_gfmmat,4},

{"lu_gfmmat",Plu_gfmmat,2},

{"mat_to_gfmmat",Pmat_to_gfmmat,2},

{"generic_gauss_elim",Pgeneric_gauss_elim,1},

{"generic_gauss_elim_mod",Pgeneric_gauss_elim_mod,2},

{"newvect",Pnewvect,-2},

{"vect",Pvect,-99999999},

{"vector",Pnewvect,-2},

{"exponent_vector",Pexponent_vector,-99999999},

{"newmat",Pnewmat,-3},

{"matrix",Pnewmat,-3},

{"mat",Pmat,-99999999},

{"matr",Pmat,-99999999},

{"matc",Pmatc,-99999999},

{"newbytearray",Pnewbytearray,-2},

{"sepmat_destructive",Psepmat_destructive,2},

{"sepvect",Psepvect,2},

Line 110 struct ftab array_tab[] = {

{"vtol",Pvtol,1},

{"size",Psize,1},

{"det",Pdet,-2},

{"invmat",Pinvmat,-2},

{"leqm",Pleqm,2},

{"leqm1",Pleqm1,2},

{"geninvm",Pgeninvm,2},

{"geninvm_swap",Pgeninvm_swap,2},

{"geninv_sf_swap",Pgeninv_sf_swap,1},

{"remainder",Premainder,2},

{"sremainder",Psremainder,2},

{"mulmat_gf2n",Pmulmat_gf2n,1},

Line 124 struct ftab array_tab[] = {

Line 126 struct ftab array_tab[] = {

{"x962_irredpoly_up2",Px962_irredpoly_up2,2},

{"irredpoly_up2",Pirredpoly_up2,2},

{"nbpoly_up2",Pnbpoly_up2,2},

{"mat_swap_row_destructive",Pmat_swap_row_destructive,3},

{"mat_swap_col_destructive",Pmat_swap_col_destructive,3},

{0,0,0},

};

int comp_obj(a,b)

int comp_obj(Obj *a,Obj *b)

Obj *a,*b;

{

return arf_comp(CO,*a,*b);

}

Line 136 Obj *a,*b;

Line 139 Obj *a,*b;

static FUNC generic_comp_obj_func;

static NODE generic_comp_obj_arg;

int generic_comp_obj(a,b)

int generic_comp_obj(Obj *a,Obj *b)

Obj *a,*b;

{

Q r;

Line 151 Obj *a,*b;

Line 153 Obj *a,*b;

}

void Pqsort(arg,rp)

void Pqsort(NODE arg,VECT *rp)

NODE arg;

VECT *rp;

{

VECT vect;

char buf[BUFSIZ];

char *fname;

NODE n;

P p;

V v;

Line 180 VECT *rp;

Line 178 VECT *rp;

*rp = vect;

}

void PNBmul_gf2n(arg,rp)

void PNBmul_gf2n(NODE arg,GF2N *rp)

NODE arg;

GF2N *rp;

{

GF2N a,b;

GF2MAT mat;

Line 219 GF2N *rp;

Line 215 GF2N *rp;

}

void Pmul_vect_mat_gf2n(arg,rp)

void Pmul_vect_mat_gf2n(NODE arg,GF2N *rp)

NODE arg;

GF2N *rp;

{

GF2N a;

GF2MAT mat;

Line 252 GF2N *rp;

Line 246 GF2N *rp;

}

void Pbconvmat_gf2n(arg,rp)

void Pbconvmat_gf2n(NODE arg,LIST *rp)

NODE arg;

LIST *rp;

{

P p0,p1;

int to;

Line 274 LIST *rp;

Line 266 LIST *rp;

MKLIST(*rp,n0);

}

void Pmulmat_gf2n(arg,rp)

void Pmulmat_gf2n(NODE arg,GF2MAT *rp)

NODE arg;

GF2MAT *rp;

{

GF2MAT m;

Line 285 GF2MAT *rp;

Line 275 GF2MAT *rp;

*rp = m;

}

void Psepmat_destructive(arg,rp)

void Psepmat_destructive(NODE arg,LIST *rp)

NODE arg;

LIST *rp;

{

MAT mat,mat1;

int i,j,row,col;

Line 318 LIST *rp;

Line 306 LIST *rp;

MKLIST(*rp,n0);

}

void Psepvect(arg,rp)

void Psepvect(NODE arg,VECT *rp)

NODE arg;

VECT *rp;

{

sepvect((VECT)ARG0(arg),QTOS((Q)ARG1(arg)),rp);

}

void sepvect(v,d,rp)

void sepvect(VECT v,int d,VECT *rp)

VECT v;

int d;

VECT *rp;

{

int i,j,k,n,q,q1,r;

pointer *pv,*pw,*pu;

Line 352 VECT *rp;

Line 335 VECT *rp;

}

void Pnewvect(arg,rp)

void Pnewvect(NODE arg,VECT *rp)

NODE arg;

VECT *rp;

{

int len,i,r;

VECT vect;

Line 381 VECT *rp;

Line 362 VECT *rp;

*rp = vect;

}

void Pexponent_vector(arg,rp)

void Pvect(NODE arg,VECT *rp) {

NODE arg;

int len,i,r;

DP *rp;

VECT vect;

pointer *vb;

NODE tn;

if ( !arg ) {

*rp =0;

return;

}

for (len = 0, tn = arg; tn; tn = NEXT(tn), len++);

if ( len == 1 ) {

if ( ARG0(arg) != 0 ) {

switch ( OID(ARG0(arg)) ) {

case O_VECT:

*rp = ARG0(arg);

return;

case O_LIST:

for ( len = 0, tn = ARG0(arg); tn; tn = NEXT(tn), len++ );

MKVECT(vect,len-1);

for ( i = 0, tn = BDY((LIST)ARG0(arg)), vb =BDY(vect);

tn; i++, tn = NEXT(tn) )

vb[i] = (pointer)BDY(tn);

*rp=vect;

return;

}

MKVECT(vect,len);

for ( i = 0, tn = arg, vb = BDY(vect); tn; i++, tn = NEXT(tn) )

vb[i] = (pointer)BDY(tn);

*rp = vect;

}

void Pexponent_vector(NODE arg,DP *rp)

{

nodetod(arg,rp);

}

void Pnewbytearray(arg,rp)

void Pnewbytearray(NODE arg,BYTEARRAY *rp)

NODE arg;

BYTEARRAY *rp;

{

int len,i,r;

BYTEARRAY array;

Line 432 BYTEARRAY *rp;

Line 444 BYTEARRAY *rp;

*rp = array;

}

void Pnewmat(arg,rp)

void Pnewmat(NODE arg,MAT *rp)

NODE arg;

MAT *rp;

{

int row,col;

int i,j,r,c;

Line 469 MAT *rp;

Line 479 MAT *rp;

*rp = m;

}

void Pvtol(arg,rp)

void Pmat(NODE arg, MAT *rp)

NODE arg;

LIST *rp;

{

int row,col;

int i;

MAT m;

pointer **mb;

pointer *ent;

NODE tn, sn;

VECT v;

if ( !arg ) {

*rp =0;

return;

}

for (row = 0, tn = arg; tn; tn = NEXT(tn), row++);

if ( OID(ARG0(arg)) == O_VECT ) {

v = ARG0(arg);

col = v->len;

} else if ( OID(ARG0(arg)) == O_LIST ) {

for (col = 0, tn = BDY((LIST)ARG0(arg)); tn ; tn = NEXT(tn), col++);

}

MKMAT(m,row,col);

for (row = 0, tn = arg, mb = BDY(m); tn; tn = NEXT(tn), row++) {

if ( BDY(tn) == 0 ) {

error("mat : invalid argument");

} else if ( OID(BDY(tn)) == O_VECT ) {

v = tn->body;

ent = BDY(v);

for (i = 0; i < v->len; i++ ) mb[row][i] = (Obj)ent[i];

} else if ( OID(BDY(tn)) == O_LIST ) {

for (col = 0, sn = BDY((LIST)BDY(tn)); sn; col++, sn = NEXT(sn) )

mb[row][col] = (pointer)BDY(sn);

} else {

error("mat : invalid argument");

}

*rp = m;

}

void Pmatc(NODE arg, MAT *rp)

{

int row,col;

int i;

MAT m;

pointer **mb;

pointer *ent;

NODE tn, sn;

VECT v;

if ( !arg ) {

*rp =0;

return;

}

for (col = 0, tn = arg; tn; tn = NEXT(tn), col++);

if ( OID(ARG0(arg)) == O_VECT ) {

v = ARG0(arg);

row = v->len;

} else if ( OID(ARG0(arg)) == O_LIST ) {

for (row = 0, tn = BDY((LIST)ARG0(arg)); tn ; tn = NEXT(tn), row++);

}

MKMAT(m,row,col);

for (col = 0, tn = arg, mb = BDY(m); tn; tn = NEXT(tn), col++) {

if ( BDY(tn) == 0 ) {

error("matc : invalid argument");

} else if ( OID(BDY(tn)) == O_VECT ) {

v = tn->body;

ent = BDY(v);

for (i = 0; i < v->len; i++ ) mb[i][col] = (Obj)ent[i];

} else if ( OID(BDY(tn)) == O_LIST ) {

for (row = 0, sn = BDY((LIST)BDY(tn)); sn; row++, sn = NEXT(sn) )

mb[row][col] = (pointer)BDY(sn);

} else {

error("matc : invalid argument");

}

*rp = m;

}

void Pvtol(NODE arg,LIST *rp)

{

NODE n,n1;

VECT v;

pointer *a;

Line 486 LIST *rp;

Line 576 LIST *rp;

MKLIST(*rp,n);

}

void Premainder(arg,rp)

void Premainder(NODE arg,Obj *rp)

NODE arg;

Obj *rp;

{

Obj a;

VECT v,w;

Line 536 Obj *rp;

Line 624 Obj *rp;

}

void Psremainder(arg,rp)

void Psremainder(NODE arg,Obj *rp)

NODE arg;

Obj *rp;

{

Obj a;

VECT v,w;

Line 587 Obj *rp;

Line 673 Obj *rp;

}

void Psize(arg,rp)

void Psize(NODE arg,LIST *rp)

NODE arg;

LIST *rp;

{

int n,m;

Line 615 LIST *rp;

Line 699 LIST *rp;

MKLIST(*rp,t);

}

void Pdet(arg,rp)

void Pdet(NODE arg,P *rp)

NODE arg;

P *rp;

{

MAT m;

int n,i,j,mod;

Line 641 P *rp;

Line 723 P *rp;

}

void Pinvmat(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

A[I] <-> A[I][0]*x_0+A[I][1]*x_1+...

Line 652 P *rp;

Line 767 P *rp;

B[I] <-> x_{R[I]}+B[I][0]x_{C[0]}+B[I][1]x_{C[1]}+...

void Pgeneric_gauss_elim_mod(arg,rp)

void Pgeneric_gauss_elim(NODE arg,LIST *rp)

NODE arg;

LIST *rp;

{

NODE n0;

MAT m,nm;

int *ri,*ci;

VECT rind,cind;

Q dn,q;

int i,j,k,l,row,col,t,rank;

asir_assert(ARG0(arg),O_MAT,"generic_gauss_elim");

m = (MAT)ARG0(arg);

row = m->row; col = m->col;

rank = generic_gauss_elim(m,&nm,&dn,&ri,&ci);

t = col-rank;

MKVECT(rind,rank);

MKVECT(cind,t);

for ( i = 0; i < rank; i++ ) {

STOQ(ri[i],q);

BDY(rind)[i] = (pointer)q;

}

for ( i = 0; i < t; i++ ) {

STOQ(ci[i],q);

BDY(cind)[i] = (pointer)q;

}

n0 = mknode(4,nm,dn,rind,cind);

MKLIST(*rp,n0);

}

input : a row x col matrix A

A[I] <-> A[I][0]*x_0+A[I][1]*x_1+...

output : [B,R,C]

B : a rank(A) x col-rank(A) matrix

R : a vector of length rank(A)

C : a vector of length col-rank(A)

B[I] <-> x_{R[I]}+B[I][0]x_{C[0]}+B[I][1]x_{C[1]}+...

void Pgeneric_gauss_elim_mod(NODE arg,LIST *rp)

{

NODE n0;

MAT m,mat;

VECT rind,cind;

Q **tmat;

Line 664 LIST *rp;

Line 816 LIST *rp;

Q *rib,*cib;

int *colstat;

Q q;

int md,i,j,k,l,row,col,t,n,rank;

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");

Line 704 LIST *rp;

Line 856 LIST *rp;

MKLIST(*rp,n0);

}

void Pleqm(arg,rp)

void Pleqm(NODE arg,VECT *rp)

NODE arg;

VECT *rp;

{

MAT m;

VECT vect;

Line 745 VECT *rp;

Line 895 VECT *rp;

}

int gauss_elim_mod(mat,row,col,md)

int gauss_elim_mod(int **mat,int row,int col,int md)

int **mat;

int row,col,md;

{

int i,j,k,inv,a,n;

int *t,*pivot;

Line 786 int row,col,md;

Line 934 int row,col,md;

struct oEGT eg_mod,eg_elim,eg_elim1,eg_elim2,eg_chrem,eg_gschk,eg_intrat,eg_symb;

int generic_gauss_elim(mat,nm,dn,rindp,cindp)

int generic_gauss_elim(MAT mat,MAT *nm,Q *dn,int **rindp,int **cindp)

MAT mat;

MAT *nm;

Q *dn;

int **rindp,**cindp;

{

int **wmat;

Q **bmat;

Line 942 RESET:

Line 1086 RESET:

}

int generic_gauss_elim_hensel(mat,nmmat,dn,rindp,cindp)

int generic_gauss_elim_hensel(MAT mat,MAT *nmmat,Q *dn,int **rindp,int **cindp)

MAT mat;

MAT *nmmat;

Q *dn;

int **rindp,**cindp;

{

MAT bmat,xmat;

Q **a0,**a,**b,**x,**nm;

Line 980 int **rindp,**cindp;

Line 1120 int **rindp,**cindp;

} else

wi[j] = 0;

rank = find_lhs_and_lu_mod(w,row,col,md,&rinfo,&cinfo);

rank = find_lhs_and_lu_mod((unsigned int **)w,row,col,md,&rinfo,&cinfo);

a = (Q **)almat_pointer(rank,rank); /* lhs mat */

MKMAT(bmat,rank,col-rank); b = (Q **)bmat->body; /* lhs mat */

for ( j = li = ri = 0; j < col; j++ )

Line 1070 int **rindp,**cindp;

Line 1210 int **rindp,**cindp;

int f4_nocheck;

int gensolve_check(mat,nm,dn,rind,cind)

int gensolve_check(MAT mat,MAT nm,Q dn,int *rind,int *cind)

MAT mat,nm;

Q dn;

int *rind,*cind;

{

int row,col,rank,clen,i,j,k,l;

Q s,t,u;

Q s,t;

Q *w;

Q *mati,*nmk;

Line 1121 int *rind,*cind;

Line 1258 int *rind,*cind;

/* assuming 0 < c < m */

int inttorat(c,m,b,sgnp,nmp,dnp)

int inttorat(N c,N m,N b,int *sgnp,N *nmp,N *dnp)

N c,m,b;

int *sgnp;

N *nmp,*dnp;

{

Q qq,t,u1,v1,r1,nm;

Q qq,t,u1,v1,r1;

N q,r,u2,v2,r2;

N q,u2,v2,r2;

u1 = 0; v1 = ONE; u2 = m; v2 = c;

while ( cmpn(v2,b) >= 0 ) {

Line 1146 N *nmp,*dnp;

Line 1280 N *nmp,*dnp;

/* mat->body = N ** */

int intmtoratm(mat,md,nm,dn)

int intmtoratm(MAT mat,N md,MAT nm,Q *dn)

MAT mat;

N md;

MAT nm;

Q *dn;

{

N t,s,b;

Q bound,dn0,dn1,nm1,q,tq;

Q dn0,dn1,nm1,q;

int i,j,k,l,row,col;

Q **rmat;

N **tmat;

Line 1203 Q *dn;

Line 1333 Q *dn;

/* mat->body = Q ** */

int intmtoratm_q(mat,md,nm,dn)

int intmtoratm_q(MAT mat,N md,MAT nm,Q *dn)

MAT mat;

N md;

MAT nm;

Q *dn;

{

N t,s,b;

Q bound,dn0,dn1,nm1,q,tq;

Q dn0,dn1,nm1,q;

int i,j,k,l,row,col;

Q **rmat;

Q **tmat;

Line 1262 Q *dn;

Line 1388 Q *dn;

#define ONE_STEP1 if ( zzz = *s ) { DMAR(zzz,hc,*tj,md,*tj) } tj++; s++;

void reduce_reducers_mod(mat,row,col,md)

void reduce_reducers_mod(int **mat,int row,int col,int md)

int **mat;

int row,col;

int md;

{

int i,j,k,l,hc,zzz;

int *t,*s,*tj,*ind;

Line 1320 int md;

Line 1443 int md;

2. reduce spolys by the reduced reducers

void pre_reduce_mod(mat,row,col,nred,md)

void pre_reduce_mod(int **mat,int row,int col,int nred,int md)

int **mat;

int row,col,nred;

int md;

{

int i,j,k,l,hc,inv;

int *t,*s,*tk,*ind;

Line 1373 int md;

Line 1493 int md;

mat[0] < mat[1] < ... < mat[nred-1] w.r.t the term order

void reduce_sp_by_red_mod(sp,redmat,ind,nred,col,md)

void reduce_sp_by_red_mod(int *sp,int **redmat,int *ind,int nred,int col,int md)

int *sp,**redmat;

int *ind;

int nred,col;

int md;

{

int i,j,k,hc,zzz;

int *t,*s,*tj;

int *s,*tj;

/* reduce the spolys by redmat */

for ( i = nred-1; i >= 0; i-- ) {

Line 1399 int md;

Line 1515 int md;

}

rlist : reducers list

mat[i] : compressed reducers (i=0,...,nred-1)

ht(BDY(rlist)) < ht(BDY(NEXT(rlist)) < ... w.r.t. the term order

mat[0] < mat[1] < ... < mat[nred-1] w.r.t the term order

void reduce_reducers_mod_compress(rlist,nred,at,col,md,redmatp,indredp)

void red_by_compress(int m,unsigned int *p,unsigned int *r,

NODE rlist;

unsigned int *ri,unsigned int hc,int len)

int nred;

DL *at;

int col,md;

CDP **redmatp;

int **indredp;

{

CDP *redmat;

unsigned int up,lo;

CDP t;

unsigned int dmy;

int *indred,*w;

unsigned int *pj;

int i,k;

NODE r;

*redmatp = redmat = (CDP *)CALLOC(nred,sizeof(CDP));

p[*ri] = 0; r++; ri++;

*indredp = indred = (int *)CALLOC(nred,sizeof(int));

for ( len--; len; len--, r++, ri++ ) {

w = (int *)CALLOC(col,sizeof(int));

pj = p+ *ri;

DMA(*r,hc,*pj,up,lo);

if ( up ) {

DSAB(m,up,lo,dmy,*pj);

} else

*pj = lo;

}

_dpmod_to_vect_compress(BDY(rlist),at,&redmat[0]);

/* p -= hc*r */

indred[0] = redmat[0]->body[0].index;

for ( i = 1, r = NEXT(rlist); i < nred; i++, r = NEXT(r) ) {

void red_by_vect(int m,unsigned int *p,unsigned int *r,unsigned int hc,int len)

bzero(w,col*sizeof(int));

{

_dpmod_to_vect(BDY(r),at,w);

reduce_sp_by_red_mod_compress(w,redmat,indred,i,col,md);

unsigned int dmy;

compress_vect(w,col,&redmat[i]);

indred[i] = redmat[i]->body[0].index;

*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;

mat[i] : compressed reducers (i=0,...,nred-1)

mat[0] < mat[1] < ... < mat[nred-1] w.r.t the term order

void reduce_sp_by_red_mod_compress (sp,redmat,ind,nred,col,md)

void reduce_sp_by_red_mod_compress (int *sp,CDP *redmat,int *ind,

int *sp;

int nred,int col,int md)

CDP *redmat;

int *ind;

int nred,col;

int md;

{

int i,j,k,hc,c,len;

int i,len;

int *tj;

CDP ri;

unsigned int hc;

unsigned int *usp;

usp = (unsigned int *)sp;

/* reduce the spolys by redmat */

for ( i = nred-1; i >= 0; i-- ) {

/* reduce sp by redmat[i] */

if ( hc = sp[ind[i]] ) {

usp[ind[i]] %= md;

if ( hc = usp[ind[i]] ) {

/* sp = sp-hc*redmat[i] */

hc = md-hc;

ri = redmat[i];

len = ri->len;

for ( k = 0; k < len; k++ ) {

red_by_compress(md,usp,psca[ri->psindex],ri->body,hc,len);

j = ri->body[k].index;

c = ri->body[k].c;

tj = sp+j;

#if 1

DMAR(c,hc,*tj,md,*tj);

#else

*tj = ((hc*c)+(*tj))%md;

#endif

}

for ( i = 0; i < col; i++ )

if ( usp[i] >= (unsigned int)md )

usp[i] %= md;

}

#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(int **mat0,int row,int col,int md,int *colstat)

int **mat;

int row,col,md;

int *colstat;

{

int i,j,k,l,inv,a,rank,zzz;

int i,j,k,l,inv,a,rank;

int *t,*pivot,*pk,*tk;

unsigned int *t,*pivot,*pk;

unsigned int **mat;

mat = (unsigned int **)mat0;

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 ) {

colstat[j] = 0;

continue;

Line 1495 int *colstat;

Line 1610 int *colstat;

inv = invm(pivot[j],md);

for ( k = j, pk = pivot+k; k < col; k++, pk++ )

if ( *pk ) {

if ( *pk >= (unsigned int)md )

*pk %= md;

DMAR(*pk,inv,0,md,*pk)

}

for ( i = rank+1; i < row; 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

}

for ( ; k > 0; k -- ) {

if ( zzz = *pk ) { DMAR(zzz,a,*tk,md,*tk) } pk++; tk++;

}

rank++;

}

Line 1532 int *colstat;

Line 1626 int *colstat;

pivot = mat[l];

for ( i = 0; i < l; 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

}

for ( ; k > 0; k -- ) {

if ( zzz = *pk ) { DMAR(zzz,a,*tk,md,*tk) } pk++; tk++;

}

l--;

}

for ( j = 0, l = 0; l < rank; j++ )

if ( colstat[j] ) {

t = mat[l];

for ( k = j; k < col; k++ )

if ( t[k] >= (unsigned int)md )

t[k] %= md;

l++;

}

return rank;

}

/* LU decomposition; a[i][i] = 1/U[i][i] */

int lu_gfmmat(mat,md,perm)

int lu_gfmmat(GFMMAT mat,unsigned int md,int *perm)

GFMMAT mat;

unsigned int md;

int *perm;

{

int row,col;

int i,j,k,l;

int i,j,k;

unsigned int *t,*pivot;

unsigned int **a;

unsigned int inv,m;

Line 1623 int *perm;

Line 1700 int *perm;

cinfo[j]=1 <=> j-th column is contained in the LU decomp.

int find_lhs_and_lu_mod(a,row,col,md,rinfo,cinfo)

int find_lhs_and_lu_mod(unsigned int **a,int row,int col,

unsigned int **a;

unsigned int md,int **rinfo,int **cinfo)

unsigned int md;

int **rinfo,**cinfo;

{

int i,j,k,l,d;

int i,j,k,d;

int *rp,*cp;

unsigned int *t,*pivot;

unsigned int inv,m;

Line 1676 int **rinfo,**cinfo;

Line 1751 int **rinfo,**cinfo;

b = a^(-1)b

void solve_by_lu_mod(a,n,md,b,l)

void solve_by_lu_mod(int **a,int n,int md,int **b,int l)

int **a;

int n;

int md;

int **b;

int l;

{

unsigned int *y,*c;

int i,j,k;

Line 1719 int l;

Line 1789 int l;

}

void Pleqm1(arg,rp)

void Pleqm1(NODE arg,VECT *rp)

NODE arg;

VECT *rp;

{

MAT m;

VECT vect;

Line 1760 VECT *rp;

Line 1828 VECT *rp;

}

gauss_elim_mod1(mat,row,col,md)

int gauss_elim_mod1(int **mat,int row,int col,int md)

int **mat;

int row,col,md;

{

int i,j,k,inv,a,n;

int *t,*pivot;

Line 1799 int row,col,md;

Line 1865 int row,col,md;

return -1;

}

void Pgeninvm(arg,rp)

void Pgeninvm(NODE arg,LIST *rp)

NODE arg;

LIST *rp;

{

MAT m;

pointer **mat;

Line 1835 LIST *rp;

Line 1899 LIST *rp;

MKMAT(mat1,col,row); MKMAT(mat2,row-col,row);

for ( i = 0, tmat = (Q **)mat1->body; i < col; i++ )

for ( j = 0; j < row; j++ )

STOQ(wmat[i][j+col],tmat[i][j]);

UTOQ(wmat[i][j+col],tmat[i][j]);

for ( tmat = (Q **)mat2->body; i < row; i++ )

for ( j = 0; j < row; j++ )

STOQ(wmat[i][j+col],tmat[i-col][j]);

UTOQ(wmat[i][j+col],tmat[i-col][j]);

MKNODE(node2,mat2,0); MKNODE(node1,mat1,node2); MKLIST(*rp,node1);

}

int gauss_elim_geninv_mod(mat,row,col,md)

int gauss_elim_geninv_mod(unsigned int **mat,int row,int col,int md)

unsigned int **mat;

int row,col,md;

{

int i,j,k,inv,a,n,m;

unsigned int *t,*pivot;

Line 1881 int row,col,md;

Line 1943 int row,col,md;

return 0;

}

void Psolve_by_lu_gfmmat(arg,rp)

void Psolve_by_lu_gfmmat(NODE arg,VECT *rp)

NODE arg;

VECT *rp;

{

GFMMAT lu;

Q *perm,*rhs,*v;

Line 1904 VECT *rp;

Line 1964 VECT *rp;

solve_by_lu_gfmmat(lu,md,b,sol);

MKVECT(r,n);

for ( i = 0, v = (Q *)r->body; i < n; i++ )

STOQ(sol[i],v[i]);

UTOQ(sol[i],v[i]);

*rp = r;

}

void solve_by_lu_gfmmat(lu,md,b,x)

void solve_by_lu_gfmmat(GFMMAT lu,unsigned int md,

GFMMAT lu;

unsigned int *b,unsigned int *x)

unsigned int md;

unsigned int *b;

unsigned int *x;

{

int n;

unsigned int **a;

Line 1944 unsigned int *x;

Line 2001 unsigned int *x;

}

void Plu_gfmmat(arg,rp)

void Plu_gfmmat(NODE arg,LIST *rp)

NODE arg;

LIST *rp;

{

MAT m;

GFMMAT mm;

Line 1976 LIST *rp;

Line 2031 LIST *rp;

MKLIST(*rp,n0);

}

void Pmat_to_gfmmat(arg,rp)

void Pmat_to_gfmmat(NODE arg,GFMMAT *rp)

NODE arg;

GFMMAT *rp;

{

MAT m;

unsigned int md;

Line 1989 GFMMAT *rp;

Line 2042 GFMMAT *rp;

mat_to_gfmmat(m,md,rp);

}

void mat_to_gfmmat(m,md,rp)

void mat_to_gfmmat(MAT m,unsigned int md,GFMMAT *rp)

MAT m;

unsigned int md;

GFMMAT *rp;

{

unsigned int **wmat;

unsigned int t;

Line 2118 int **indexp;

Line 2168 int **indexp;

return 0;

}

void Pgeninv_sf_swap(NODE arg,LIST *rp)

{

MAT m;

GFS **mat,**tmat;

Q *tvect;

GFS q;

int **wmat,**invmat;

int *index;

unsigned int t;

int i,j,row,col,status;

MAT mat1;

VECT vect1;

NODE node1,node2;

asir_assert(ARG0(arg),O_MAT,"geninv_sf_swap");

m = (MAT)ARG0(arg);

row = m->row; col = m->col; mat = (GFS **)m->body;

wmat = (int **)almat(row,col+row);

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] )

wmat[i][j] = FTOIF(CONT(q));

wmat[i][col+i] = _onesf();

}

status = gauss_elim_geninv_sf_swap(wmat,row,col,&invmat,&index);

if ( status > 0 )

*rp = 0;

else {

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] ) {

MKGFS(IFTOF(t),tmat[i][j]);

}

MKVECT(vect1,row);

for ( i = 0, tvect = (Q *)vect1->body; i < row; i++ )

STOQ(index[i],tvect[i]);

MKNODE(node2,vect1,0); MKNODE(node1,mat1,node2); MKLIST(*rp,node1);

}

int gauss_elim_geninv_sf_swap(int **mat,int row,int col,

int ***invmatp,int **indexp)

{

int i,j,k,inv,a,n,m,u;

int *t,*pivot,*s;

int *index;

int **invmat;

n = col; m = row+col;

*indexp = index = (int *)MALLOC_ATOMIC(row*sizeof(int));

for ( i = 0; i < row; i++ )

index[i] = i;

for ( j = 0; j < n; j++ ) {

for ( i = j; i < row && !mat[i][j]; i++ );

if ( i == row ) {

*indexp = 0; *invmatp = 0; return 1;

}

if ( i != j ) {

t = mat[i]; mat[i] = mat[j]; mat[j] = t;

k = index[i]; index[i] = index[j]; index[j] = k;

}

pivot = mat[j];

inv = _invsf(pivot[j]);

for ( k = j; k < m; k++ )

if ( pivot[k] )

pivot[k] = _mulsf(pivot[k],inv);

for ( i = j+1; i < row; i++ ) {

t = mat[i];

if ( a = t[j] )

for ( k = j, a = _chsgnsf(a); k < m; k++ )

if ( pivot[k] ) {

u = _mulsf(pivot[k],a);

t[k] = _addsf(u,t[k]);

}

for ( j = n-1; j >= 0; j-- ) {

pivot = mat[j];

for ( i = j-1; i >= 0; i-- ) {

t = mat[i];

if ( a = t[j] )

for ( k = j, a = _chsgnsf(a); k < m; k++ )

if ( pivot[k] ) {

u = _mulsf(pivot[k],a);

t[k] = _addsf(u,t[k]);

}

*invmatp = invmat = (int **)almat(col,col);

for ( i = 0; i < col; i++ )

for ( j = 0, s = invmat[i], t = mat[i]; j < col; j++ )

s[j] = t[col+index[j]];

return 0;

}

void _addn(N,N,N);

int _subn(N,N,N);

void _muln(N,N,N);

void inner_product_int(a,b,n,r)

void inner_product_int(Q *a,Q *b,int n,Q *r)

Q *a,*b;

int n;

Q *r;

{

int la,lb,i;

int sgn,sgn1;

Line 2178 Q *r;

Line 2322 Q *r;

/* (k,l) element of a*b where a: .x n matrix, b: n x . integer matrix */

void inner_product_mat_int_mod(a,b,n,k,l,r)

void inner_product_mat_int_mod(Q **a,int **b,int n,int k,int l,Q *r)

Q **a;

int **b;

int n,k,l;

Q *r;

{

int la,lb,i;

int sgn,sgn1;

Line 2238 Q *r;

Line 2378 Q *r;

NTOQ(sum,sgn,*r);

}

void Pmul_mat_vect_int(arg,rp)

void Pmul_mat_vect_int(NODE arg,VECT *rp)

NODE arg;

VECT *rp;

{

MAT mat;

VECT vect,r;

Line 2251 VECT *rp;

Line 2389 VECT *rp;

row = mat->row;

col = mat->col;

MKVECT(r,row);

for ( i = 0; i < row; i++ )

for ( i = 0; i < row; i++ ) {

inner_product_int(mat->body[i],vect->body,col,&r->body[i]);

inner_product_int((Q *)mat->body[i],(Q *)vect->body,col,(Q *)&r->body[i]);

}

*rp = r;

}

void Pnbpoly_up2(arg,rp)

void Pnbpoly_up2(NODE arg,GF2N *rp)

NODE arg;

GF2N *rp;

{

int m,type,ret;

UP2 r;

Line 2272 GF2N *rp;

Line 2409 GF2N *rp;

*rp = 0;

}

void Px962_irredpoly_up2(arg,rp)

void Px962_irredpoly_up2(NODE arg,GF2N *rp)

NODE arg;

GF2N *rp;

{

int m,type,ret,w;

int m,ret,w;

GF2N prev;

UP2 r;

Line 2292 GF2N *rp;

Line 2427 GF2N *rp;

bzero((char *)r->b,w*sizeof(unsigned int));

}

ret = _generate_irreducible_polynomial(r,m,type);

ret = _generate_irreducible_polynomial(r,m);

if ( ret == 0 )

MKGF2N(r,*rp);

else

*rp = 0;

}

void Pirredpoly_up2(arg,rp)

void Pirredpoly_up2(NODE arg,GF2N *rp)

NODE arg;

GF2N *rp;

{

int m,type,ret,w;

int m,ret,w;

GF2N prev;

UP2 r;

Line 2319 GF2N *rp;

Line 2452 GF2N *rp;

bzero((char *)r->b,w*sizeof(unsigned int));

}

ret = _generate_good_irreducible_polynomial(r,m,type);

ret = _generate_good_irreducible_polynomial(r,m);

if ( ret == 0 )

MKGF2N(r,*rp);

else

*rp = 0;

}

void Pmat_swap_row_destructive(NODE arg, MAT *m)

{

int i1,i2;

pointer *t;

MAT mat;

asir_assert(ARG0(arg),O_MAT,"mat_swap_row_destructive");

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));

if ( i1 < 0 || i2 < 0 || i1 >= mat->row || i2 >= mat->row )

error("mat_swap_row_destructive : Out of range");

t = mat->body[i1];

mat->body[i1] = mat->body[i2];

mat->body[i2] = t;

*m = mat;

}

void Pmat_swap_col_destructive(NODE arg, MAT *m)

{

int j1,j2,i,n;

pointer *mi;

pointer t;

MAT mat;

asir_assert(ARG0(arg),O_MAT,"mat_swap_col_destructive");

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));

if ( j1 < 0 || j2 < 0 || j1 >= mat->col || j2 >= mat->col )

error("mat_swap_col_destructive : Out of range");

n = mat->row;

for ( i = 0; i < n; i++ ) {

mi = mat->body[i];

t = mi[j1]; mi[j1] = mi[j2]; mi[j2] = t;

}

*m = mat;

}

* f = type 'type' normal polynomial of degree m if exists

* IEEE P1363 A.7.2

Line 2599 PENTA:

Line 2774 PENTA:

return 1;

}

printqmat(mat,row,col)

void printqmat(Q **mat,int row,int col)

Q **mat;

int row,col;

{

int i,j;

Line 2613 int row,col;

Line 2786 int row,col;

}

printimat(mat,row,col)

void printimat(int **mat,int row,int col)

int **mat;

int row,col;

{

int i,j;

FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>