===================================================================
RCS file: /home/cvs/OpenXM_contrib2/asir2000/builtin/array.c,v
retrieving revision 1.5
retrieving revision 1.25
diff -u -p -r1.5 -r1.25
--- OpenXM_contrib2/asir2000/builtin/array.c 2000/06/05 02:26:47 1.5
+++ OpenXM_contrib2/asir2000/builtin/array.c 2001/12/20 08:18:26 1.25
@@ -1,4 +1,52 @@
-/* $OpenXM: OpenXM_contrib2/asir2000/builtin/array.c,v 1.4 2000/05/29 08:54:44 noro Exp $ */
+/*
+ * Copyright (c) 1994-2000 FUJITSU LABORATORIES LIMITED
+ * All rights reserved.
+ *
+ * FUJITSU LABORATORIES LIMITED ("FLL") hereby grants you a limited,
+ * non-exclusive and royalty-free license to use, copy, modify and
+ * redistribute, solely for non-commercial and non-profit purposes, the
+ * computer program, "Risa/Asir" ("SOFTWARE"), subject to the terms and
+ * conditions of this Agreement. For the avoidance of doubt, you acquire
+ * only a limited right to use the SOFTWARE hereunder, and FLL or any
+ * third party developer retains all rights, including but not limited to
+ * copyrights, in and to the SOFTWARE.
+ *
+ * (1) FLL does not grant you a license in any way for commercial
+ * purposes. You may use the SOFTWARE only for non-commercial and
+ * non-profit purposes only, such as academic, research and internal
+ * business use.
+ * (2) The SOFTWARE is protected by the Copyright Law of Japan and
+ * international copyright treaties. If you make copies of the SOFTWARE,
+ * with or without modification, as permitted hereunder, you shall affix
+ * to all such copies of the SOFTWARE the above copyright notice.
+ * (3) An explicit reference to this SOFTWARE and its copyright owner
+ * shall be made on your publication or presentation in any form of the
+ * results obtained by use of the SOFTWARE.
+ * (4) In the event that you modify the SOFTWARE, you shall notify FLL by
+ * e-mail at risa-admin@sec.flab.fujitsu.co.jp of the detailed specification
+ * for such modification or the source code of the modified part of the
+ * SOFTWARE.
+ *
+ * THE SOFTWARE IS PROVIDED AS IS WITHOUT ANY WARRANTY OF ANY KIND. FLL
+ * MAKES ABSOLUTELY NO WARRANTIES, EXPRESSED, IMPLIED OR STATUTORY, AND
+ * EXPRESSLY DISCLAIMS ANY IMPLIED WARRANTY OF MERCHANTABILITY, FITNESS
+ * FOR A PARTICULAR PURPOSE OR NONINFRINGEMENT OF THIRD PARTIES'
+ * RIGHTS. NO FLL DEALER, AGENT, EMPLOYEES IS AUTHORIZED TO MAKE ANY
+ * MODIFICATIONS, EXTENSIONS, OR ADDITIONS TO THIS WARRANTY.
+ * UNDER NO CIRCUMSTANCES AND UNDER NO LEGAL THEORY, TORT, CONTRACT,
+ * OR OTHERWISE, SHALL FLL BE LIABLE TO YOU OR ANY OTHER PERSON FOR ANY
+ * DIRECT, INDIRECT, SPECIAL, INCIDENTAL, PUNITIVE OR CONSEQUENTIAL
+ * DAMAGES OF ANY CHARACTER, INCLUDING, WITHOUT LIMITATION, DAMAGES
+ * ARISING OUT OF OR RELATING TO THE SOFTWARE OR THIS AGREEMENT, DAMAGES
+ * FOR LOSS OF GOODWILL, WORK STOPPAGE, OR LOSS OF DATA, OR FOR ANY
+ * DAMAGES, EVEN IF FLL SHALL HAVE BEEN INFORMED OF THE POSSIBILITY OF
+ * SUCH DAMAGES, OR FOR ANY CLAIM BY ANY OTHER PARTY. EVEN IF A PART
+ * OF THE SOFTWARE HAS BEEN DEVELOPED BY A THIRD PARTY, THE THIRD PARTY
+ * DEVELOPER SHALL HAVE NO LIABILITY IN CONNECTION WITH THE USE,
+ * PERFORMANCE OR NON-PERFORMANCE OF THE SOFTWARE.
+ *
+ * $OpenXM: OpenXM_contrib2/asir2000/builtin/array.c,v 1.24 2001/10/09 01:36:05 noro Exp $
+*/
#include "ca.h"
#include "base.h"
#include "parse.h"
@@ -9,23 +57,14 @@
#define DMAR(a1,a2,a3,d,r) (r)=dmar(a1,a2,a3,d);
#endif
-extern int Print; /* XXX */
+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 **);
-
-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();
@@ -41,20 +80,27 @@ void Px962_irredpoly_up2();
void Pirredpoly_up2();
void Pnbpoly_up2();
void Pqsort();
+void Pexponent_vector();
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},
+ {"vector",Pnewvect,-2},
+ {"exponent_vector",Pexponent_vector,-99999999},
{"newmat",Pnewmat,-3},
+ {"matrix",Pnewmat,-3},
+ {"newbytearray",Pnewbytearray,-2},
{"sepmat_destructive",Psepmat_destructive,2},
{"sepvect",Psepvect,2},
{"qsort",Pqsort,-2},
{"vtol",Pvtol,1},
{"size",Psize,1},
{"det",Pdet,-2},
+ {"invmat",Pinvmat,-2},
{"leqm",Pleqm,2},
{"leqm1",Pleqm1,2},
{"geninvm",Pgeninvm,2},
@@ -72,8 +118,7 @@ struct ftab array_tab[] = {
{0,0,0},
};
-int comp_obj(a,b)
-Obj *a,*b;
+int comp_obj(Obj *a,Obj *b)
{
return arf_comp(CO,*a,*b);
}
@@ -81,8 +126,7 @@ Obj *a,*b;
static FUNC generic_comp_obj_func;
static NODE generic_comp_obj_arg;
-int generic_comp_obj(a,b)
-Obj *a,*b;
+int generic_comp_obj(Obj *a,Obj *b)
{
Q r;
@@ -96,13 +140,9 @@ Obj *a,*b;
}
-void Pqsort(arg,rp)
-NODE arg;
-VECT *rp;
+void Pqsort(NODE arg,VECT *rp)
{
VECT vect;
- char buf[BUFSIZ];
- char *fname;
NODE n;
P p;
V v;
@@ -125,9 +165,7 @@ VECT *rp;
*rp = vect;
}
-void PNBmul_gf2n(arg,rp)
-NODE arg;
-GF2N *rp;
+void PNBmul_gf2n(NODE arg,GF2N *rp)
{
GF2N a,b;
GF2MAT mat;
@@ -164,9 +202,7 @@ GF2N *rp;
}
}
-void Pmul_vect_mat_gf2n(arg,rp)
-NODE arg;
-GF2N *rp;
+void Pmul_vect_mat_gf2n(NODE arg,GF2N *rp)
{
GF2N a;
GF2MAT mat;
@@ -197,9 +233,7 @@ GF2N *rp;
}
}
-void Pbconvmat_gf2n(arg,rp)
-NODE arg;
-LIST *rp;
+void Pbconvmat_gf2n(NODE arg,LIST *rp)
{
P p0,p1;
int to;
@@ -219,9 +253,7 @@ LIST *rp;
MKLIST(*rp,n0);
}
-void Pmulmat_gf2n(arg,rp)
-NODE arg;
-GF2MAT *rp;
+void Pmulmat_gf2n(NODE arg,GF2MAT *rp)
{
GF2MAT m;
@@ -230,9 +262,7 @@ GF2MAT *rp;
*rp = m;
}
-void Psepmat_destructive(arg,rp)
-NODE arg;
-LIST *rp;
+void Psepmat_destructive(NODE arg,LIST *rp)
{
MAT mat,mat1;
int i,j,row,col;
@@ -263,17 +293,12 @@ LIST *rp;
MKLIST(*rp,n0);
}
-void Psepvect(arg,rp)
-NODE arg;
-VECT *rp;
+void Psepvect(NODE arg,VECT *rp)
{
sepvect((VECT)ARG0(arg),QTOS((Q)ARG1(arg)),rp);
}
-void sepvect(v,d,rp)
-VECT v;
-int d;
-VECT *rp;
+void sepvect(VECT v,int d,VECT *rp)
{
int i,j,k,n,q,q1,r;
pointer *pv,*pw,*pu;
@@ -297,9 +322,7 @@ VECT *rp;
}
}
-void Pnewvect(arg,rp)
-NODE arg;
-VECT *rp;
+void Pnewvect(NODE arg,VECT *rp)
{
int len,i,r;
VECT vect;
@@ -326,10 +349,55 @@ VECT *rp;
*rp = vect;
}
-void Pnewmat(arg,rp)
-NODE arg;
-MAT *rp;
+void Pexponent_vector(NODE arg,DP *rp)
{
+ nodetod(arg,rp);
+}
+
+void Pnewbytearray(NODE arg,BYTEARRAY *rp)
+{
+ int len,i,r;
+ BYTEARRAY array;
+ unsigned char *vb;
+ char *str;
+ LIST list;
+ NODE tn;
+
+ asir_assert(ARG0(arg),O_N,"newbytearray");
+ len = QTOS((Q)ARG0(arg));
+ if ( len < 0 )
+ error("newbytearray : invalid size");
+ MKBYTEARRAY(array,len);
+ if ( argc(arg) == 2 ) {
+ if ( !ARG1(arg) )
+ error("newbytearray : invalid initialization");
+ switch ( OID((Obj)ARG1(arg)) ) {
+ case O_LIST:
+ list = (LIST)ARG1(arg);
+ asir_assert(list,O_LIST,"newbytearray");
+ for ( r = 0, tn = BDY(list); tn; r++, tn = NEXT(tn) );
+ if ( r <= len ) {
+ for ( i = 0, tn = BDY(list), vb = BDY(array); tn;
+ i++, tn = NEXT(tn) )
+ vb[i] = (unsigned char)QTOS((Q)BDY(tn));
+ }
+ break;
+ case O_STR:
+ str = BDY((STRING)ARG1(arg));
+ r = strlen(str);
+ if ( r <= len )
+ bcopy(str,BDY(array),r);
+ break;
+ default:
+ if ( !ARG1(arg) )
+ error("newbytearray : invalid initialization");
+ }
+ }
+ *rp = array;
+}
+
+void Pnewmat(NODE arg,MAT *rp)
+{
int row,col;
int i,j,r,c;
NODE tn,sn;
@@ -363,9 +431,7 @@ MAT *rp;
*rp = m;
}
-void Pvtol(arg,rp)
-NODE arg;
-LIST *rp;
+void Pvtol(NODE arg,LIST *rp)
{
NODE n,n1;
VECT v;
@@ -380,9 +446,7 @@ LIST *rp;
MKLIST(*rp,n);
}
-void Premainder(arg,rp)
-NODE arg;
-Obj *rp;
+void Premainder(NODE arg,Obj *rp)
{
Obj a;
VECT v,w;
@@ -430,9 +494,7 @@ Obj *rp;
}
}
-void Psremainder(arg,rp)
-NODE arg;
-Obj *rp;
+void Psremainder(NODE arg,Obj *rp)
{
Obj a;
VECT v,w;
@@ -481,9 +543,7 @@ Obj *rp;
}
}
-void Psize(arg,rp)
-NODE arg;
-LIST *rp;
+void Psize(NODE arg,LIST *rp)
{
int n,m;
@@ -509,9 +569,7 @@ LIST *rp;
MKLIST(*rp,t);
}
-void Pdet(arg,rp)
-NODE arg;
-P *rp;
+void Pdet(NODE arg,P *rp)
{
MAT m;
int n,i,j,mod;
@@ -535,6 +593,39 @@ 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+...
@@ -546,11 +637,48 @@ 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)
-NODE arg;
-LIST *rp;
+void Pgeneric_gauss_elim(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;
@@ -558,7 +686,7 @@ 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");
@@ -598,9 +726,7 @@ LIST *rp;
MKLIST(*rp,n0);
}
-void Pleqm(arg,rp)
-NODE arg;
-VECT *rp;
+void Pleqm(NODE arg,VECT *rp)
{
MAT m;
VECT vect;
@@ -639,9 +765,7 @@ VECT *rp;
}
}
-int gauss_elim_mod(mat,row,col,md)
-int **mat;
-int row,col,md;
+int gauss_elim_mod(int **mat,int row,int col,int md)
{
int i,j,k,inv,a,n;
int *t,*pivot;
@@ -664,8 +788,10 @@ int row,col,md;
t = mat[i];
if ( i != j && (a = t[j]) )
for ( k = j, a = md - a; k <= n; k++ ) {
+ unsigned int tk;
/* t[k] = dmar(pivot[k],a,t[k],md); */
- DMAR(pivot[k],a,t[k],md,t[k])
+ DMAR(pivot[k],a,t[k],md,tk)
+ t[k] = tk;
}
}
}
@@ -678,11 +804,7 @@ 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)
-MAT mat;
-MAT *nm;
-Q *dn;
-int **rindp,**cindp;
+int generic_gauss_elim(MAT mat,MAT *nm,Q *dn,int **rindp,int **cindp)
{
int **wmat;
Q **bmat;
@@ -709,10 +831,10 @@ int **rindp,**cindp;
colstat = (int *)MALLOC_ATOMIC(col*sizeof(int));
wcolstat = (int *)MALLOC_ATOMIC(col*sizeof(int));
for ( ind = 0; ; ind++ ) {
- if ( Print ) {
+ if ( DP_Print ) {
fprintf(asir_out,"."); fflush(asir_out);
}
- md = lprime[ind];
+ md = get_lprime(ind);
get_eg(&tmp0);
for ( i = 0; i < row; i++ )
for ( j = 0, bmi = bmat[i], wmi = wmat[i]; j < col; j++ )
@@ -746,13 +868,13 @@ RESET:
}
} else {
if ( rank < rank0 ) {
- if ( Print ) {
+ if ( DP_Print ) {
fprintf(asir_out,"lower rank matrix; continuing...\n");
fflush(asir_out);
}
continue;
} else if ( rank > rank0 ) {
- if ( Print ) {
+ if ( DP_Print ) {
fprintf(asir_out,"higher rank matrix; resetting...\n");
fflush(asir_out);
}
@@ -760,7 +882,7 @@ RESET:
} else {
for ( j = 0; (j
row; col = mat->col;
w = (int **)almat(row,col);
for ( ind = 0; ; ind++ ) {
- md = lprime[ind];
+ md = get_lprime(ind);
STOQ(md,mdq);
for ( i = 0; i < row; i++ )
for ( j = 0, ai = a0[i], wi = w[i]; j < col; j++ )
@@ -939,7 +1060,7 @@ int **rindp,**cindp;
add_eg(&eg_mul,&tmp0,&tmp1);
/* q = q*md */
mulq(q,mdq,&u); q = u;
- if ( !(count % 2) && intmtoratm_q(xmat,NM(q),*nmmat,dn) ) {
+ if ( !(count % 16) && intmtoratm_q(xmat,NM(q),*nmmat,dn) ) {
for ( j = k = l = 0; j < col; j++ )
if ( cinfo[j] )
rind[k++] = j;
@@ -959,13 +1080,10 @@ int **rindp,**cindp;
int f4_nocheck;
-int gensolve_check(mat,nm,dn,rind,cind)
-MAT mat,nm;
-Q dn;
-int *rind,*cind;
+int gensolve_check(MAT mat,MAT nm,Q dn,int *rind,int *cind)
{
int row,col,rank,clen,i,j,k,l;
- Q s,t,u;
+ Q s,t;
Q *w;
Q *mati,*nmk;
@@ -1010,13 +1128,10 @@ int *rind,*cind;
/* assuming 0 < c < m */
-int inttorat(c,m,b,sgnp,nmp,dnp)
-N c,m,b;
-int *sgnp;
-N *nmp,*dnp;
+int inttorat(N c,N m,N b,int *sgnp,N *nmp,N *dnp)
{
- Q qq,t,u1,v1,r1,nm;
- N q,r,u2,v2,r2;
+ Q qq,t,u1,v1,r1;
+ N q,u2,v2,r2;
u1 = 0; v1 = ONE; u2 = m; v2 = c;
while ( cmpn(v2,b) >= 0 ) {
@@ -1035,14 +1150,10 @@ N *nmp,*dnp;
/* mat->body = N ** */
-int intmtoratm(mat,md,nm,dn)
-MAT mat;
-N md;
-MAT nm;
-Q *dn;
+int intmtoratm(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;
@@ -1092,14 +1203,10 @@ Q *dn;
/* mat->body = Q ** */
-int intmtoratm_q(mat,md,nm,dn)
-MAT mat;
-N md;
-MAT nm;
-Q *dn;
+int intmtoratm_q(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;
@@ -1151,10 +1258,7 @@ Q *dn;
#define ONE_STEP1 if ( zzz = *s ) { DMAR(zzz,hc,*tj,md,*tj) } tj++; s++;
-void reduce_reducers_mod(mat,row,col,md)
-int **mat;
-int row,col;
-int md;
+void reduce_reducers_mod(int **mat,int row,int col,int md)
{
int i,j,k,l,hc,zzz;
int *t,*s,*tj,*ind;
@@ -1193,7 +1297,7 @@ int md;
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++;
}
}
@@ -1209,10 +1313,7 @@ int md;
2. reduce spolys by the reduced reducers
*/
-void pre_reduce_mod(mat,row,col,nred,md)
-int **mat;
-int row,col,nred;
-int md;
+void pre_reduce_mod(int **mat,int row,int col,int nred,int md)
{
int i,j,k,l,hc,inv;
int *t,*s,*tk,*ind;
@@ -1262,14 +1363,10 @@ 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)
-int *sp,**redmat;
-int *ind;
-int nred,col;
-int md;
+void reduce_sp_by_red_mod(int *sp,int **redmat,int *ind,int nred,int 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-- ) {
@@ -1280,25 +1377,97 @@ int md;
hc = md-hc;
s = redmat[i]+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++;
}
}
}
}
+/*
+ mat[i] : compressed reducers (i=0,...,nred-1)
+ mat[0] < mat[1] < ... < mat[nred-1] w.r.t the term order
+*/
+
+void red_by_compress(int m,unsigned int *p,unsigned int *r,
+ unsigned int *ri,unsigned int hc,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 */
+
+void red_by_vect(int m,unsigned int *p,unsigned int *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 (int *sp,CDP *redmat,int *ind,
+ int nred,int col,int md)
+{
+ int i,len;
+ 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] */
+ 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] >= (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 **mat;
-int row,col,md;
-int *colstat;
+int generic_gauss_elim_mod(int **mat0,int row,int col,int md,int *colstat)
{
- int i,j,k,l,inv,a,rank,zzz;
- int *t,*pivot,*pk,*tk;
+ int i,j,k,l,inv,a,rank;
+ 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;
@@ -1311,35 +1480,14 @@ 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] ) {
- a = md - a; pk = pivot+j; tk = t+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++;
- }
- }
+ if ( a = t[j] )
+ red_by_vect(md,t+j,pivot+j,md-a,col-j);
}
rank++;
}
@@ -1348,46 +1496,29 @@ int *colstat;
pivot = mat[l];
for ( i = 0; i < l; i++ ) {
t = mat[i];
- if ( a = t[j] ) {
- a = md-a; pk = pivot+j; tk = t+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++;
- }
- }
+ t[j] %= md;
+ if ( a = t[j] )
+ red_by_vect(md,t+j,pivot+j,md-a,col-j);
}
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)
-GFMMAT mat;
-unsigned int md;
-int *perm;
+int lu_gfmmat(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;
@@ -1414,7 +1545,10 @@ int *perm;
DMAR(inv,m,0,md,t[k])
for ( j = k+1, m = md - t[k]; j < col; j++ )
if ( pivot[j] ) {
- DMAR(m,pivot[j],t[j],md,t[j])
+ unsigned int tj;
+
+ DMAR(m,pivot[j],t[j],md,tj)
+ t[j] = tj;
}
}
}
@@ -1436,12 +1570,10 @@ 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)
-unsigned int **a;
-unsigned int md;
-int **rinfo,**cinfo;
+int find_lhs_and_lu_mod(unsigned int **a,int row,int col,
+ unsigned int md,int **rinfo,int **cinfo)
{
- int i,j,k,l,d;
+ int i,j,k,d;
int *rp,*cp;
unsigned int *t,*pivot;
unsigned int inv,m;
@@ -1469,7 +1601,9 @@ int **rinfo,**cinfo;
DMAR(inv,m,0,md,t[k])
for ( j = k+1, m = md - t[k]; j < col; j++ )
if ( pivot[j] ) {
- DMAR(m,pivot[j],t[j],md,t[j])
+ unsigned int tj;
+ DMAR(m,pivot[j],t[j],md,tj)
+ t[j] = tj;
}
}
}
@@ -1487,12 +1621,7 @@ int **rinfo,**cinfo;
b = a^(-1)b
*/
-void solve_by_lu_mod(a,n,md,b,l)
-int **a;
-int n;
-int md;
-int **b;
-int l;
+void solve_by_lu_mod(int **a,int n,int md,int **b,int l)
{
unsigned int *y,*c;
int i,j,k;
@@ -1530,9 +1659,7 @@ int l;
}
}
-void Pleqm1(arg,rp)
-NODE arg;
-VECT *rp;
+void Pleqm1(NODE arg,VECT *rp)
{
MAT m;
VECT vect;
@@ -1571,9 +1698,7 @@ VECT *rp;
}
}
-gauss_elim_mod1(mat,row,col,md)
-int **mat;
-int row,col,md;
+int gauss_elim_mod1(int **mat,int row,int col,int md)
{
int i,j,k,inv,a,n;
int *t,*pivot;
@@ -1610,9 +1735,7 @@ int row,col,md;
return -1;
}
-void Pgeninvm(arg,rp)
-NODE arg;
-LIST *rp;
+void Pgeninvm(NODE arg,LIST *rp)
{
MAT m;
pointer **mat;
@@ -1646,17 +1769,15 @@ 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)
-unsigned int **mat;
-int row,col,md;
+int gauss_elim_geninv_mod(unsigned int **mat,int row,int col,int md)
{
int i,j,k,inv,a,n,m;
unsigned int *t,*pivot;
@@ -1692,9 +1813,7 @@ int row,col,md;
return 0;
}
-void Psolve_by_lu_gfmmat(arg,rp)
-NODE arg;
-VECT *rp;
+void Psolve_by_lu_gfmmat(NODE arg,VECT *rp)
{
GFMMAT lu;
Q *perm,*rhs,*v;
@@ -1715,15 +1834,12 @@ 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)
-GFMMAT lu;
-unsigned int md;
-unsigned int *b;
-unsigned int *x;
+void solve_by_lu_gfmmat(GFMMAT lu,unsigned int md,
+ unsigned int *b,unsigned int *x)
{
int n;
unsigned int **a;
@@ -1755,9 +1871,7 @@ unsigned int *x;
}
}
-void Plu_gfmmat(arg,rp)
-NODE arg;
-LIST *rp;
+void Plu_gfmmat(NODE arg,LIST *rp)
{
MAT m;
GFMMAT mm;
@@ -1787,9 +1901,7 @@ LIST *rp;
MKLIST(*rp,n0);
}
-void Pmat_to_gfmmat(arg,rp)
-NODE arg;
-GFMMAT *rp;
+void Pmat_to_gfmmat(NODE arg,GFMMAT *rp)
{
MAT m;
unsigned int md;
@@ -1800,10 +1912,7 @@ GFMMAT *rp;
mat_to_gfmmat(m,md,rp);
}
-void mat_to_gfmmat(m,md,rp)
-MAT m;
-unsigned int md;
-GFMMAT *rp;
+void mat_to_gfmmat(MAT m,unsigned int md,GFMMAT *rp)
{
unsigned int **wmat;
unsigned int t;
@@ -1826,9 +1935,7 @@ GFMMAT *rp;
TOGFMMAT(row,col,wmat,*rp);
}
-void Pgeninvm_swap(arg,rp)
-NODE arg;
-LIST *rp;
+void Pgeninvm_swap(NODE arg,LIST *rp)
{
MAT m;
pointer **mat;
@@ -1874,12 +1981,8 @@ LIST *rp;
}
}
-gauss_elim_geninv_mod_swap(mat,row,col,md,invmatp,indexp)
-unsigned int **mat;
-int row,col;
-unsigned int md;
-unsigned int ***invmatp;
-int **indexp;
+int gauss_elim_geninv_mod_swap(unsigned int **mat,int row,int col,
+ unsigned int md,unsigned int ***invmatp,int **indexp)
{
int i,j,k,inv,a,n,m;
unsigned int *t,*pivot,*s;
@@ -1933,10 +2036,7 @@ void _addn(N,N,N);
int _subn(N,N,N);
void _muln(N,N,N);
-void inner_product_int(a,b,n,r)
-Q *a,*b;
-int n;
-Q *r;
+void inner_product_int(Q *a,Q *b,int n,Q *r)
{
int la,lb,i;
int sgn,sgn1;
@@ -1989,11 +2089,7 @@ 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)
-Q **a;
-int **b;
-int n,k,l;
-Q *r;
+void inner_product_mat_int_mod(Q **a,int **b,int n,int k,int l,Q *r)
{
int la,lb,i;
int sgn,sgn1;
@@ -2049,9 +2145,7 @@ Q *r;
NTOQ(sum,sgn,*r);
}
-void Pmul_mat_vect_int(arg,rp)
-NODE arg;
-VECT *rp;
+void Pmul_mat_vect_int(NODE arg,VECT *rp)
{
MAT mat;
VECT vect,r;
@@ -2062,14 +2156,13 @@ VECT *rp;
row = mat->row;
col = mat->col;
MKVECT(r,row);
- for ( i = 0; i < row; i++ )
- inner_product_int(mat->body[i],vect->body,col,&r->body[i]);
+ for ( i = 0; i < row; i++ ) {
+ inner_product_int((Q *)mat->body[i],(Q *)vect->body,col,(Q *)&r->body[i]);
+ }
*rp = r;
}
-void Pnbpoly_up2(arg,rp)
-NODE arg;
-GF2N *rp;
+void Pnbpoly_up2(NODE arg,GF2N *rp)
{
int m,type,ret;
UP2 r;
@@ -2083,11 +2176,9 @@ GF2N *rp;
*rp = 0;
}
-void Px962_irredpoly_up2(arg,rp)
-NODE arg;
-GF2N *rp;
+void Px962_irredpoly_up2(NODE arg,GF2N *rp)
{
- int m,type,ret,w;
+ int m,ret,w;
GF2N prev;
UP2 r;
@@ -2103,18 +2194,16 @@ 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)
-NODE arg;
-GF2N *rp;
+void Pirredpoly_up2(NODE arg,GF2N *rp)
{
- int m,type,ret,w;
+ int m,ret,w;
GF2N prev;
UP2 r;
@@ -2130,7 +2219,7 @@ 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
@@ -2410,23 +2499,19 @@ PENTA:
return 1;
}
-printqmat(mat,row,col)
-Q **mat;
-int row,col;
+void printqmat(Q **mat,int row,int col)
{
int i,j;
for ( i = 0; i < row; i++ ) {
for ( j = 0; j < col; j++ ) {
- printnum(mat[i][j]); printf(" ");
+ printnum((Num)mat[i][j]); printf(" ");
}
printf("\n");
}
}
-printimat(mat,row,col)
-int **mat;
-int row,col;
+void printimat(int **mat,int row,int col)
{
int i,j;