=================================================================== RCS file: /home/cvs/OpenXM_contrib2/asir2000/builtin/array.c,v retrieving revision 1.1 retrieving revision 1.24 diff -u -p -r1.1 -r1.24 --- OpenXM_contrib2/asir2000/builtin/array.c 1999/12/03 07:39:07 1.1 +++ OpenXM_contrib2/asir2000/builtin/array.c 2001/10/09 01:36:05 1.24 @@ -1,27 +1,68 @@ -/* $OpenXM: OpenXM/src/asir99/builtin/array.c,v 1.2 1999/11/23 07:14:14 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.23 2001/10/01 01:58:01 noro Exp $ +*/ #include "ca.h" #include "base.h" #include "parse.h" #include "inline.h" -/* + +#if 0 #undef DMAR #define DMAR(a1,a2,a3,d,r) (r)=dmar(a1,a2,a3,d); -*/ +#endif -extern int Print; /* XXX */ +extern int DP_Print; /* XXX */ -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_mod(); @@ -38,6 +79,7 @@ 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}, @@ -45,13 +87,18 @@ struct ftab array_tab[] = { {"mat_to_gfmmat",Pmat_to_gfmmat,2}, {"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}, @@ -69,8 +116,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); } @@ -78,8 +124,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; @@ -93,13 +138,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; @@ -122,9 +163,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; @@ -161,9 +200,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; @@ -194,9 +231,7 @@ GF2N *rp; } } -void Pbconvmat_gf2n(arg,rp) -NODE arg; -LIST *rp; +void Pbconvmat_gf2n(NODE arg,LIST *rp) { P p0,p1; int to; @@ -216,9 +251,7 @@ LIST *rp; MKLIST(*rp,n0); } -void Pmulmat_gf2n(arg,rp) -NODE arg; -GF2MAT *rp; +void Pmulmat_gf2n(NODE arg,GF2MAT *rp) { GF2MAT m; @@ -227,9 +260,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; @@ -260,17 +291,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; @@ -294,9 +320,7 @@ VECT *rp; } } -void Pnewvect(arg,rp) -NODE arg; -VECT *rp; +void Pnewvect(NODE arg,VECT *rp) { int len,i,r; VECT vect; @@ -306,7 +330,7 @@ VECT *rp; asir_assert(ARG0(arg),O_N,"newvect"); len = QTOS((Q)ARG0(arg)); - if ( len <= 0 ) + if ( len < 0 ) error("newvect : invalid size"); MKVECT(vect,len); if ( argc(arg) == 2 ) { @@ -323,10 +347,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; @@ -337,7 +406,7 @@ MAT *rp; asir_assert(ARG0(arg),O_N,"newmat"); asir_assert(ARG1(arg),O_N,"newmat"); row = QTOS((Q)ARG0(arg)); col = QTOS((Q)ARG1(arg)); - if ( row <= 0 || col <= 0 ) + if ( row < 0 || col < 0 ) error("newmat : invalid size"); MKMAT(m,row,col); if ( argc(arg) == 3 ) { @@ -360,9 +429,7 @@ MAT *rp; *rp = m; } -void Pvtol(arg,rp) -NODE arg; -LIST *rp; +void Pvtol(NODE arg,LIST *rp) { NODE n,n1; VECT v; @@ -377,9 +444,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; @@ -427,9 +492,7 @@ Obj *rp; } } -void Psremainder(arg,rp) -NODE arg; -Obj *rp; +void Psremainder(NODE arg,Obj *rp) { Obj a; VECT v,w; @@ -478,9 +541,7 @@ Obj *rp; } } -void Psize(arg,rp) -NODE arg; -LIST *rp; +void Psize(NODE arg,LIST *rp) { int n,m; @@ -506,9 +567,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; @@ -532,6 +591,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+... @@ -543,9 +635,7 @@ 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_mod(NODE arg,LIST *rp) { NODE n0; MAT m,mat; @@ -555,7 +645,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"); @@ -595,9 +685,7 @@ LIST *rp; MKLIST(*rp,n0); } -void Pleqm(arg,rp) -NODE arg; -VECT *rp; +void Pleqm(NODE arg,VECT *rp) { MAT m; VECT vect; @@ -636,9 +724,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; @@ -661,8 +747,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; } } } @@ -673,13 +761,9 @@ int row,col,md; return -1; } -struct oEGT eg_mod,eg_elim,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; -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; @@ -706,9 +790,10 @@ int **rindp,**cindp; colstat = (int *)MALLOC_ATOMIC(col*sizeof(int)); wcolstat = (int *)MALLOC_ATOMIC(col*sizeof(int)); for ( ind = 0; ; ind++ ) { - if ( Print ) - fprintf(asir_out,"."); - md = lprime[ind]; + if ( DP_Print ) { + fprintf(asir_out,"."); fflush(asir_out); + } + 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++ ) @@ -742,18 +827,24 @@ 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); + } goto RESET; } else { for ( j = 0; (jbody; + row = mat->row; col = mat->col; + w = (int **)almat(row,col); + for ( ind = 0; ; 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++ ) + if ( q = (Q)ai[j] ) { + t = rem(NM(q),md); + if ( t && SGN(q) < 0 ) + t = (md - t) % md; + wi[j] = t; + } else + wi[j] = 0; + + rank = find_lhs_and_lu_mod(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++ ) + if ( cinfo[j] ) { + /* the column is in lhs */ + for ( i = 0; i < rank; i++ ) { + w[i][li] = w[i][j]; + a[i][li] = a0[rinfo[i]][j]; + } + li++; + } else { + /* the column is in rhs */ + for ( i = 0; i < rank; i++ ) + b[i][ri] = a0[rinfo[i]][j]; + ri++; + } + + /* solve Ax+B=0; A: rank x rank, B: rank x ri */ + MKMAT(xmat,rank,ri); x = (Q **)(xmat)->body; + MKMAT(*nmmat,rank,ri); nm = (Q **)(*nmmat)->body; + /* use the right part of w as work area */ + /* ri = col - rank */ + wc = (int **)almat(rank,ri); + for ( i = 0; i < rank; i++ ) + wc[i] = w[i]+rank; + *rindp = rind = (int *)MALLOC_ATOMIC(rank*sizeof(int)); + *cindp = cind = (int *)MALLOC_ATOMIC((ri)*sizeof(int)); + + init_eg(&eg_mul); init_eg(&eg_inv); + for ( q = ONE, count = 0; ; count++ ) { + fprintf(stderr,"."); + /* wc = -b mod md */ + for ( i = 0; i < rank; i++ ) + for ( j = 0, bi = b[i], wi = wc[i]; j < ri; j++ ) + if ( u = (Q)bi[j] ) { + t = rem(NM(u),md); + if ( t && SGN(u) > 0 ) + t = (md - t) % md; + wi[j] = t; + } else + wi[j] = 0; + /* wc = A^(-1)wc; wc is normalized */ + get_eg(&tmp0); + solve_by_lu_mod(w,rank,md,wc,ri); + get_eg(&tmp1); + add_eg(&eg_inv,&tmp0,&tmp1); + /* x = x-q*wc */ + for ( i = 0; i < rank; i++ ) + for ( j = 0, xi = x[i], wi = wc[i]; j < ri; j++ ) { + STOQ(wi[j],u); mulq(q,u,&s); + subq(xi[j],s,&u); xi[j] = u; + } + get_eg(&tmp0); + for ( i = 0; i < rank; i++ ) + for ( j = 0; j < ri; j++ ) { + inner_product_mat_int_mod(a,wc,rank,i,j,&u); + addq(b[i][j],u,&s); + if ( s ) { + t = divin(NM(s),md,&tn); + if ( t ) + error("generic_gauss_elim_hensel:incosistent"); + NTOQ(tn,SGN(s),b[i][j]); + } else + b[i][j] = 0; + } + get_eg(&tmp1); + add_eg(&eg_mul,&tmp0,&tmp1); + /* q = q*md */ + mulq(q,mdq,&u); q = u; + if ( !(count % 16) && intmtoratm_q(xmat,NM(q),*nmmat,dn) ) { + for ( j = k = l = 0; j < col; j++ ) + if ( cinfo[j] ) + rind[k++] = j; + else + cind[l++] = j; + if ( gensolve_check(mat,*nmmat,*dn,rind,cind) ) { + fprintf(stderr,"\n"); + print_eg("INV",&eg_inv); + print_eg("MUL",&eg_mul); + fflush(asir_out); + return rank; + } + } + } + } +} + 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; @@ -873,13 +1087,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 ) { @@ -898,14 +1109,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; @@ -914,6 +1121,8 @@ Q *dn; N u,unm,udn; int sgn,ret; + if ( UNIN(md) ) + return 0; row = mat->row; col = mat->col; bshiftn(md,1,&t); isqrt(t,&s); @@ -951,16 +1160,273 @@ Q *dn; return 1; } -int generic_gauss_elim_mod(mat,row,col,md,colstat) -int **mat; -int row,col,md; -int *colstat; +/* mat->body = Q ** */ + +int intmtoratm_q(MAT mat,N md,MAT nm,Q *dn) { + N t,s,b; + Q dn0,dn1,nm1,q; + int i,j,k,l,row,col; + Q **rmat; + Q **tmat; + Q *tmi; + Q *nmk; + N u,unm,udn; + int sgn,ret; + + if ( UNIN(md) ) + return 0; + row = mat->row; col = mat->col; + bshiftn(md,1,&t); + isqrt(t,&s); + bshiftn(s,64,&b); + if ( !b ) + b = ONEN; + dn0 = ONE; + tmat = (Q **)mat->body; + rmat = (Q **)nm->body; + for ( i = 0; i < row; i++ ) + for ( j = 0, tmi = tmat[i]; j < col; j++ ) + if ( tmi[j] ) { + muln(NM(tmi[j]),NM(dn0),&s); + remn(s,md,&u); + ret = inttorat(u,md,b,&sgn,&unm,&udn); + if ( !ret ) + return 0; + else { + if ( SGN(tmi[j])<0 ) + sgn = -sgn; + NTOQ(unm,sgn,nm1); + NTOQ(udn,1,dn1); + if ( !UNIQ(dn1) ) { + for ( k = 0; k < i; k++ ) + for ( l = 0, nmk = rmat[k]; l < col; l++ ) { + mulq(nmk[l],dn1,&q); nmk[l] = q; + } + for ( l = 0, nmk = rmat[i]; l < j; l++ ) { + mulq(nmk[l],dn1,&q); nmk[l] = q; + } + } + rmat[i][j] = nm1; + mulq(dn0,dn1,&q); dn0 = q; + } + } + *dn = dn0; + return 1; +} + +#define ONE_STEP1 if ( zzz = *s ) { DMAR(zzz,hc,*tj,md,*tj) } tj++; s++; + +void reduce_reducers_mod(int **mat,int row,int col,int md) +{ + int i,j,k,l,hc,zzz; + int *t,*s,*tj,*ind; + + /* reduce the reducers */ + ind = (int *)ALLOCA(row*sizeof(int)); + for ( i = 0; i < row; i++ ) { + t = mat[i]; + for ( j = 0; j < col && !t[j]; j++ ); + /* register the position of the head term */ + ind[i] = j; + for ( l = i-1; l >= 0; l-- ) { + /* reduce mat[i] by mat[l] */ + if ( hc = t[ind[l]] ) { + /* mat[i] = mat[i]-hc*mat[l] */ + j = ind[l]; + s = mat[l]+j; + tj = t+j; + hc = md-hc; + k = col-j; + for ( ; k >= 64; k -= 64 ) { + 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 + 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 + 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 + 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-- ) { + if ( zzz = *s ) { DMAR(zzz,hc,*tj,md,*tj) } tj++; s++; + } + } + } + } +} + +/* + mat[i] : reducers (i=0,...,nred-1) + spolys (i=nred,...,row-1) + mat[0] < mat[1] < ... < mat[nred-1] w.r.t the term order + 1. reduce the reducers + 2. reduce spolys by the reduced reducers +*/ + +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; + +#if 1 + /* reduce the reducers */ + ind = (int *)ALLOCA(row*sizeof(int)); + for ( i = 0; i < nred; i++ ) { + /* make mat[i] monic and mat[i] by mat[0],...,mat[i-1] */ + t = mat[i]; + for ( j = 0; j < col && !t[j]; j++ ); + /* register the position of the head term */ + ind[i] = j; + inv = invm(t[j],md); + for ( k = j; k < col; k++ ) + if ( t[k] ) + DMAR(t[k],inv,0,md,t[k]) + for ( l = i-1; l >= 0; l-- ) { + /* reduce mat[i] by mat[l] */ + if ( hc = t[ind[l]] ) { + /* mat[i] = mat[i]-hc*mat[l] */ + for ( k = ind[l], hc = md-hc, s = mat[l]+k, tk = t+k; + k < col; k++, tk++, s++ ) + if ( *s ) + DMAR(*s,hc,*tk,md,*tk) + } + } + } + /* reduce the spolys */ + for ( i = nred; i < row; i++ ) { + t = mat[i]; + for ( l = nred-1; l >= 0; l-- ) { + /* reduce mat[i] by mat[l] */ + if ( hc = t[ind[l]] ) { + /* mat[i] = mat[i]-hc*mat[l] */ + for ( k = ind[l], hc = md-hc, s = mat[l]+k, tk = t+k; + k < col; k++, tk++, s++ ) + if ( *s ) + DMAR(*s,hc,*tk,md,*tk) + } + } + } +#endif +} +/* + mat[i] : reducers (i=0,...,nred-1) + mat[0] < mat[1] < ... < mat[nred-1] w.r.t the term order +*/ + +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 *s,*tj; + + /* reduce the spolys by redmat */ + for ( i = nred-1; i >= 0; i-- ) { + /* reduce sp by redmat[i] */ + if ( hc = sp[ind[i]] ) { + /* sp = sp-hc*redmat[i] */ + j = ind[i]; + hc = md-hc; + s = redmat[i]+j; + tj = sp+j; + 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(int **mat0,int row,int col,int md,int *colstat) +{ int i,j,k,l,inv,a,rank; - int *t,*pivot; + 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; @@ -971,17 +1437,16 @@ int *colstat; } pivot = mat[rank]; inv = invm(pivot[j],md); - for ( k = j; k < col; k++ ) - if ( pivot[k] ) { - DMAR(pivot[k],inv,0,md,pivot[k]) + 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] ) - for ( k = j, a = md - a; k < col; k++ ) - if ( pivot[k] ) { - DMAR(pivot[k],a,t[k],md,t[k]) - } + red_by_vect(md,t+j,pivot+j,md-a,col-j); } rank++; } @@ -990,26 +1455,29 @@ int *colstat; pivot = mat[l]; for ( i = 0; i < l; i++ ) { t = mat[i]; + t[j] %= md; if ( a = t[j] ) - for ( k = j, a = md-a; k < col; k++ ) - if ( pivot[k] ) { - DMAR(pivot[k],a,t[k],md,t[k]) - } + 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; @@ -1036,7 +1504,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; } } } @@ -1044,10 +1515,111 @@ int *perm; return 1; } -void Pleqm1(arg,rp) -NODE arg; -VECT *rp; +/* + Input + a: a row x col matrix + md : a modulus + + Output: + return : d = the rank of mat + a[0..(d-1)][0..(d-1)] : LU decomposition (a[i][i] = 1/U[i][i]) + rinfo: array of length row + cinfo: array of length col + i-th row in new a <-> rinfo[i]-th row in old a + cinfo[j]=1 <=> j-th column is contained in the LU decomp. +*/ + +int find_lhs_and_lu_mod(unsigned int **a,int row,int col, + unsigned int md,int **rinfo,int **cinfo) { + int i,j,k,d; + int *rp,*cp; + unsigned int *t,*pivot; + unsigned int 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 = invm(pivot[k],md); + for ( i = d+1; i < row; i++ ) { + t = a[i]; + if ( m = t[k] ) { + DMAR(inv,m,0,md,t[k]) + for ( j = k+1, m = md - t[k]; j < col; j++ ) + if ( pivot[j] ) { + unsigned int tj; + DMAR(m,pivot[j],t[j],md,tj) + t[j] = tj; + } + } + } + d++; + } + return d; +} + +/* + 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_mod(int **a,int n,int md,int **b,int l) +{ + unsigned int *y,*c; + int i,j,k; + unsigned int t,m,m2; + + y = (int *)MALLOC_ATOMIC(n*sizeof(int)); + c = (int *)MALLOC_ATOMIC(n*sizeof(int)); + m2 = md>>1; + for ( k = 0; k < l; k++ ) { + /* copy b[.][k] to c */ + for ( i = 0; i < n; i++ ) + c[i] = (unsigned int)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]; + DMAR(m,y[j],t,md,t) + } + 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]; + DMAR(m,c[j],t,md,t) + } + /* a[i][i] = 1/U[i][i] */ + DMAR(t,a[i][i],0,md,c[i]) + } + /* copy c to b[.][k] with normalization */ + for ( i = 0; i < n; i++ ) + b[i][k] = (int)(c[i]>m2 ? c[i]-md : c[i]); + } +} + +void Pleqm1(NODE arg,VECT *rp) +{ MAT m; VECT vect; pointer **mat; @@ -1085,9 +1657,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; @@ -1124,9 +1694,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; @@ -1160,17 +1728,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; @@ -1206,9 +1772,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; @@ -1229,15 +1793,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; @@ -1269,9 +1830,7 @@ unsigned int *x; } } -void Plu_gfmmat(arg,rp) -NODE arg; -LIST *rp; +void Plu_gfmmat(NODE arg,LIST *rp) { MAT m; GFMMAT mm; @@ -1301,9 +1860,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; @@ -1314,10 +1871,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; @@ -1340,9 +1894,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; @@ -1388,12 +1940,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; @@ -1447,10 +1995,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; @@ -1501,10 +2046,66 @@ Q *r; NTOQ(sum,sgn,*r); } -void Pmul_mat_vect_int(arg,rp) -NODE arg; -VECT *rp; +/* (k,l) element of a*b where a: .x n matrix, b: n x . integer matrix */ + +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; + N wm,wma,sum,t; + Q aki; + int bil,bilsgn; + struct oN tn; + + for ( la = 0, i = 0; i < n; i++ ) { + if ( aki = a[k][i] ) + if ( DN(aki) ) + error("inner_product_int : invalid argument"); + else + la = MAX(PL(NM(aki)),la); + } + lb = 1; + sgn = 0; + sum= NALLOC(la+lb+2); + bzero((char *)sum,(la+lb+3)*sizeof(unsigned int)); + wm = NALLOC(la+lb+2); + wma = NALLOC(la+lb+2); + for ( i = 0; i < n; i++ ) { + if ( !(aki = a[k][i]) || !(bil = b[i][l]) ) + continue; + tn.p = 1; + if ( bil > 0 ) { + tn.b[0] = bil; bilsgn = 1; + } else { + tn.b[0] = -bil; bilsgn = -1; + } + _muln(NM(aki),&tn,wm); + sgn1 = SGN(aki)*bilsgn; + if ( !sgn ) { + sgn = sgn1; + t = wm; wm = sum; sum = t; + } else if ( sgn == sgn1 ) { + _addn(sum,wm,wma); + if ( !PL(wma) ) + sgn = 0; + t = wma; wma = sum; sum = t; + } else { + /* sgn*sum+sgn1*wm = sgn*(sum-wm) */ + sgn *= _subn(sum,wm,wma); + t = wma; wma = sum; sum = t; + } + } + GC_free(wm); + GC_free(wma); + if ( !sgn ) { + GC_free(sum); + *r = 0; + } else + NTOQ(sum,sgn,*r); +} + +void Pmul_mat_vect_int(NODE arg,VECT *rp) +{ MAT mat; VECT vect,r; int row,col,i; @@ -1514,14 +2115,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; @@ -1535,11 +2135,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; @@ -1555,18 +2153,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; @@ -1582,7 +2178,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 @@ -1860,4 +2456,28 @@ PENTA: } /* exhausted */ return 1; +} + +void printqmat(Q **mat,int row,int col) +{ + int i,j; + + for ( i = 0; i < row; i++ ) { + for ( j = 0; j < col; j++ ) { + printnum((Num)mat[i][j]); printf(" "); + } + printf("\n"); + } +} + +void printimat(int **mat,int row,int col) +{ + int i,j; + + for ( i = 0; i < row; i++ ) { + for ( j = 0; j < col; j++ ) { + printf("%d ",mat[i][j]); + } + printf("\n"); + } }