/* $OpenXM: OpenXM_contrib2/asir2000/builtin/array.c,v 1.5 2000/06/05 02:26:47 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 */ 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 Pgeneric_gauss_elim_mod(); void Pmat_to_gfmmat(),Plu_gfmmat(),Psolve_by_lu_gfmmat(); void Pgeninvm_swap(), Premainder(), Psremainder(), Pvtol(); void sepvect(); void Pmulmat_gf2n(); void Pbconvmat_gf2n(); void Pmul_vect_mat_gf2n(); void PNBmul_gf2n(); void Pmul_mat_vect_int(); void Psepmat_destructive(); void Px962_irredpoly_up2(); void Pirredpoly_up2(); void Pnbpoly_up2(); void Pqsort(); 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_mod",Pgeneric_gauss_elim_mod,2}, {"newvect",Pnewvect,-2}, {"newmat",Pnewmat,-3}, {"sepmat_destructive",Psepmat_destructive,2}, {"sepvect",Psepvect,2}, {"qsort",Pqsort,-2}, {"vtol",Pvtol,1}, {"size",Psize,1}, {"det",Pdet,-2}, {"leqm",Pleqm,2}, {"leqm1",Pleqm1,2}, {"geninvm",Pgeninvm,2}, {"geninvm_swap",Pgeninvm_swap,2}, {"remainder",Premainder,2}, {"sremainder",Psremainder,2}, {"mulmat_gf2n",Pmulmat_gf2n,1}, {"bconvmat_gf2n",Pbconvmat_gf2n,-4}, {"mul_vect_mat_gf2n",Pmul_vect_mat_gf2n,2}, {"mul_mat_vect_int",Pmul_mat_vect_int,2}, {"nbmul_gf2n",PNBmul_gf2n,3}, {"x962_irredpoly_up2",Px962_irredpoly_up2,2}, {"irredpoly_up2",Pirredpoly_up2,2}, {"nbpoly_up2",Pnbpoly_up2,2}, {0,0,0}, }; int comp_obj(a,b) Obj *a,*b; { return arf_comp(CO,*a,*b); } static FUNC generic_comp_obj_func; static NODE generic_comp_obj_arg; int generic_comp_obj(a,b) Obj *a,*b; { Q r; BDY(generic_comp_obj_arg)=(pointer)(*a); BDY(NEXT(generic_comp_obj_arg))=(pointer)(*b); r = (Q)bevalf(generic_comp_obj_func,generic_comp_obj_arg); if ( !r ) return 0; else return SGN(r)>0?1:-1; } void Pqsort(arg,rp) NODE arg; VECT *rp; { VECT vect; char buf[BUFSIZ]; char *fname; NODE n; P p; V v; asir_assert(ARG0(arg),O_VECT,"qsort"); vect = (VECT)ARG0(arg); if ( argc(arg) == 1 ) qsort(BDY(vect),vect->len,sizeof(Obj),(int (*)(const void *,const void *))comp_obj); else { p = (P)ARG1(arg); if ( !p || OID(p)!=2 ) error("qsort : invalid argument"); v = VR(p); if ( (int)v->attr != V_SR ) error("qsort : no such function"); generic_comp_obj_func = (FUNC)v->priv; MKNODE(n,0,0); MKNODE(generic_comp_obj_arg,0,n); qsort(BDY(vect),vect->len,sizeof(Obj),(int (*)(const void *,const void *))generic_comp_obj); } *rp = vect; } void PNBmul_gf2n(arg,rp) NODE arg; GF2N *rp; { GF2N a,b; GF2MAT mat; int n,w; unsigned int *ab,*bb; UP2 r; a = (GF2N)ARG0(arg); b = (GF2N)ARG1(arg); mat = (GF2MAT)ARG2(arg); if ( !a || !b ) *rp = 0; else { n = mat->row; w = (n+BSH-1)/BSH; ab = (unsigned int *)ALLOCA(w*sizeof(unsigned int)); bzero((char *)ab,w*sizeof(unsigned int)); bcopy(a->body->b,ab,(a->body->w)*sizeof(unsigned int)); bb = (unsigned int *)ALLOCA(w*sizeof(unsigned int)); bzero((char *)bb,w*sizeof(unsigned int)); bcopy(b->body->b,bb,(b->body->w)*sizeof(unsigned int)); NEWUP2(r,w); bzero((char *)r->b,w*sizeof(unsigned int)); mul_nb(mat,ab,bb,r->b); r->w = w; _adjup2(r); if ( !r->w ) *rp = 0; else MKGF2N(r,*rp); } } void Pmul_vect_mat_gf2n(arg,rp) NODE arg; GF2N *rp; { GF2N a; GF2MAT mat; int n,w; unsigned int *b; UP2 r; a = (GF2N)ARG0(arg); mat = (GF2MAT)ARG1(arg); if ( !a ) *rp = 0; else { n = mat->row; w = (n+BSH-1)/BSH; b = (unsigned int *)ALLOCA(w*sizeof(unsigned int)); bzero((char *)b,w*sizeof(unsigned int)); bcopy(a->body->b,b,(a->body->w)*sizeof(unsigned int)); NEWUP2(r,w); bzero((char *)r->b,w*sizeof(unsigned int)); mulgf2vectmat(mat->row,b,mat->body,r->b); r->w = w; _adjup2(r); if ( !r->w ) *rp = 0; else { MKGF2N(r,*rp); } } } void Pbconvmat_gf2n(arg,rp) NODE arg; LIST *rp; { P p0,p1; int to; GF2MAT p01,p10; GF2N root; NODE n0,n1; p0 = (P)ARG0(arg); p1 = (P)ARG1(arg); to = ARG2(arg)?1:0; if ( argc(arg) == 4 ) { root = (GF2N)ARG3(arg); compute_change_of_basis_matrix_with_root(p0,p1,to,root,&p01,&p10); } else compute_change_of_basis_matrix(p0,p1,to,&p01,&p10); MKNODE(n1,p10,0); MKNODE(n0,p01,n1); MKLIST(*rp,n0); } void Pmulmat_gf2n(arg,rp) NODE arg; GF2MAT *rp; { GF2MAT m; if ( !compute_multiplication_matrix((P)ARG0(arg),&m) ) error("mulmat_gf2n : input is not a normal polynomial"); *rp = m; } void Psepmat_destructive(arg,rp) NODE arg; LIST *rp; { MAT mat,mat1; int i,j,row,col; Q **a,**a1; Q ent; N nm,mod,rem,quo; int sgn; NODE n0,n1; mat = (MAT)ARG0(arg); mod = NM((Q)ARG1(arg)); row = mat->row; col = mat->col; MKMAT(mat1,row,col); a = (Q **)mat->body; a1 = (Q **)mat1->body; for ( i = 0; i < row; i++ ) for ( j = 0; j < col; j++ ) { ent = a[i][j]; if ( !ent ) continue; nm = NM(ent); sgn = SGN(ent); divn(nm,mod,&quo,&rem); /* if ( quo != nm && rem != nm ) */ /* GC_free(nm); */ /* GC_free(ent); */ NTOQ(rem,sgn,a[i][j]); NTOQ(quo,sgn,a1[i][j]); } MKNODE(n1,mat1,0); MKNODE(n0,mat,n1); MKLIST(*rp,n0); } void Psepvect(arg,rp) NODE arg; VECT *rp; { sepvect((VECT)ARG0(arg),QTOS((Q)ARG1(arg)),rp); } void sepvect(v,d,rp) VECT v; int d; VECT *rp; { int i,j,k,n,q,q1,r; pointer *pv,*pw,*pu; VECT w,u; n = v->len; if ( d > n ) d = n; q = n/d; r = n%d; q1 = q+1; MKVECT(w,d); *rp = w; pv = BDY(v); pw = BDY(w); k = 0; for ( i = 0; i < r; i++ ) { MKVECT(u,q1); pw[i] = (pointer)u; for ( pu = BDY(u), j = 0; j < q1; j++, k++ ) pu[j] = pv[k]; } for ( ; i < d; i++ ) { MKVECT(u,q); pw[i] = (pointer)u; for ( pu = BDY(u), j = 0; j < q; j++, k++ ) pu[j] = pv[k]; } } void Pnewvect(arg,rp) NODE arg; VECT *rp; { int len,i,r; VECT vect; pointer *vb; LIST list; NODE tn; asir_assert(ARG0(arg),O_N,"newvect"); len = QTOS((Q)ARG0(arg)); if ( len < 0 ) error("newvect : invalid size"); MKVECT(vect,len); if ( argc(arg) == 2 ) { list = (LIST)ARG1(arg); asir_assert(list,O_LIST,"newvect"); for ( r = 0, tn = BDY(list); tn; r++, tn = NEXT(tn) ); if ( r > len ) { *rp = vect; return; } for ( i = 0, tn = BDY(list), vb = BDY(vect); tn; i++, tn = NEXT(tn) ) vb[i] = (pointer)BDY(tn); } *rp = vect; } void Pnewmat(arg,rp) NODE arg; MAT *rp; { int row,col; int i,j,r,c; NODE tn,sn; MAT m; pointer **mb; LIST list; 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 ) error("newmat : invalid size"); MKMAT(m,row,col); if ( argc(arg) == 3 ) { list = (LIST)ARG2(arg); asir_assert(list,O_LIST,"newmat"); for ( r = 0, c = 0, tn = BDY(list); tn; r++, tn = NEXT(tn) ) { for ( j = 0, sn = BDY((LIST)BDY(tn)); sn; j++, sn = NEXT(sn) ); c = MAX(c,j); } if ( (r > row) || (c > col) ) { *rp = m; return; } for ( i = 0, tn = BDY(list), mb = BDY(m); tn; i++, tn = NEXT(tn) ) { asir_assert(BDY(tn),O_LIST,"newmat"); for ( j = 0, sn = BDY((LIST)BDY(tn)); sn; j++, sn = NEXT(sn) ) mb[i][j] = (pointer)BDY(sn); } } *rp = m; } void Pvtol(arg,rp) NODE arg; LIST *rp; { NODE n,n1; VECT v; pointer *a; int len,i; asir_assert(ARG0(arg),O_VECT,"vtol"); v = (VECT)ARG0(arg); len = v->len; a = BDY(v); for ( i = len - 1, n = 0; i >= 0; i-- ) { MKNODE(n1,a[i],n); n = n1; } MKLIST(*rp,n); } void Premainder(arg,rp) NODE arg; Obj *rp; { Obj a; VECT v,w; MAT m,l; pointer *vb,*wb; pointer **mb,**lb; int id,i,j,n,row,col,t,smd,sgn; Q md,q; a = (Obj)ARG0(arg); md = (Q)ARG1(arg); if ( !a ) *rp = 0; else { id = OID(a); switch ( id ) { case O_N: case O_P: cmp(md,(P)a,(P *)rp); break; case O_VECT: smd = QTOS(md); v = (VECT)a; n = v->len; vb = v->body; MKVECT(w,n); wb = w->body; for ( i = 0; i < n; i++ ) { if ( q = (Q)vb[i] ) { sgn = SGN(q); t = rem(NM(q),smd); STOQ(t,q); if ( q ) SGN(q) = sgn; } wb[i] = (pointer)q; } *rp = (Obj)w; break; case O_MAT: m = (MAT)a; row = m->row; col = m->col; mb = m->body; MKMAT(l,row,col); lb = l->body; for ( i = 0; i < row; i++ ) for ( j = 0, vb = mb[i], wb = lb[i]; j < col; j++ ) cmp(md,(P)vb[j],(P *)&wb[j]); *rp = (Obj)l; break; default: error("remainder : invalid argument"); } } } void Psremainder(arg,rp) NODE arg; Obj *rp; { Obj a; VECT v,w; MAT m,l; pointer *vb,*wb; pointer **mb,**lb; unsigned int t,smd; int id,i,j,n,row,col; Q md,q; a = (Obj)ARG0(arg); md = (Q)ARG1(arg); if ( !a ) *rp = 0; else { id = OID(a); switch ( id ) { case O_N: case O_P: cmp(md,(P)a,(P *)rp); break; case O_VECT: smd = QTOS(md); v = (VECT)a; n = v->len; vb = v->body; MKVECT(w,n); wb = w->body; for ( i = 0; i < n; i++ ) { if ( q = (Q)vb[i] ) { t = (unsigned int)rem(NM(q),smd); if ( SGN(q) < 0 ) t = (smd - t) % smd; UTOQ(t,q); } wb[i] = (pointer)q; } *rp = (Obj)w; break; case O_MAT: m = (MAT)a; row = m->row; col = m->col; mb = m->body; MKMAT(l,row,col); lb = l->body; for ( i = 0; i < row; i++ ) for ( j = 0, vb = mb[i], wb = lb[i]; j < col; j++ ) cmp(md,(P)vb[j],(P *)&wb[j]); *rp = (Obj)l; break; default: error("remainder : invalid argument"); } } } void Psize(arg,rp) NODE arg; LIST *rp; { int n,m; Q q; NODE t,s; if ( !ARG0(arg) ) t = 0; else { switch (OID(ARG0(arg))) { case O_VECT: n = ((VECT)ARG0(arg))->len; STOQ(n,q); MKNODE(t,q,0); break; case O_MAT: n = ((MAT)ARG0(arg))->row; m = ((MAT)ARG0(arg))->col; STOQ(m,q); MKNODE(s,q,0); STOQ(n,q); MKNODE(t,q,s); break; default: error("size : invalid argument"); break; } } MKLIST(*rp,t); } void Pdet(arg,rp) NODE arg; P *rp; { MAT m; int n,i,j,mod; P d; P **mat,**w; m = (MAT)ARG0(arg); asir_assert(m,O_MAT,"det"); if ( m->row != m->col ) error("det : non-square matrix"); else if ( argc(arg) == 1 ) detp(CO,(P **)BDY(m),m->row,rp); 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]); detmp(CO,mod,w,n,&d); mptop(d,rp); } } /* 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(arg,rp) NODE arg; LIST *rp; { NODE n0; MAT m,mat; VECT rind,cind; Q **tmat; int **wmat; Q *rib,*cib; int *colstat; Q q; int md,i,j,k,l,row,col,t,n,rank; asir_assert(ARG0(arg),O_MAT,"generic_gauss_elim_mod"); asir_assert(ARG1(arg),O_N,"generic_gauss_elim_mod"); m = (MAT)ARG0(arg); md = QTOS((Q)ARG1(arg)); row = m->row; col = m->col; tmat = (Q **)m->body; wmat = (int **)almat(row,col); colstat = (int *)MALLOC_ATOMIC(col*sizeof(int)); for ( i = 0; i < row; i++ ) for ( j = 0; j < col; j++ ) if ( q = (Q)tmat[i][j] ) { t = rem(NM(q),md); if ( t && SGN(q) < 0 ) t = (md - t) % md; wmat[i][j] = t; } else wmat[i][j] = 0; rank = generic_gauss_elim_mod(wmat,row,col,md,colstat); MKMAT(mat,rank,col-rank); tmat = (Q **)mat->body; for ( i = 0; i < rank; i++ ) for ( j = k = 0; j < col; j++ ) if ( !colstat[j] ) { UTOQ(wmat[i][j],tmat[i][k]); k++; } MKVECT(rind,rank); MKVECT(cind,col-rank); rib = (Q *)rind->body; cib = (Q *)cind->body; for ( j = k = l = 0; j < col; j++ ) if ( colstat[j] ) { STOQ(j,rib[k]); k++; } else { STOQ(j,cib[l]); l++; } n0 = mknode(3,mat,rind,cind); MKLIST(*rp,n0); } void Pleqm(arg,rp) NODE arg; VECT *rp; { MAT m; VECT vect; pointer **mat; Q *v; Q q; int **wmat; int md,i,j,row,col,t,n,status; asir_assert(ARG0(arg),O_MAT,"leqm"); asir_assert(ARG1(arg),O_N,"leqm"); m = (MAT)ARG0(arg); md = QTOS((Q)ARG1(arg)); row = m->row; col = m->col; mat = m->body; wmat = (int **)almat(row,col); for ( i = 0; i < row; i++ ) for ( j = 0; j < col; j++ ) if ( q = (Q)mat[i][j] ) { t = rem(NM(q),md); if ( SGN(q) < 0 ) t = (md - t) % md; wmat[i][j] = t; } else wmat[i][j] = 0; status = gauss_elim_mod(wmat,row,col,md); if ( status < 0 ) *rp = 0; else if ( status > 0 ) *rp = (VECT)ONE; else { n = col - 1; MKVECT(vect,n); for ( i = 0, v = (Q *)vect->body; i < n; i++ ) { t = (md-wmat[i][n])%md; STOQ(t,v[i]); } *rp = vect; } } int gauss_elim_mod(mat,row,col,md) int **mat; int row,col,md; { int i,j,k,inv,a,n; int *t,*pivot; n = col - 1; for ( j = 0; j < n; j++ ) { for ( i = j; i < row && !mat[i][j]; i++ ); if ( i == row ) return 1; if ( i != j ) { t = mat[i]; mat[i] = mat[j]; mat[j] = t; } pivot = mat[j]; inv = invm(pivot[j],md); for ( k = j; k <= n; k++ ) { /* pivot[k] = dmar(pivot[k],inv,0,md); */ DMAR(pivot[k],inv,0,md,pivot[k]) } for ( i = 0; i < row; i++ ) { t = mat[i]; if ( i != j && (a = t[j]) ) for ( k = j, a = md - a; k <= n; k++ ) { /* t[k] = dmar(pivot[k],a,t[k],md); */ DMAR(pivot[k],a,t[k],md,t[k]) } } } for ( i = n; i < row && !mat[i][n]; i++ ); if ( i == row ) return 0; else return -1; } 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 **wmat; Q **bmat; N **tmat; Q *bmi; N *tmi; Q q; int *wmi; int *colstat,*wcolstat,*rind,*cind; int row,col,ind,md,i,j,k,l,t,t1,rank,rank0,inv; N m1,m2,m3,s,u; MAT r,crmat; struct oEGT tmp0,tmp1; struct oEGT eg_mod_split,eg_elim_split,eg_chrem_split; struct oEGT eg_intrat_split,eg_gschk_split; int ret; init_eg(&eg_mod_split); init_eg(&eg_chrem_split); init_eg(&eg_elim_split); init_eg(&eg_intrat_split); init_eg(&eg_gschk_split); bmat = (Q **)mat->body; row = mat->row; col = mat->col; wmat = (int **)almat(row,col); colstat = (int *)MALLOC_ATOMIC(col*sizeof(int)); wcolstat = (int *)MALLOC_ATOMIC(col*sizeof(int)); for ( ind = 0; ; ind++ ) { if ( Print ) { fprintf(asir_out,"."); fflush(asir_out); } md = lprime[ind]; get_eg(&tmp0); for ( i = 0; i < row; i++ ) for ( j = 0, bmi = bmat[i], wmi = wmat[i]; j < col; j++ ) if ( q = (Q)bmi[j] ) { t = rem(NM(q),md); if ( t && SGN(q) < 0 ) t = (md - t) % md; wmi[j] = t; } else wmi[j] = 0; get_eg(&tmp1); add_eg(&eg_mod,&tmp0,&tmp1); add_eg(&eg_mod_split,&tmp0,&tmp1); get_eg(&tmp0); rank = generic_gauss_elim_mod(wmat,row,col,md,wcolstat); get_eg(&tmp1); add_eg(&eg_elim,&tmp0,&tmp1); add_eg(&eg_elim_split,&tmp0,&tmp1); if ( !ind ) { RESET: UTON(md,m1); rank0 = rank; bcopy(wcolstat,colstat,col*sizeof(int)); MKMAT(crmat,rank,col-rank); MKMAT(r,rank,col-rank); *nm = r; tmat = (N **)crmat->body; for ( i = 0; i < rank; i++ ) for ( j = k = 0, tmi = tmat[i], wmi = wmat[i]; j < col; j++ ) if ( !colstat[j] ) { UTON(wmi[j],tmi[k]); k++; } } else { if ( rank < rank0 ) { if ( Print ) { fprintf(asir_out,"lower rank matrix; continuing...\n"); fflush(asir_out); } continue; } else if ( rank > rank0 ) { if ( Print ) { fprintf(asir_out,"higher rank matrix; resetting...\n"); fflush(asir_out); } goto RESET; } else { for ( j = 0; (j= t ) t = wmi[j]-t; else t = md-(t-wmi[j]); DMAR(t,inv,0,md,t1) UTON(t1,u); muln(m1,u,&s); addn(tmi[k],s,&u); tmi[k] = u; } else if ( wmi[j] ) { /* f3 = m1*(m1 mod m2)^(-1)*f2 */ DMAR(wmi[j],inv,0,md,t) UTON(t,u); muln(m1,u,&s); tmi[k] = s; } k++; } m1 = m3; get_eg(&tmp1); add_eg(&eg_chrem,&tmp0,&tmp1); add_eg(&eg_chrem_split,&tmp0,&tmp1); get_eg(&tmp0); ret = intmtoratm(crmat,m1,*nm,dn); get_eg(&tmp1); add_eg(&eg_intrat,&tmp0,&tmp1); add_eg(&eg_intrat_split,&tmp0,&tmp1); if ( ret ) { *rindp = rind = (int *)MALLOC_ATOMIC(rank*sizeof(int)); *cindp = cind = (int *)MALLOC_ATOMIC((col-rank)*sizeof(int)); for ( j = k = l = 0; j < col; j++ ) if ( colstat[j] ) rind[k++] = j; else cind[l++] = j; get_eg(&tmp0); if ( gensolve_check(mat,*nm,*dn,rind,cind) ) { get_eg(&tmp1); add_eg(&eg_gschk,&tmp0,&tmp1); add_eg(&eg_gschk_split,&tmp0,&tmp1); if ( Print ) { print_eg("Mod",&eg_mod_split); print_eg("Elim",&eg_elim_split); print_eg("ChRem",&eg_chrem_split); print_eg("IntRat",&eg_intrat_split); print_eg("Check",&eg_gschk_split); fflush(asir_out); } return rank; } } } } } int generic_gauss_elim_hensel(mat,nmmat,dn,rindp,cindp) MAT mat; MAT *nmmat; Q *dn; int **rindp,**cindp; { MAT bmat,xmat; Q **a0,**a,**b,**x,**nm; Q *ai,*bi,*xi; int row,col; int **w; int *wi; int **wc; Q mdq,q,s,u; N tn; int ind,md,i,j,k,l,li,ri,rank; unsigned int t; int *cinfo,*rinfo; int *rind,*cind; int count; struct oEGT eg_mul,eg_inv,tmp0,tmp1; a0 = (Q **)mat->body; row = mat->row; col = mat->col; w = (int **)almat(row,col); for ( ind = 0; ; ind++ ) { md = 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 % 2) && 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 row,col,rank,clen,i,j,k,l; Q s,t,u; Q *w; Q *mati,*nmk; if ( f4_nocheck ) return 1; row = mat->row; col = mat->col; rank = nm->row; clen = nm->col; w = (Q *)MALLOC(clen*sizeof(Q)); for ( i = 0; i < row; i++ ) { mati = (Q *)mat->body[i]; #if 1 bzero(w,clen*sizeof(Q)); for ( k = 0; k < rank; k++ ) for ( l = 0, nmk = (Q *)nm->body[k]; l < clen; l++ ) { mulq(mati[rind[k]],nmk[l],&t); addq(w[l],t,&s); w[l] = s; } for ( j = 0; j < clen; j++ ) { mulq(dn,mati[cind[j]],&t); if ( cmpq(w[j],t) ) break; } #else for ( j = 0; j < clen; j++ ) { for ( k = 0, s = 0; k < rank; k++ ) { mulq(mati[rind[k]],nm->body[k][j],&t); addq(s,t,&u); s = u; } mulq(dn,mati[cind[j]],&t); if ( cmpq(s,t) ) break; } #endif if ( j != clen ) break; } if ( i != row ) return 0; else return 1; } /* assuming 0 < c < m */ int inttorat(c,m,b,sgnp,nmp,dnp) N c,m,b; int *sgnp; N *nmp,*dnp; { Q qq,t,u1,v1,r1,nm; N q,r,u2,v2,r2; u1 = 0; v1 = ONE; u2 = m; v2 = c; while ( cmpn(v2,b) >= 0 ) { divn(u2,v2,&q,&r2); u2 = v2; v2 = r2; NTOQ(q,1,qq); mulq(qq,v1,&t); subq(u1,t,&r1); u1 = v1; v1 = r1; } if ( cmpn(NM(v1),b) >= 0 ) return 0; else { *nmp = v2; *dnp = NM(v1); *sgnp = SGN(v1); return 1; } } /* mat->body = N ** */ int intmtoratm(mat,md,nm,dn) MAT mat; N md; MAT nm; Q *dn; { N t,s,b; Q bound,dn0,dn1,nm1,q,tq; int i,j,k,l,row,col; Q **rmat; N **tmat; N *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 = (N **)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(tmi[j],NM(dn0),&s); remn(s,md,&u); ret = inttorat(u,md,b,&sgn,&unm,&udn); if ( !ret ) return 0; else { 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; } /* mat->body = Q ** */ int intmtoratm_q(mat,md,nm,dn) MAT mat; N md; MAT nm; Q *dn; { N t,s,b; Q bound,dn0,dn1,nm1,q,tq; 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(mat,row,col,md) int **mat; int row,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(mat,row,col,nred,md) int **mat; int row,col,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(sp,redmat,ind,nred,col,md) int *sp,**redmat; int *ind; int nred,col; int md; { int i,j,k,hc,zzz; int *t,*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++; } } } } #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 i,j,k,l,inv,a,rank,zzz; int *t,*pivot,*pk,*tk; for ( rank = 0, j = 0; j < col; j++ ) { for ( i = rank; i < row && !mat[i][j]; i++ ); if ( i == row ) { colstat[j] = 0; continue; } else colstat[j] = 1; if ( i != rank ) { t = mat[i]; mat[i] = mat[rank]; mat[rank] = t; } pivot = mat[rank]; inv = invm(pivot[j],md); for ( k = j, pk = pivot+k; k < col; k++, pk++ ) if ( *pk ) { 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++; } } } rank++; } for ( j = col-1, l = rank-1; j >= 0; j-- ) if ( colstat[j] ) { pivot = mat[l]; for ( i = 0; i < l; i++ ) { t = mat[i]; if ( a = t[j] ) { 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++; } } } 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 row,col; int i,j,k,l; unsigned int *t,*pivot; unsigned int **a; unsigned int inv,m; row = mat->row; col = mat->col; a = mat->body; bzero(perm,row*sizeof(int)); for ( i = 0; i < row; i++ ) perm[i] = i; for ( k = 0; k < col; k++ ) { for ( i = k; i < row && !a[i][k]; i++ ); if ( i == row ) return 0; if ( i != k ) { j = perm[i]; perm[i] = perm[k]; perm[k] = j; t = a[i]; a[i] = a[k]; a[k] = t; } pivot = a[k]; pivot[k] = inv = invm(pivot[k],md); for ( i = k+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] ) { DMAR(m,pivot[j],t[j],md,t[j]) } } } } return 1; } /* 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(a,row,col,md,rinfo,cinfo) unsigned int **a; unsigned int md; int **rinfo,**cinfo; { int i,j,k,l,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] ) { DMAR(m,pivot[j],t[j],md,t[j]) } } } 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(a,n,md,b,l) 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(arg,rp) NODE arg; VECT *rp; { MAT m; VECT vect; pointer **mat; Q *v; Q q; int **wmat; int md,i,j,row,col,t,n,status; asir_assert(ARG0(arg),O_MAT,"leqm1"); asir_assert(ARG1(arg),O_N,"leqm1"); m = (MAT)ARG0(arg); md = QTOS((Q)ARG1(arg)); row = m->row; col = m->col; mat = m->body; wmat = (int **)almat(row,col); for ( i = 0; i < row; i++ ) for ( j = 0; j < col; j++ ) if ( q = (Q)mat[i][j] ) { t = rem(NM(q),md); if ( SGN(q) < 0 ) t = (md - t) % md; wmat[i][j] = t; } else wmat[i][j] = 0; status = gauss_elim_mod1(wmat,row,col,md); if ( status < 0 ) *rp = 0; else if ( status > 0 ) *rp = (VECT)ONE; else { n = col - 1; MKVECT(vect,n); for ( i = 0, v = (Q *)vect->body; i < n; i++ ) { t = (md-wmat[i][n])%md; STOQ(t,v[i]); } *rp = vect; } } gauss_elim_mod1(mat,row,col,md) int **mat; int row,col,md; { int i,j,k,inv,a,n; int *t,*pivot; n = col - 1; for ( j = 0; j < n; j++ ) { for ( i = j; i < row && !mat[i][j]; i++ ); if ( i == row ) return 1; if ( i != j ) { t = mat[i]; mat[i] = mat[j]; mat[j] = t; } pivot = mat[j]; inv = invm(pivot[j],md); for ( k = j; k <= n; k++ ) pivot[k] = dmar(pivot[k],inv,0,md); for ( i = j+1; i < row; i++ ) { t = mat[i]; if ( i != j && (a = t[j]) ) for ( k = j, a = md - a; k <= n; k++ ) t[k] = dmar(pivot[k],a,t[k],md); } } for ( i = n; i < row && !mat[i][n]; i++ ); if ( i == row ) { for ( j = n-1; j >= 0; j-- ) { for ( i = j-1, a = (md-mat[j][n])%md; i >= 0; i-- ) { mat[i][n] = dmar(mat[i][j],a,mat[i][n],md); mat[i][j] = 0; } } return 0; } else return -1; } void Pgeninvm(arg,rp) NODE arg; LIST *rp; { MAT m; pointer **mat; Q **tmat; Q q; unsigned int **wmat; int md,i,j,row,col,t,status; MAT mat1,mat2; NODE node1,node2; asir_assert(ARG0(arg),O_MAT,"leqm1"); asir_assert(ARG1(arg),O_N,"leqm1"); m = (MAT)ARG0(arg); md = QTOS((Q)ARG1(arg)); row = m->row; col = m->col; mat = m->body; wmat = (unsigned 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 = (Q)mat[i][j] ) { t = rem(NM(q),md); if ( SGN(q) < 0 ) t = (md - t) % md; wmat[i][j] = t; } wmat[i][col+i] = 1; } status = gauss_elim_geninv_mod(wmat,row,col,md); if ( status > 0 ) *rp = 0; else { 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]); for ( tmat = (Q **)mat2->body; i < row; i++ ) for ( j = 0; j < row; j++ ) STOQ(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 i,j,k,inv,a,n,m; unsigned int *t,*pivot; n = col; m = row+col; for ( j = 0; j < n; j++ ) { for ( i = j; i < row && !mat[i][j]; i++ ); if ( i == row ) return 1; if ( i != j ) { t = mat[i]; mat[i] = mat[j]; mat[j] = t; } pivot = mat[j]; inv = invm(pivot[j],md); for ( k = j; k < m; k++ ) pivot[k] = dmar(pivot[k],inv,0,md); for ( i = j+1; i < row; i++ ) { t = mat[i]; if ( a = t[j] ) for ( k = j, a = md - a; k < m; k++ ) t[k] = dmar(pivot[k],a,t[k],md); } } 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 = md - a; k < m; k++ ) t[k] = dmar(pivot[k],a,t[k],md); } } return 0; } void Psolve_by_lu_gfmmat(arg,rp) NODE arg; VECT *rp; { GFMMAT lu; Q *perm,*rhs,*v; int n,i; unsigned int md; unsigned int *b,*sol; VECT r; lu = (GFMMAT)ARG0(arg); perm = (Q *)BDY((VECT)ARG1(arg)); rhs = (Q *)BDY((VECT)ARG2(arg)); md = (unsigned int)QTOS((Q)ARG3(arg)); n = lu->col; b = (unsigned int *)MALLOC_ATOMIC(n*sizeof(int)); sol = (unsigned int *)MALLOC_ATOMIC(n*sizeof(int)); for ( i = 0; i < n; i++ ) b[i] = QTOS(rhs[QTOS(perm[i])]); 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]); *rp = r; } void solve_by_lu_gfmmat(lu,md,b,x) GFMMAT lu; unsigned int md; unsigned int *b; unsigned int *x; { int n; unsigned int **a; unsigned int *y; int i,j; unsigned int t,m; n = lu->col; a = lu->body; y = (unsigned int *)MALLOC_ATOMIC(n*sizeof(int)); /* solve Ly=b */ for ( i = 0; i < n; i++ ) { for ( t = b[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 Ux=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,x[j],t,md,t) } /* a[i][i] = 1/U[i][i] */ DMAR(t,a[i][i],0,md,x[i]) } } void Plu_gfmmat(arg,rp) NODE arg; LIST *rp; { MAT m; GFMMAT mm; unsigned int md; int i,row,col,status; int *iperm; Q *v; VECT perm; NODE n0; asir_assert(ARG0(arg),O_MAT,"mat_to_gfmmat"); asir_assert(ARG1(arg),O_N,"mat_to_gfmmat"); m = (MAT)ARG0(arg); md = (unsigned int)QTOS((Q)ARG1(arg)); mat_to_gfmmat(m,md,&mm); row = m->row; col = m->col; iperm = (int *)MALLOC_ATOMIC(row*sizeof(int)); status = lu_gfmmat(mm,md,iperm); if ( !status ) n0 = 0; else { MKVECT(perm,row); for ( i = 0, v = (Q *)perm->body; i < row; i++ ) STOQ(iperm[i],v[i]); n0 = mknode(2,mm,perm); } MKLIST(*rp,n0); } void Pmat_to_gfmmat(arg,rp) NODE arg; GFMMAT *rp; { MAT m; unsigned int md; asir_assert(ARG0(arg),O_MAT,"mat_to_gfmmat"); asir_assert(ARG1(arg),O_N,"mat_to_gfmmat"); m = (MAT)ARG0(arg); md = (unsigned int)QTOS((Q)ARG1(arg)); mat_to_gfmmat(m,md,rp); } void mat_to_gfmmat(m,md,rp) MAT m; unsigned int md; GFMMAT *rp; { unsigned int **wmat; unsigned int t; Q **mat; Q q; int i,j,row,col; row = m->row; col = m->col; mat = (Q **)m->body; wmat = (unsigned int **)almat(row,col); for ( i = 0; i < row; i++ ) { bzero((char *)wmat[i],col*sizeof(unsigned int)); for ( j = 0; j < col; j++ ) if ( q = mat[i][j] ) { t = (unsigned int)rem(NM(q),md); if ( SGN(q) < 0 ) t = (md - t) % md; wmat[i][j] = t; } } TOGFMMAT(row,col,wmat,*rp); } void Pgeninvm_swap(arg,rp) NODE arg; LIST *rp; { MAT m; pointer **mat; Q **tmat; Q *tvect; Q q; unsigned int **wmat,**invmat; int *index; unsigned int t,md; int i,j,row,col,status; MAT mat1; VECT vect1; NODE node1,node2; asir_assert(ARG0(arg),O_MAT,"geninvm_swap"); asir_assert(ARG1(arg),O_N,"geninvm_swap"); m = (MAT)ARG0(arg); md = QTOS((Q)ARG1(arg)); row = m->row; col = m->col; mat = m->body; wmat = (unsigned 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 = (Q)mat[i][j] ) { t = (unsigned int)rem(NM(q),md); if ( SGN(q) < 0 ) t = (md - t) % md; wmat[i][j] = t; } wmat[i][col+i] = 1; } status = gauss_elim_geninv_mod_swap(wmat,row,col,md,&invmat,&index); if ( status > 0 ) *rp = 0; else { MKMAT(mat1,col,col); for ( i = 0, tmat = (Q **)mat1->body; i < col; i++ ) for ( j = 0; j < col; j++ ) UTOQ(invmat[i][j],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); } } 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 i,j,k,inv,a,n,m; unsigned int *t,*pivot,*s; int *index; unsigned 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 = (unsigned int)invm(pivot[j],md); for ( k = j; k < m; k++ ) if ( pivot[k] ) pivot[k] = (unsigned int)dmar(pivot[k],inv,0,md); for ( i = j+1; i < row; i++ ) { t = mat[i]; if ( a = t[j] ) for ( k = j, a = md - a; k < m; k++ ) if ( pivot[k] ) t[k] = dmar(pivot[k],a,t[k],md); } } 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 = md - a; k < m; k++ ) if ( pivot[k] ) t[k] = dmar(pivot[k],a,t[k],md); } } *invmatp = invmat = (unsigned 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) Q *a,*b; int n; Q *r; { int la,lb,i; int sgn,sgn1; N wm,wma,sum,t; for ( la = lb = 0, i = 0; i < n; i++ ) { if ( a[i] ) if ( DN(a[i]) ) error("inner_product_int : invalid argument"); else la = MAX(PL(NM(a[i])),la); if ( b[i] ) if ( DN(b[i]) ) error("inner_product_int : invalid argument"); else lb = MAX(PL(NM(b[i])),lb); } 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 ( !a[i] || !b[i] ) continue; _muln(NM(a[i]),NM(b[i]),wm); sgn1 = SGN(a[i])*SGN(b[i]); 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); } /* (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; { 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(arg,rp) NODE arg; VECT *rp; { MAT mat; VECT vect,r; int row,col,i; mat = (MAT)ARG0(arg); vect = (VECT)ARG1(arg); 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]); *rp = r; } void Pnbpoly_up2(arg,rp) NODE arg; GF2N *rp; { int m,type,ret; UP2 r; m = QTOS((Q)ARG0(arg)); type = QTOS((Q)ARG1(arg)); ret = generate_ONB_polynomial(&r,m,type); if ( ret == 0 ) MKGF2N(r,*rp); else *rp = 0; } void Px962_irredpoly_up2(arg,rp) NODE arg; GF2N *rp; { int m,type,ret,w; GF2N prev; UP2 r; m = QTOS((Q)ARG0(arg)); prev = (GF2N)ARG1(arg); if ( !prev ) { w = (m>>5)+1; NEWUP2(r,w); r->w = 0; bzero((char *)r->b,w*sizeof(unsigned int)); } else { r = prev->body; if ( degup2(r) != m ) { w = (m>>5)+1; NEWUP2(r,w); r->w = 0; bzero((char *)r->b,w*sizeof(unsigned int)); } } ret = _generate_irreducible_polynomial(r,m,type); if ( ret == 0 ) MKGF2N(r,*rp); else *rp = 0; } void Pirredpoly_up2(arg,rp) NODE arg; GF2N *rp; { int m,type,ret,w; GF2N prev; UP2 r; m = QTOS((Q)ARG0(arg)); prev = (GF2N)ARG1(arg); if ( !prev ) { w = (m>>5)+1; NEWUP2(r,w); r->w = 0; bzero((char *)r->b,w*sizeof(unsigned int)); } else { r = prev->body; if ( degup2(r) != m ) { w = (m>>5)+1; NEWUP2(r,w); r->w = 0; bzero((char *)r->b,w*sizeof(unsigned int)); } } ret = _generate_good_irreducible_polynomial(r,m,type); if ( ret == 0 ) MKGF2N(r,*rp); else *rp = 0; } /* * f = type 'type' normal polynomial of degree m if exists * IEEE P1363 A.7.2 * * return value : 0 --- exists * 1 --- does not exist * -1 --- failure (memory allocation error) */ int generate_ONB_polynomial(UP2 *rp,int m,int type) { int i,r; int w; UP2 f,f0,f1,f2,t; w = (m>>5)+1; switch ( type ) { case 1: if ( !TypeT_NB_check(m,1) ) return 1; NEWUP2(f,w); *rp = f; f->w = w; /* set all the bits */ for ( i = 0; i < w; i++ ) f->b[i] = 0xffffffff; /* mask the top word if necessary */ if ( r = (m+1)&31 ) f->b[w-1] &= (1<w = 1; f0->b[0] = 1; f1->w = 1; f1->b[0] = 3; for ( i = 2; i <= m; i++ ) { /* f2 = t*f1+f0 */ _bshiftup2(f1,-1,f2); _addup2_destructive(f2,f0); /* cyclic change of the variables */ t = f0; f0 = f1; f1 = f2; f2 = t; } _copyup2(f1,f); return 0; break; default: return -1; break; } } /* * f = an irreducible trinomial or pentanomial of degree d 'after' f * return value : 0 --- exists * 1 --- does not exist (exhaustion) */ int _generate_irreducible_polynomial(UP2 f,int d) { int ret,i,j,k,nz,i0,j0,k0; int w; unsigned int *fd; /* * if f = x^d+x^i+1 then i0 <- i, j0 <- 0, k0 <-0. * if f = x^d+x^k+x^j+x^i+1 (k>j>i) then i0 <- i, j0 <- j, k0 <-k. * otherwise i0,j0,k0 is set to 0. */ fd = f->b; w = (d>>5)+1; if ( f->w && (d==degup2(f)) ) { for ( nz = 0, i = d; i >= 0; i-- ) if ( fd[i>>5]&(1<<(i&31)) ) nz++; switch ( nz ) { case 3: for ( i0 = 1; !(fd[i0>>5]&(1<<(i0&31))) ; i0++ ); /* reset i0-th bit */ fd[i0>>5] &= ~(1<<(i0&31)); j0 = k0 = 0; break; case 5: for ( i0 = 1; !(fd[i0>>5]&(1<<(i0&31))) ; i0++ ); /* reset i0-th bit */ fd[i0>>5] &= ~(1<<(i0&31)); for ( j0 = i0+1; !(fd[j0>>5]&(1<<(j0&31))) ; j0++ ); /* reset j0-th bit */ fd[j0>>5] &= ~(1<<(j0&31)); for ( k0 = j0+1; !(fd[k0>>5]&(1<<(k0&31))) ; k0++ ); /* reset k0-th bit */ fd[k0>>5] &= ~(1<<(k0&31)); break; default: f->w = 0; break; } } else f->w = 0; if ( !f->w ) { fd = f->b; f->w = w; fd[0] |= 1; fd[d>>5] |= (1<<(d&31)); i0 = j0 = k0 = 0; } /* if j0 > 0 then f is already a pentanomial */ if ( j0 > 0 ) goto PENTA; /* searching for an irreducible trinomial */ for ( i = 1; 2*i <= d; i++ ) { /* skip the polynomials 'before' f */ if ( i < i0 ) continue; if ( i == i0 ) { i0 = 0; continue; } /* set i-th bit */ fd[i>>5] |= (1<<(i&31)); ret = irredcheck_dddup2(f); if ( ret == 1 ) return 0; /* reset i-th bit */ fd[i>>5] &= ~(1<<(i&31)); } /* searching for an irreducible pentanomial */ PENTA: for ( i = 1; i < d; i++ ) { /* skip the polynomials 'before' f */ if ( i < i0 ) continue; if ( i == i0 ) i0 = 0; /* set i-th bit */ fd[i>>5] |= (1<<(i&31)); for ( j = i+1; j < d; j++ ) { /* skip the polynomials 'before' f */ if ( j < j0 ) continue; if ( j == j0 ) j0 = 0; /* set j-th bit */ fd[j>>5] |= (1<<(j&31)); for ( k = j+1; k < d; k++ ) { /* skip the polynomials 'before' f */ if ( k < k0 ) continue; else if ( k == k0 ) { k0 = 0; continue; } /* set k-th bit */ fd[k>>5] |= (1<<(k&31)); ret = irredcheck_dddup2(f); if ( ret == 1 ) return 0; /* reset k-th bit */ fd[k>>5] &= ~(1<<(k&31)); } /* reset j-th bit */ fd[j>>5] &= ~(1<<(j&31)); } /* reset i-th bit */ fd[i>>5] &= ~(1<<(i&31)); } /* exhausted */ return 1; } /* * f = an irreducible trinomial or pentanomial of degree d 'after' f * * searching strategy: * trinomial x^d+x^i+1: * i is as small as possible. * trinomial x^d+x^i+x^j+x^k+1: * i is as small as possible. * For such i, j is as small as possible. * For such i and j, 'k' is as small as possible. * * return value : 0 --- exists * 1 --- does not exist (exhaustion) */ int _generate_good_irreducible_polynomial(UP2 f,int d) { int ret,i,j,k,nz,i0,j0,k0; int w; unsigned int *fd; /* * if f = x^d+x^i+1 then i0 <- i, j0 <- 0, k0 <-0. * if f = x^d+x^k+x^j+x^i+1 (k>j>i) then i0 <- i, j0 <- j, k0 <-k. * otherwise i0,j0,k0 is set to 0. */ fd = f->b; w = (d>>5)+1; if ( f->w && (d==degup2(f)) ) { for ( nz = 0, i = d; i >= 0; i-- ) if ( fd[i>>5]&(1<<(i&31)) ) nz++; switch ( nz ) { case 3: for ( i0 = 1; !(fd[i0>>5]&(1<<(i0&31))) ; i0++ ); /* reset i0-th bit */ fd[i0>>5] &= ~(1<<(i0&31)); j0 = k0 = 0; break; case 5: for ( i0 = 1; !(fd[i0>>5]&(1<<(i0&31))) ; i0++ ); /* reset i0-th bit */ fd[i0>>5] &= ~(1<<(i0&31)); for ( j0 = i0+1; !(fd[j0>>5]&(1<<(j0&31))) ; j0++ ); /* reset j0-th bit */ fd[j0>>5] &= ~(1<<(j0&31)); for ( k0 = j0+1; !(fd[k0>>5]&(1<<(k0&31))) ; k0++ ); /* reset k0-th bit */ fd[k0>>5] &= ~(1<<(k0&31)); break; default: f->w = 0; break; } } else f->w = 0; if ( !f->w ) { fd = f->b; f->w = w; fd[0] |= 1; fd[d>>5] |= (1<<(d&31)); i0 = j0 = k0 = 0; } /* if j0 > 0 then f is already a pentanomial */ if ( j0 > 0 ) goto PENTA; /* searching for an irreducible trinomial */ for ( i = 1; 2*i <= d; i++ ) { /* skip the polynomials 'before' f */ if ( i < i0 ) continue; if ( i == i0 ) { i0 = 0; continue; } /* set i-th bit */ fd[i>>5] |= (1<<(i&31)); ret = irredcheck_dddup2(f); if ( ret == 1 ) return 0; /* reset i-th bit */ fd[i>>5] &= ~(1<<(i&31)); } /* searching for an irreducible pentanomial */ PENTA: for ( i = 3; i < d; i++ ) { /* skip the polynomials 'before' f */ if ( i < i0 ) continue; if ( i == i0 ) i0 = 0; /* set i-th bit */ fd[i>>5] |= (1<<(i&31)); for ( j = 2; j < i; j++ ) { /* skip the polynomials 'before' f */ if ( j < j0 ) continue; if ( j == j0 ) j0 = 0; /* set j-th bit */ fd[j>>5] |= (1<<(j&31)); for ( k = 1; k < j; k++ ) { /* skip the polynomials 'before' f */ if ( k < k0 ) continue; else if ( k == k0 ) { k0 = 0; continue; } /* set k-th bit */ fd[k>>5] |= (1<<(k&31)); ret = irredcheck_dddup2(f); if ( ret == 1 ) return 0; /* reset k-th bit */ fd[k>>5] &= ~(1<<(k&31)); } /* reset j-th bit */ fd[j>>5] &= ~(1<<(j&31)); } /* reset i-th bit */ fd[i>>5] &= ~(1<<(i&31)); } /* exhausted */ return 1; } printqmat(mat,row,col) Q **mat; int row,col; { int i,j; for ( i = 0; i < row; i++ ) { for ( j = 0; j < col; j++ ) { printnum(mat[i][j]); printf(" "); } printf("\n"); } } printimat(mat,row,col) int **mat; int row,col; { int i,j; for ( i = 0; i < row; i++ ) { for ( j = 0; j < col; j++ ) { printf("%d ",mat[i][j]); } printf("\n"); } }