![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to array.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib2 / asir2000 / builtin|
Return to array.c CVS log | Up to [local] / OpenXM_contrib2 / asir2000 / builtin |
1.4 ! noro 1: /* $OpenXM: OpenXM_contrib2/asir2000/builtin/array.c,v 1.3 2000/04/20 02:20:15 noro Exp $ */
1.1 noro 2: #include "ca.h"
3: #include "base.h"
4: #include "parse.h"
5: #include "inline.h"
1.4 ! noro 6:
! 7: #if 0
1.1 noro 8: #undef DMAR
9: #define DMAR(a1,a2,a3,d,r) (r)=dmar(a1,a2,a3,d);
1.4 ! noro 10: #endif
1.1 noro 11:
12: extern int Print; /* XXX */
13:
1.3 noro 14: void inner_product_mat_int_mod(Q **,int **,int,int,int,Q *);
15: void solve_by_lu_mod(int **,int,int,int **,int);
1.1 noro 16: void solve_by_lu_gfmmat(GFMMAT,unsigned int,unsigned int *,unsigned int *);
17: int lu_gfmmat(GFMMAT,unsigned int,int *);
18: void mat_to_gfmmat(MAT,unsigned int,GFMMAT *);
19:
20: int generic_gauss_elim_mod(int **,int,int,int,int *);
21: int generic_gauss_elim(MAT ,MAT *,Q *,int **,int **);
22:
23: int gauss_elim_mod(int **,int,int,int);
24: int gauss_elim_mod1(int **,int,int,int);
25: int gauss_elim_geninv_mod(unsigned int **,int,int,int);
26: int gauss_elim_geninv_mod_swap(unsigned int **,int,int,unsigned int,unsigned int ***,int **);
27: void Pnewvect(), Pnewmat(), Psepvect(), Psize(), Pdet(), Pleqm(), Pleqm1(), Pgeninvm();
28:
29: void Pgeneric_gauss_elim_mod();
30:
31: void Pmat_to_gfmmat(),Plu_gfmmat(),Psolve_by_lu_gfmmat();
32: void Pgeninvm_swap(), Premainder(), Psremainder(), Pvtol();
33: void sepvect();
34: void Pmulmat_gf2n();
35: void Pbconvmat_gf2n();
36: void Pmul_vect_mat_gf2n();
37: void PNBmul_gf2n();
38: void Pmul_mat_vect_int();
39: void Psepmat_destructive();
40: void Px962_irredpoly_up2();
41: void Pirredpoly_up2();
42: void Pnbpoly_up2();
43: void Pqsort();
44:
45: struct ftab array_tab[] = {
46: {"solve_by_lu_gfmmat",Psolve_by_lu_gfmmat,4},
47: {"lu_gfmmat",Plu_gfmmat,2},
48: {"mat_to_gfmmat",Pmat_to_gfmmat,2},
49: {"generic_gauss_elim_mod",Pgeneric_gauss_elim_mod,2},
50: {"newvect",Pnewvect,-2},
51: {"newmat",Pnewmat,-3},
52: {"sepmat_destructive",Psepmat_destructive,2},
53: {"sepvect",Psepvect,2},
54: {"qsort",Pqsort,-2},
55: {"vtol",Pvtol,1},
56: {"size",Psize,1},
57: {"det",Pdet,-2},
58: {"leqm",Pleqm,2},
59: {"leqm1",Pleqm1,2},
60: {"geninvm",Pgeninvm,2},
61: {"geninvm_swap",Pgeninvm_swap,2},
62: {"remainder",Premainder,2},
63: {"sremainder",Psremainder,2},
64: {"mulmat_gf2n",Pmulmat_gf2n,1},
65: {"bconvmat_gf2n",Pbconvmat_gf2n,-4},
66: {"mul_vect_mat_gf2n",Pmul_vect_mat_gf2n,2},
67: {"mul_mat_vect_int",Pmul_mat_vect_int,2},
68: {"nbmul_gf2n",PNBmul_gf2n,3},
69: {"x962_irredpoly_up2",Px962_irredpoly_up2,2},
70: {"irredpoly_up2",Pirredpoly_up2,2},
71: {"nbpoly_up2",Pnbpoly_up2,2},
72: {0,0,0},
73: };
74:
75: int comp_obj(a,b)
76: Obj *a,*b;
77: {
78: return arf_comp(CO,*a,*b);
79: }
80:
81: static FUNC generic_comp_obj_func;
82: static NODE generic_comp_obj_arg;
83:
84: int generic_comp_obj(a,b)
85: Obj *a,*b;
86: {
87: Q r;
88:
89: BDY(generic_comp_obj_arg)=(pointer)(*a);
90: BDY(NEXT(generic_comp_obj_arg))=(pointer)(*b);
91: r = (Q)bevalf(generic_comp_obj_func,generic_comp_obj_arg);
92: if ( !r )
93: return 0;
94: else
95: return SGN(r)>0?1:-1;
96: }
97:
98:
99: void Pqsort(arg,rp)
100: NODE arg;
101: VECT *rp;
102: {
103: VECT vect;
104: char buf[BUFSIZ];
105: char *fname;
106: NODE n;
107: P p;
108: V v;
109:
110: asir_assert(ARG0(arg),O_VECT,"qsort");
111: vect = (VECT)ARG0(arg);
112: if ( argc(arg) == 1 )
113: qsort(BDY(vect),vect->len,sizeof(Obj),(int (*)(const void *,const void *))comp_obj);
114: else {
115: p = (P)ARG1(arg);
116: if ( !p || OID(p)!=2 )
117: error("qsort : invalid argument");
118: v = VR(p);
119: if ( (int)v->attr != V_SR )
120: error("qsort : no such function");
121: generic_comp_obj_func = (FUNC)v->priv;
122: MKNODE(n,0,0); MKNODE(generic_comp_obj_arg,0,n);
123: qsort(BDY(vect),vect->len,sizeof(Obj),(int (*)(const void *,const void *))generic_comp_obj);
124: }
125: *rp = vect;
126: }
127:
128: void PNBmul_gf2n(arg,rp)
129: NODE arg;
130: GF2N *rp;
131: {
132: GF2N a,b;
133: GF2MAT mat;
134: int n,w;
135: unsigned int *ab,*bb;
136: UP2 r;
137:
138: a = (GF2N)ARG0(arg);
139: b = (GF2N)ARG1(arg);
140: mat = (GF2MAT)ARG2(arg);
141: if ( !a || !b )
142: *rp = 0;
143: else {
144: n = mat->row;
145: w = (n+BSH-1)/BSH;
146:
147: ab = (unsigned int *)ALLOCA(w*sizeof(unsigned int));
148: bzero((char *)ab,w*sizeof(unsigned int));
149: bcopy(a->body->b,ab,(a->body->w)*sizeof(unsigned int));
150:
151: bb = (unsigned int *)ALLOCA(w*sizeof(unsigned int));
152: bzero((char *)bb,w*sizeof(unsigned int));
153: bcopy(b->body->b,bb,(b->body->w)*sizeof(unsigned int));
154:
155: NEWUP2(r,w);
156: bzero((char *)r->b,w*sizeof(unsigned int));
157: mul_nb(mat,ab,bb,r->b);
158: r->w = w;
159: _adjup2(r);
160: if ( !r->w )
161: *rp = 0;
162: else
163: MKGF2N(r,*rp);
164: }
165: }
166:
167: void Pmul_vect_mat_gf2n(arg,rp)
168: NODE arg;
169: GF2N *rp;
170: {
171: GF2N a;
172: GF2MAT mat;
173: int n,w;
174: unsigned int *b;
175: UP2 r;
176:
177: a = (GF2N)ARG0(arg);
178: mat = (GF2MAT)ARG1(arg);
179: if ( !a )
180: *rp = 0;
181: else {
182: n = mat->row;
183: w = (n+BSH-1)/BSH;
184: b = (unsigned int *)ALLOCA(w*sizeof(unsigned int));
185: bzero((char *)b,w*sizeof(unsigned int));
186: bcopy(a->body->b,b,(a->body->w)*sizeof(unsigned int));
187: NEWUP2(r,w);
188: bzero((char *)r->b,w*sizeof(unsigned int));
189: mulgf2vectmat(mat->row,b,mat->body,r->b);
190: r->w = w;
191: _adjup2(r);
192: if ( !r->w )
193: *rp = 0;
194: else {
195: MKGF2N(r,*rp);
196: }
197: }
198: }
199:
200: void Pbconvmat_gf2n(arg,rp)
201: NODE arg;
202: LIST *rp;
203: {
204: P p0,p1;
205: int to;
206: GF2MAT p01,p10;
207: GF2N root;
208: NODE n0,n1;
209:
210: p0 = (P)ARG0(arg);
211: p1 = (P)ARG1(arg);
212: to = ARG2(arg)?1:0;
213: if ( argc(arg) == 4 ) {
214: root = (GF2N)ARG3(arg);
215: compute_change_of_basis_matrix_with_root(p0,p1,to,root,&p01,&p10);
216: } else
217: compute_change_of_basis_matrix(p0,p1,to,&p01,&p10);
218: MKNODE(n1,p10,0); MKNODE(n0,p01,n1);
219: MKLIST(*rp,n0);
220: }
221:
222: void Pmulmat_gf2n(arg,rp)
223: NODE arg;
224: GF2MAT *rp;
225: {
226: GF2MAT m;
227:
228: if ( !compute_multiplication_matrix((P)ARG0(arg),&m) )
229: error("mulmat_gf2n : input is not a normal polynomial");
230: *rp = m;
231: }
232:
233: void Psepmat_destructive(arg,rp)
234: NODE arg;
235: LIST *rp;
236: {
237: MAT mat,mat1;
238: int i,j,row,col;
239: Q **a,**a1;
240: Q ent;
241: N nm,mod,rem,quo;
242: int sgn;
243: NODE n0,n1;
244:
245: mat = (MAT)ARG0(arg); mod = NM((Q)ARG1(arg));
246: row = mat->row; col = mat->col;
247: MKMAT(mat1,row,col);
248: a = (Q **)mat->body; a1 = (Q **)mat1->body;
249: for ( i = 0; i < row; i++ )
250: for ( j = 0; j < col; j++ ) {
251: ent = a[i][j];
252: if ( !ent )
253: continue;
254: nm = NM(ent);
255: sgn = SGN(ent);
256: divn(nm,mod,&quo,&rem);
257: /* if ( quo != nm && rem != nm ) */
258: /* GC_free(nm); */
259: /* GC_free(ent); */
260: NTOQ(rem,sgn,a[i][j]); NTOQ(quo,sgn,a1[i][j]);
261: }
262: MKNODE(n1,mat1,0); MKNODE(n0,mat,n1);
263: MKLIST(*rp,n0);
264: }
265:
266: void Psepvect(arg,rp)
267: NODE arg;
268: VECT *rp;
269: {
270: sepvect((VECT)ARG0(arg),QTOS((Q)ARG1(arg)),rp);
271: }
272:
273: void sepvect(v,d,rp)
274: VECT v;
275: int d;
276: VECT *rp;
277: {
278: int i,j,k,n,q,q1,r;
279: pointer *pv,*pw,*pu;
280: VECT w,u;
281:
282: n = v->len;
283: if ( d > n )
284: d = n;
285: q = n/d; r = n%d; q1 = q+1;
286: MKVECT(w,d); *rp = w;
287: pv = BDY(v); pw = BDY(w); k = 0;
288: for ( i = 0; i < r; i++ ) {
289: MKVECT(u,q1); pw[i] = (pointer)u;
290: for ( pu = BDY(u), j = 0; j < q1; j++, k++ )
291: pu[j] = pv[k];
292: }
293: for ( ; i < d; i++ ) {
294: MKVECT(u,q); pw[i] = (pointer)u;
295: for ( pu = BDY(u), j = 0; j < q; j++, k++ )
296: pu[j] = pv[k];
297: }
298: }
299:
300: void Pnewvect(arg,rp)
301: NODE arg;
302: VECT *rp;
303: {
304: int len,i,r;
305: VECT vect;
306: pointer *vb;
307: LIST list;
308: NODE tn;
309:
310: asir_assert(ARG0(arg),O_N,"newvect");
311: len = QTOS((Q)ARG0(arg));
312: if ( len <= 0 )
313: error("newvect : invalid size");
314: MKVECT(vect,len);
315: if ( argc(arg) == 2 ) {
316: list = (LIST)ARG1(arg);
317: asir_assert(list,O_LIST,"newvect");
318: for ( r = 0, tn = BDY(list); tn; r++, tn = NEXT(tn) );
319: if ( r > len ) {
320: *rp = vect;
321: return;
322: }
323: for ( i = 0, tn = BDY(list), vb = BDY(vect); tn; i++, tn = NEXT(tn) )
324: vb[i] = (pointer)BDY(tn);
325: }
326: *rp = vect;
327: }
328:
329: void Pnewmat(arg,rp)
330: NODE arg;
331: MAT *rp;
332: {
333: int row,col;
334: int i,j,r,c;
335: NODE tn,sn;
336: MAT m;
337: pointer **mb;
338: LIST list;
339:
340: asir_assert(ARG0(arg),O_N,"newmat");
341: asir_assert(ARG1(arg),O_N,"newmat");
342: row = QTOS((Q)ARG0(arg)); col = QTOS((Q)ARG1(arg));
343: if ( row <= 0 || col <= 0 )
344: error("newmat : invalid size");
345: MKMAT(m,row,col);
346: if ( argc(arg) == 3 ) {
347: list = (LIST)ARG2(arg);
348: asir_assert(list,O_LIST,"newmat");
349: for ( r = 0, c = 0, tn = BDY(list); tn; r++, tn = NEXT(tn) ) {
350: for ( j = 0, sn = BDY((LIST)BDY(tn)); sn; j++, sn = NEXT(sn) );
351: c = MAX(c,j);
352: }
353: if ( (r > row) || (c > col) ) {
354: *rp = m;
355: return;
356: }
357: for ( i = 0, tn = BDY(list), mb = BDY(m); tn; i++, tn = NEXT(tn) ) {
358: asir_assert(BDY(tn),O_LIST,"newmat");
359: for ( j = 0, sn = BDY((LIST)BDY(tn)); sn; j++, sn = NEXT(sn) )
360: mb[i][j] = (pointer)BDY(sn);
361: }
362: }
363: *rp = m;
364: }
365:
366: void Pvtol(arg,rp)
367: NODE arg;
368: LIST *rp;
369: {
370: NODE n,n1;
371: VECT v;
372: pointer *a;
373: int len,i;
374:
375: asir_assert(ARG0(arg),O_VECT,"vtol");
376: v = (VECT)ARG0(arg); len = v->len; a = BDY(v);
377: for ( i = len - 1, n = 0; i >= 0; i-- ) {
378: MKNODE(n1,a[i],n); n = n1;
379: }
380: MKLIST(*rp,n);
381: }
382:
383: void Premainder(arg,rp)
384: NODE arg;
385: Obj *rp;
386: {
387: Obj a;
388: VECT v,w;
389: MAT m,l;
390: pointer *vb,*wb;
391: pointer **mb,**lb;
392: int id,i,j,n,row,col,t,smd,sgn;
393: Q md,q;
394:
395: a = (Obj)ARG0(arg); md = (Q)ARG1(arg);
396: if ( !a )
397: *rp = 0;
398: else {
399: id = OID(a);
400: switch ( id ) {
401: case O_N:
402: case O_P:
403: cmp(md,(P)a,(P *)rp); break;
404: case O_VECT:
405: smd = QTOS(md);
406: v = (VECT)a; n = v->len; vb = v->body;
407: MKVECT(w,n); wb = w->body;
408: for ( i = 0; i < n; i++ ) {
409: if ( q = (Q)vb[i] ) {
410: sgn = SGN(q); t = rem(NM(q),smd);
411: STOQ(t,q);
412: if ( q )
413: SGN(q) = sgn;
414: }
415: wb[i] = (pointer)q;
416: }
417: *rp = (Obj)w;
418: break;
419: case O_MAT:
420: m = (MAT)a; row = m->row; col = m->col; mb = m->body;
421: MKMAT(l,row,col); lb = l->body;
422: for ( i = 0; i < row; i++ )
423: for ( j = 0, vb = mb[i], wb = lb[i]; j < col; j++ )
424: cmp(md,(P)vb[j],(P *)&wb[j]);
425: *rp = (Obj)l;
426: break;
427: default:
428: error("remainder : invalid argument");
429: }
430: }
431: }
432:
433: void Psremainder(arg,rp)
434: NODE arg;
435: Obj *rp;
436: {
437: Obj a;
438: VECT v,w;
439: MAT m,l;
440: pointer *vb,*wb;
441: pointer **mb,**lb;
442: unsigned int t,smd;
443: int id,i,j,n,row,col;
444: Q md,q;
445:
446: a = (Obj)ARG0(arg); md = (Q)ARG1(arg);
447: if ( !a )
448: *rp = 0;
449: else {
450: id = OID(a);
451: switch ( id ) {
452: case O_N:
453: case O_P:
454: cmp(md,(P)a,(P *)rp); break;
455: case O_VECT:
456: smd = QTOS(md);
457: v = (VECT)a; n = v->len; vb = v->body;
458: MKVECT(w,n); wb = w->body;
459: for ( i = 0; i < n; i++ ) {
460: if ( q = (Q)vb[i] ) {
461: t = (unsigned int)rem(NM(q),smd);
462: if ( SGN(q) < 0 )
463: t = (smd - t) % smd;
464: UTOQ(t,q);
465: }
466: wb[i] = (pointer)q;
467: }
468: *rp = (Obj)w;
469: break;
470: case O_MAT:
471: m = (MAT)a; row = m->row; col = m->col; mb = m->body;
472: MKMAT(l,row,col); lb = l->body;
473: for ( i = 0; i < row; i++ )
474: for ( j = 0, vb = mb[i], wb = lb[i]; j < col; j++ )
475: cmp(md,(P)vb[j],(P *)&wb[j]);
476: *rp = (Obj)l;
477: break;
478: default:
479: error("remainder : invalid argument");
480: }
481: }
482: }
483:
484: void Psize(arg,rp)
485: NODE arg;
486: LIST *rp;
487: {
488:
489: int n,m;
490: Q q;
491: NODE t,s;
492:
493: if ( !ARG0(arg) )
494: t = 0;
495: else {
496: switch (OID(ARG0(arg))) {
497: case O_VECT:
498: n = ((VECT)ARG0(arg))->len;
499: STOQ(n,q); MKNODE(t,q,0);
500: break;
501: case O_MAT:
502: n = ((MAT)ARG0(arg))->row; m = ((MAT)ARG0(arg))->col;
503: STOQ(m,q); MKNODE(s,q,0); STOQ(n,q); MKNODE(t,q,s);
504: break;
505: default:
506: error("size : invalid argument"); break;
507: }
508: }
509: MKLIST(*rp,t);
510: }
511:
512: void Pdet(arg,rp)
513: NODE arg;
514: P *rp;
515: {
516: MAT m;
517: int n,i,j,mod;
518: P d;
519: P **mat,**w;
520:
521: m = (MAT)ARG0(arg);
522: asir_assert(m,O_MAT,"det");
523: if ( m->row != m->col )
524: error("det : non-square matrix");
525: else if ( argc(arg) == 1 )
526: detp(CO,(P **)BDY(m),m->row,rp);
527: else {
528: n = m->row; mod = QTOS((Q)ARG1(arg)); mat = (P **)BDY(m);
529: w = (P **)almat_pointer(n,n);
530: for ( i = 0; i < n; i++ )
531: for ( j = 0; j < n; j++ )
532: ptomp(mod,mat[i][j],&w[i][j]);
533: detmp(CO,mod,w,n,&d);
534: mptop(d,rp);
535: }
536: }
537:
538: /*
539: input : a row x col matrix A
540: A[I] <-> A[I][0]*x_0+A[I][1]*x_1+...
541:
542: output : [B,R,C]
543: B : a rank(A) x col-rank(A) matrix
544: R : a vector of length rank(A)
545: C : a vector of length col-rank(A)
546: B[I] <-> x_{R[I]}+B[I][0]x_{C[0]}+B[I][1]x_{C[1]}+...
547: */
548:
549: void Pgeneric_gauss_elim_mod(arg,rp)
550: NODE arg;
551: LIST *rp;
552: {
553: NODE n0;
554: MAT m,mat;
555: VECT rind,cind;
556: Q **tmat;
557: int **wmat;
558: Q *rib,*cib;
559: int *colstat;
560: Q q;
561: int md,i,j,k,l,row,col,t,n,rank;
562:
563: asir_assert(ARG0(arg),O_MAT,"generic_gauss_elim_mod");
564: asir_assert(ARG1(arg),O_N,"generic_gauss_elim_mod");
565: m = (MAT)ARG0(arg); md = QTOS((Q)ARG1(arg));
566: row = m->row; col = m->col; tmat = (Q **)m->body;
567: wmat = (int **)almat(row,col);
568: colstat = (int *)MALLOC_ATOMIC(col*sizeof(int));
569: for ( i = 0; i < row; i++ )
570: for ( j = 0; j < col; j++ )
571: if ( q = (Q)tmat[i][j] ) {
572: t = rem(NM(q),md);
573: if ( t && SGN(q) < 0 )
574: t = (md - t) % md;
575: wmat[i][j] = t;
576: } else
577: wmat[i][j] = 0;
578: rank = generic_gauss_elim_mod(wmat,row,col,md,colstat);
579:
580: MKMAT(mat,rank,col-rank);
581: tmat = (Q **)mat->body;
582: for ( i = 0; i < rank; i++ )
583: for ( j = k = 0; j < col; j++ )
584: if ( !colstat[j] ) {
585: UTOQ(wmat[i][j],tmat[i][k]); k++;
586: }
587:
588: MKVECT(rind,rank);
589: MKVECT(cind,col-rank);
590: rib = (Q *)rind->body; cib = (Q *)cind->body;
591: for ( j = k = l = 0; j < col; j++ )
592: if ( colstat[j] ) {
593: STOQ(j,rib[k]); k++;
594: } else {
595: STOQ(j,cib[l]); l++;
596: }
597: n0 = mknode(3,mat,rind,cind);
598: MKLIST(*rp,n0);
599: }
600:
601: void Pleqm(arg,rp)
602: NODE arg;
603: VECT *rp;
604: {
605: MAT m;
606: VECT vect;
607: pointer **mat;
608: Q *v;
609: Q q;
610: int **wmat;
611: int md,i,j,row,col,t,n,status;
612:
613: asir_assert(ARG0(arg),O_MAT,"leqm");
614: asir_assert(ARG1(arg),O_N,"leqm");
615: m = (MAT)ARG0(arg); md = QTOS((Q)ARG1(arg));
616: row = m->row; col = m->col; mat = m->body;
617: wmat = (int **)almat(row,col);
618: for ( i = 0; i < row; i++ )
619: for ( j = 0; j < col; j++ )
620: if ( q = (Q)mat[i][j] ) {
621: t = rem(NM(q),md);
622: if ( SGN(q) < 0 )
623: t = (md - t) % md;
624: wmat[i][j] = t;
625: } else
626: wmat[i][j] = 0;
627: status = gauss_elim_mod(wmat,row,col,md);
628: if ( status < 0 )
629: *rp = 0;
630: else if ( status > 0 )
631: *rp = (VECT)ONE;
632: else {
633: n = col - 1;
634: MKVECT(vect,n);
635: for ( i = 0, v = (Q *)vect->body; i < n; i++ ) {
636: t = (md-wmat[i][n])%md; STOQ(t,v[i]);
637: }
638: *rp = vect;
639: }
640: }
641:
642: int gauss_elim_mod(mat,row,col,md)
643: int **mat;
644: int row,col,md;
645: {
646: int i,j,k,inv,a,n;
647: int *t,*pivot;
648:
649: n = col - 1;
650: for ( j = 0; j < n; j++ ) {
651: for ( i = j; i < row && !mat[i][j]; i++ );
652: if ( i == row )
653: return 1;
654: if ( i != j ) {
655: t = mat[i]; mat[i] = mat[j]; mat[j] = t;
656: }
657: pivot = mat[j];
658: inv = invm(pivot[j],md);
659: for ( k = j; k <= n; k++ ) {
660: /* pivot[k] = dmar(pivot[k],inv,0,md); */
661: DMAR(pivot[k],inv,0,md,pivot[k])
662: }
663: for ( i = 0; i < row; i++ ) {
664: t = mat[i];
665: if ( i != j && (a = t[j]) )
666: for ( k = j, a = md - a; k <= n; k++ ) {
667: /* t[k] = dmar(pivot[k],a,t[k],md); */
668: DMAR(pivot[k],a,t[k],md,t[k])
669: }
670: }
671: }
672: for ( i = n; i < row && !mat[i][n]; i++ );
673: if ( i == row )
674: return 0;
675: else
676: return -1;
677: }
678:
1.4 ! noro 679: struct oEGT eg_mod,eg_elim,eg_elim1,eg_elim2,eg_chrem,eg_gschk,eg_intrat,eg_symb;
1.1 noro 680:
681: int generic_gauss_elim(mat,nm,dn,rindp,cindp)
682: MAT mat;
683: MAT *nm;
684: Q *dn;
685: int **rindp,**cindp;
686: {
687: int **wmat;
688: Q **bmat;
689: N **tmat;
690: Q *bmi;
691: N *tmi;
692: Q q;
693: int *wmi;
694: int *colstat,*wcolstat,*rind,*cind;
695: int row,col,ind,md,i,j,k,l,t,t1,rank,rank0,inv;
696: N m1,m2,m3,s,u;
697: MAT r,crmat;
698: struct oEGT tmp0,tmp1;
699: struct oEGT eg_mod_split,eg_elim_split,eg_chrem_split;
700: struct oEGT eg_intrat_split,eg_gschk_split;
701: int ret;
702:
703: init_eg(&eg_mod_split); init_eg(&eg_chrem_split);
704: init_eg(&eg_elim_split); init_eg(&eg_intrat_split);
705: init_eg(&eg_gschk_split);
706: bmat = (Q **)mat->body;
707: row = mat->row; col = mat->col;
708: wmat = (int **)almat(row,col);
709: colstat = (int *)MALLOC_ATOMIC(col*sizeof(int));
710: wcolstat = (int *)MALLOC_ATOMIC(col*sizeof(int));
711: for ( ind = 0; ; ind++ ) {
1.2 noro 712: if ( Print ) {
713: fprintf(asir_out,"."); fflush(asir_out);
714: }
1.1 noro 715: md = lprime[ind];
716: get_eg(&tmp0);
717: for ( i = 0; i < row; i++ )
718: for ( j = 0, bmi = bmat[i], wmi = wmat[i]; j < col; j++ )
719: if ( q = (Q)bmi[j] ) {
720: t = rem(NM(q),md);
721: if ( t && SGN(q) < 0 )
722: t = (md - t) % md;
723: wmi[j] = t;
724: } else
725: wmi[j] = 0;
726: get_eg(&tmp1);
727: add_eg(&eg_mod,&tmp0,&tmp1);
728: add_eg(&eg_mod_split,&tmp0,&tmp1);
729: get_eg(&tmp0);
730: rank = generic_gauss_elim_mod(wmat,row,col,md,wcolstat);
731: get_eg(&tmp1);
732: add_eg(&eg_elim,&tmp0,&tmp1);
733: add_eg(&eg_elim_split,&tmp0,&tmp1);
734: if ( !ind ) {
735: RESET:
736: UTON(md,m1);
737: rank0 = rank;
738: bcopy(wcolstat,colstat,col*sizeof(int));
739: MKMAT(crmat,rank,col-rank);
740: MKMAT(r,rank,col-rank); *nm = r;
741: tmat = (N **)crmat->body;
742: for ( i = 0; i < rank; i++ )
743: for ( j = k = 0, tmi = tmat[i], wmi = wmat[i]; j < col; j++ )
744: if ( !colstat[j] ) {
745: UTON(wmi[j],tmi[k]); k++;
746: }
747: } else {
748: if ( rank < rank0 ) {
1.2 noro 749: if ( Print ) {
1.1 noro 750: fprintf(asir_out,"lower rank matrix; continuing...\n");
1.2 noro 751: fflush(asir_out);
752: }
1.1 noro 753: continue;
754: } else if ( rank > rank0 ) {
1.2 noro 755: if ( Print ) {
1.1 noro 756: fprintf(asir_out,"higher rank matrix; resetting...\n");
1.2 noro 757: fflush(asir_out);
758: }
1.1 noro 759: goto RESET;
760: } else {
761: for ( j = 0; (j<col) && (colstat[j]==wcolstat[j]); j++ );
762: if ( j < col ) {
1.2 noro 763: if ( Print ) {
1.1 noro 764: fprintf(asir_out,"inconsitent colstat; resetting...\n");
1.2 noro 765: fflush(asir_out);
766: }
1.1 noro 767: goto RESET;
768: }
769: }
770:
771: get_eg(&tmp0);
772: inv = invm(rem(m1,md),md);
773: UTON(md,m2); muln(m1,m2,&m3);
774: for ( i = 0; i < rank; i++ )
775: for ( j = k = 0, tmi = tmat[i], wmi = wmat[i]; j < col; j++ )
776: if ( !colstat[j] ) {
777: if ( tmi[k] ) {
778: /* f3 = f1+m1*(m1 mod m2)^(-1)*(f2 - f1 mod m2) */
779: t = rem(tmi[k],md);
780: if ( wmi[j] >= t )
781: t = wmi[j]-t;
782: else
783: t = md-(t-wmi[j]);
784: DMAR(t,inv,0,md,t1)
785: UTON(t1,u);
786: muln(m1,u,&s);
787: addn(tmi[k],s,&u); tmi[k] = u;
788: } else if ( wmi[j] ) {
789: /* f3 = m1*(m1 mod m2)^(-1)*f2 */
790: DMAR(wmi[j],inv,0,md,t)
791: UTON(t,u);
792: muln(m1,u,&s); tmi[k] = s;
793: }
794: k++;
795: }
796: m1 = m3;
797: get_eg(&tmp1);
798: add_eg(&eg_chrem,&tmp0,&tmp1);
799: add_eg(&eg_chrem_split,&tmp0,&tmp1);
800:
801: get_eg(&tmp0);
802: ret = intmtoratm(crmat,m1,*nm,dn);
803: get_eg(&tmp1);
804: add_eg(&eg_intrat,&tmp0,&tmp1);
805: add_eg(&eg_intrat_split,&tmp0,&tmp1);
806: if ( ret ) {
807: *rindp = rind = (int *)MALLOC_ATOMIC(rank*sizeof(int));
808: *cindp = cind = (int *)MALLOC_ATOMIC((col-rank)*sizeof(int));
809: for ( j = k = l = 0; j < col; j++ )
810: if ( colstat[j] )
811: rind[k++] = j;
812: else
813: cind[l++] = j;
814: get_eg(&tmp0);
1.3 noro 815: if ( gensolve_check(mat,*nm,*dn,rind,cind) ) {
816: get_eg(&tmp1);
817: add_eg(&eg_gschk,&tmp0,&tmp1);
818: add_eg(&eg_gschk_split,&tmp0,&tmp1);
819: if ( Print ) {
820: print_eg("Mod",&eg_mod_split);
821: print_eg("Elim",&eg_elim_split);
822: print_eg("ChRem",&eg_chrem_split);
823: print_eg("IntRat",&eg_intrat_split);
824: print_eg("Check",&eg_gschk_split);
825: fflush(asir_out);
826: }
827: return rank;
828: }
829: }
830: }
831: }
832: }
833:
834: int generic_gauss_elim_hensel(mat,nmmat,dn,rindp,cindp)
835: MAT mat;
836: MAT *nmmat;
837: Q *dn;
838: int **rindp,**cindp;
839: {
840: MAT bmat,xmat;
841: Q **a0,**a,**b,**x,**nm;
842: Q *ai,*bi,*xi;
843: int row,col;
844: int **w;
845: int *wi;
846: int **wc;
847: Q mdq,q,s,u;
848: N tn;
849: int ind,md,i,j,k,l,li,ri,rank;
850: unsigned int t;
851: int *cinfo,*rinfo;
852: int *rind,*cind;
853: int count;
854: struct oEGT eg_mul,eg_inv,tmp0,tmp1;
855:
856: a0 = (Q **)mat->body;
857: row = mat->row; col = mat->col;
858: w = (int **)almat(row,col);
859: for ( ind = 0; ; ind++ ) {
860: md = lprime[ind];
861: STOQ(md,mdq);
862: for ( i = 0; i < row; i++ )
863: for ( j = 0, ai = a0[i], wi = w[i]; j < col; j++ )
864: if ( q = (Q)ai[j] ) {
865: t = rem(NM(q),md);
866: if ( t && SGN(q) < 0 )
867: t = (md - t) % md;
868: wi[j] = t;
869: } else
870: wi[j] = 0;
871:
872: rank = find_lhs_and_lu_mod(w,row,col,md,&rinfo,&cinfo);
873: a = (Q **)almat_pointer(rank,rank); /* lhs mat */
874: MKMAT(bmat,rank,col-rank); b = (Q **)bmat->body; /* lhs mat */
875: for ( j = li = ri = 0; j < col; j++ )
876: if ( cinfo[j] ) {
877: /* the column is in lhs */
878: for ( i = 0; i < rank; i++ ) {
879: w[i][li] = w[i][j];
880: a[i][li] = a0[rinfo[i]][j];
881: }
882: li++;
883: } else {
884: /* the column is in rhs */
885: for ( i = 0; i < rank; i++ )
886: b[i][ri] = a0[rinfo[i]][j];
887: ri++;
888: }
889:
890: /* solve Ax+B=0; A: rank x rank, B: rank x ri */
891: MKMAT(xmat,rank,ri); x = (Q **)(xmat)->body;
892: MKMAT(*nmmat,rank,ri); nm = (Q **)(*nmmat)->body;
893: /* use the right part of w as work area */
894: /* ri = col - rank */
895: wc = (int **)almat(rank,ri);
896: for ( i = 0; i < rank; i++ )
897: wc[i] = w[i]+rank;
898: *rindp = rind = (int *)MALLOC_ATOMIC(rank*sizeof(int));
899: *cindp = cind = (int *)MALLOC_ATOMIC((ri)*sizeof(int));
900:
901: init_eg(&eg_mul); init_eg(&eg_inv);
902: for ( q = ONE, count = 0; ; count++ ) {
903: fprintf(stderr,".");
904: /* wc = -b mod md */
905: for ( i = 0; i < rank; i++ )
906: for ( j = 0, bi = b[i], wi = wc[i]; j < ri; j++ )
907: if ( u = (Q)bi[j] ) {
908: t = rem(NM(u),md);
909: if ( t && SGN(u) > 0 )
910: t = (md - t) % md;
911: wi[j] = t;
912: } else
913: wi[j] = 0;
914: /* wc = A^(-1)wc; wc is normalized */
915: get_eg(&tmp0);
916: solve_by_lu_mod(w,rank,md,wc,ri);
1.1 noro 917: get_eg(&tmp1);
1.3 noro 918: add_eg(&eg_inv,&tmp0,&tmp1);
919: /* x = x-q*wc */
920: for ( i = 0; i < rank; i++ )
921: for ( j = 0, xi = x[i], wi = wc[i]; j < ri; j++ ) {
922: STOQ(wi[j],u); mulq(q,u,&s);
923: subq(xi[j],s,&u); xi[j] = u;
924: }
925: get_eg(&tmp0);
926: for ( i = 0; i < rank; i++ )
927: for ( j = 0; j < ri; j++ ) {
928: inner_product_mat_int_mod(a,wc,rank,i,j,&u);
929: addq(b[i][j],u,&s);
930: if ( s ) {
931: t = divin(NM(s),md,&tn);
932: if ( t )
933: error("generic_gauss_elim_hensel:incosistent");
934: NTOQ(tn,SGN(s),b[i][j]);
935: } else
936: b[i][j] = 0;
937: }
938: get_eg(&tmp1);
939: add_eg(&eg_mul,&tmp0,&tmp1);
940: /* q = q*md */
941: mulq(q,mdq,&u); q = u;
942: if ( !(count % 2) && intmtoratm_q(xmat,NM(q),*nmmat,dn) ) {
943: for ( j = k = l = 0; j < col; j++ )
944: if ( cinfo[j] )
945: rind[k++] = j;
946: else
947: cind[l++] = j;
948: if ( gensolve_check(mat,*nmmat,*dn,rind,cind) ) {
949: fprintf(stderr,"\n");
950: print_eg("INV",&eg_inv);
951: print_eg("MUL",&eg_mul);
952: fflush(asir_out);
953: return rank;
954: }
1.1 noro 955: }
956: }
957: }
958: }
959:
960: int f4_nocheck;
961:
962: int gensolve_check(mat,nm,dn,rind,cind)
963: MAT mat,nm;
964: Q dn;
965: int *rind,*cind;
966: {
967: int row,col,rank,clen,i,j,k,l;
968: Q s,t,u;
969: Q *w;
970: Q *mati,*nmk;
971:
972: if ( f4_nocheck )
973: return 1;
974: row = mat->row; col = mat->col;
975: rank = nm->row; clen = nm->col;
976: w = (Q *)MALLOC(clen*sizeof(Q));
977: for ( i = 0; i < row; i++ ) {
978: mati = (Q *)mat->body[i];
979: #if 1
980: bzero(w,clen*sizeof(Q));
981: for ( k = 0; k < rank; k++ )
982: for ( l = 0, nmk = (Q *)nm->body[k]; l < clen; l++ ) {
983: mulq(mati[rind[k]],nmk[l],&t);
984: addq(w[l],t,&s); w[l] = s;
985: }
986: for ( j = 0; j < clen; j++ ) {
987: mulq(dn,mati[cind[j]],&t);
988: if ( cmpq(w[j],t) )
989: break;
990: }
991: #else
992: for ( j = 0; j < clen; j++ ) {
993: for ( k = 0, s = 0; k < rank; k++ ) {
994: mulq(mati[rind[k]],nm->body[k][j],&t);
995: addq(s,t,&u); s = u;
996: }
997: mulq(dn,mati[cind[j]],&t);
998: if ( cmpq(s,t) )
999: break;
1000: }
1001: #endif
1002: if ( j != clen )
1003: break;
1004: }
1005: if ( i != row )
1006: return 0;
1007: else
1008: return 1;
1009: }
1010:
1011: /* assuming 0 < c < m */
1012:
1013: int inttorat(c,m,b,sgnp,nmp,dnp)
1014: N c,m,b;
1015: int *sgnp;
1016: N *nmp,*dnp;
1017: {
1018: Q qq,t,u1,v1,r1,nm;
1019: N q,r,u2,v2,r2;
1020:
1021: u1 = 0; v1 = ONE; u2 = m; v2 = c;
1022: while ( cmpn(v2,b) >= 0 ) {
1023: divn(u2,v2,&q,&r2); u2 = v2; v2 = r2;
1024: NTOQ(q,1,qq); mulq(qq,v1,&t); subq(u1,t,&r1); u1 = v1; v1 = r1;
1025: }
1026: if ( cmpn(NM(v1),b) >= 0 )
1027: return 0;
1028: else {
1029: *nmp = v2;
1030: *dnp = NM(v1);
1031: *sgnp = SGN(v1);
1032: return 1;
1033: }
1034: }
1035:
1036: /* mat->body = N ** */
1037:
1038: int intmtoratm(mat,md,nm,dn)
1039: MAT mat;
1040: N md;
1041: MAT nm;
1042: Q *dn;
1043: {
1044: N t,s,b;
1045: Q bound,dn0,dn1,nm1,q,tq;
1046: int i,j,k,l,row,col;
1047: Q **rmat;
1048: N **tmat;
1049: N *tmi;
1050: Q *nmk;
1051: N u,unm,udn;
1052: int sgn,ret;
1053:
1.3 noro 1054: if ( UNIN(md) )
1055: return 0;
1.1 noro 1056: row = mat->row; col = mat->col;
1057: bshiftn(md,1,&t);
1058: isqrt(t,&s);
1059: bshiftn(s,64,&b);
1060: if ( !b )
1061: b = ONEN;
1062: dn0 = ONE;
1063: tmat = (N **)mat->body;
1064: rmat = (Q **)nm->body;
1065: for ( i = 0; i < row; i++ )
1066: for ( j = 0, tmi = tmat[i]; j < col; j++ )
1067: if ( tmi[j] ) {
1068: muln(tmi[j],NM(dn0),&s);
1069: remn(s,md,&u);
1070: ret = inttorat(u,md,b,&sgn,&unm,&udn);
1071: if ( !ret )
1072: return 0;
1073: else {
1074: NTOQ(unm,sgn,nm1);
1075: NTOQ(udn,1,dn1);
1076: if ( !UNIQ(dn1) ) {
1077: for ( k = 0; k < i; k++ )
1078: for ( l = 0, nmk = rmat[k]; l < col; l++ ) {
1079: mulq(nmk[l],dn1,&q); nmk[l] = q;
1080: }
1081: for ( l = 0, nmk = rmat[i]; l < j; l++ ) {
1082: mulq(nmk[l],dn1,&q); nmk[l] = q;
1083: }
1084: }
1085: rmat[i][j] = nm1;
1086: mulq(dn0,dn1,&q); dn0 = q;
1087: }
1088: }
1089: *dn = dn0;
1090: return 1;
1091: }
1092:
1.3 noro 1093: /* mat->body = Q ** */
1094:
1095: int intmtoratm_q(mat,md,nm,dn)
1096: MAT mat;
1097: N md;
1098: MAT nm;
1099: Q *dn;
1100: {
1101: N t,s,b;
1102: Q bound,dn0,dn1,nm1,q,tq;
1103: int i,j,k,l,row,col;
1104: Q **rmat;
1105: Q **tmat;
1106: Q *tmi;
1107: Q *nmk;
1108: N u,unm,udn;
1109: int sgn,ret;
1110:
1111: if ( UNIN(md) )
1112: return 0;
1113: row = mat->row; col = mat->col;
1114: bshiftn(md,1,&t);
1115: isqrt(t,&s);
1116: bshiftn(s,64,&b);
1117: if ( !b )
1118: b = ONEN;
1119: dn0 = ONE;
1120: tmat = (Q **)mat->body;
1121: rmat = (Q **)nm->body;
1122: for ( i = 0; i < row; i++ )
1123: for ( j = 0, tmi = tmat[i]; j < col; j++ )
1124: if ( tmi[j] ) {
1125: muln(NM(tmi[j]),NM(dn0),&s);
1126: remn(s,md,&u);
1127: ret = inttorat(u,md,b,&sgn,&unm,&udn);
1128: if ( !ret )
1129: return 0;
1130: else {
1131: if ( SGN(tmi[j])<0 )
1132: sgn = -sgn;
1133: NTOQ(unm,sgn,nm1);
1134: NTOQ(udn,1,dn1);
1135: if ( !UNIQ(dn1) ) {
1136: for ( k = 0; k < i; k++ )
1137: for ( l = 0, nmk = rmat[k]; l < col; l++ ) {
1138: mulq(nmk[l],dn1,&q); nmk[l] = q;
1139: }
1140: for ( l = 0, nmk = rmat[i]; l < j; l++ ) {
1141: mulq(nmk[l],dn1,&q); nmk[l] = q;
1142: }
1143: }
1144: rmat[i][j] = nm1;
1145: mulq(dn0,dn1,&q); dn0 = q;
1146: }
1147: }
1148: *dn = dn0;
1149: return 1;
1150: }
1151:
1.4 ! noro 1152: #define ONE_STEP1 if ( zzz = *s ) { DMAR(zzz,hc,*tj,md,*tj) } tj++; s++;
! 1153:
! 1154: void reduce_reducers_mod(mat,row,col,md)
! 1155: int **mat;
! 1156: int row,col;
! 1157: int md;
! 1158: {
! 1159: int i,j,k,l,hc,zzz;
! 1160: int *t,*s,*tj,*ind;
! 1161:
! 1162: /* reduce the reducers */
! 1163: ind = (int *)ALLOCA(row*sizeof(int));
! 1164: for ( i = 0; i < row; i++ ) {
! 1165: t = mat[i];
! 1166: for ( j = 0; j < col && !t[j]; j++ );
! 1167: /* register the position of the head term */
! 1168: ind[i] = j;
! 1169: for ( l = i-1; l >= 0; l-- ) {
! 1170: /* reduce mat[i] by mat[l] */
! 1171: if ( hc = t[ind[l]] ) {
! 1172: /* mat[i] = mat[i]-hc*mat[l] */
! 1173: j = ind[l];
! 1174: s = mat[l]+j;
! 1175: tj = t+j;
! 1176: hc = md-hc;
! 1177: k = col-j;
! 1178: for ( ; k >= 64; k -= 64 ) {
! 1179: ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
! 1180: ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
! 1181: ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
! 1182: ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
! 1183: ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
! 1184: ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
! 1185: ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
! 1186: ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
! 1187: ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
! 1188: ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
! 1189: ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
! 1190: ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
! 1191: ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
! 1192: ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
! 1193: ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
! 1194: ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
! 1195: }
! 1196: for ( ; k >= 0; k-- ) {
! 1197: if ( zzz = *s ) { DMAR(zzz,hc,*tj,md,*tj) } tj++; s++;
! 1198: }
! 1199: }
! 1200: }
! 1201: }
! 1202: }
! 1203:
! 1204: /*
! 1205: mat[i] : reducers (i=0,...,nred-1)
! 1206: spolys (i=nred,...,row-1)
! 1207: mat[0] < mat[1] < ... < mat[nred-1] w.r.t the term order
! 1208: 1. reduce the reducers
! 1209: 2. reduce spolys by the reduced reducers
! 1210: */
! 1211:
! 1212: void pre_reduce_mod(mat,row,col,nred,md)
! 1213: int **mat;
! 1214: int row,col,nred;
! 1215: int md;
! 1216: {
! 1217: int i,j,k,l,hc,inv;
! 1218: int *t,*s,*tk,*ind;
! 1219:
! 1220: #if 1
! 1221: /* reduce the reducers */
! 1222: ind = (int *)ALLOCA(row*sizeof(int));
! 1223: for ( i = 0; i < nred; i++ ) {
! 1224: /* make mat[i] monic and mat[i] by mat[0],...,mat[i-1] */
! 1225: t = mat[i];
! 1226: for ( j = 0; j < col && !t[j]; j++ );
! 1227: /* register the position of the head term */
! 1228: ind[i] = j;
! 1229: inv = invm(t[j],md);
! 1230: for ( k = j; k < col; k++ )
! 1231: if ( t[k] )
! 1232: DMAR(t[k],inv,0,md,t[k])
! 1233: for ( l = i-1; l >= 0; l-- ) {
! 1234: /* reduce mat[i] by mat[l] */
! 1235: if ( hc = t[ind[l]] ) {
! 1236: /* mat[i] = mat[i]-hc*mat[l] */
! 1237: for ( k = ind[l], hc = md-hc, s = mat[l]+k, tk = t+k;
! 1238: k < col; k++, tk++, s++ )
! 1239: if ( *s )
! 1240: DMAR(*s,hc,*tk,md,*tk)
! 1241: }
! 1242: }
! 1243: }
! 1244: /* reduce the spolys */
! 1245: for ( i = nred; i < row; i++ ) {
! 1246: t = mat[i];
! 1247: for ( l = nred-1; l >= 0; l-- ) {
! 1248: /* reduce mat[i] by mat[l] */
! 1249: if ( hc = t[ind[l]] ) {
! 1250: /* mat[i] = mat[i]-hc*mat[l] */
! 1251: for ( k = ind[l], hc = md-hc, s = mat[l]+k, tk = t+k;
! 1252: k < col; k++, tk++, s++ )
! 1253: if ( *s )
! 1254: DMAR(*s,hc,*tk,md,*tk)
! 1255: }
! 1256: }
! 1257: }
! 1258: #endif
! 1259: }
! 1260: /*
! 1261: mat[i] : reducers (i=0,...,nred-1)
! 1262: mat[0] < mat[1] < ... < mat[nred-1] w.r.t the term order
! 1263: */
! 1264:
! 1265: void reduce_sp_by_red_mod(sp,redmat,ind,nred,col,md)
! 1266: int *sp,**redmat;
! 1267: int *ind;
! 1268: int nred,col;
! 1269: int md;
! 1270: {
! 1271: int i,j,k,hc,zzz;
! 1272: int *t,*s,*tj;
! 1273:
! 1274: /* reduce the spolys by redmat */
! 1275: for ( i = nred-1; i >= 0; i-- ) {
! 1276: /* reduce sp by redmat[i] */
! 1277: if ( hc = sp[ind[i]] ) {
! 1278: /* sp = sp-hc*redmat[i] */
! 1279: j = ind[i];
! 1280: hc = md-hc;
! 1281: s = redmat[i]+j;
! 1282: tj = sp+j;
! 1283: for ( k = col-j; k >= 0; k-- ) {
! 1284: if ( zzz = *s ) { DMAR(zzz,hc,*tj,md,*tj) } tj++; s++;
! 1285: }
! 1286: }
! 1287: }
! 1288: }
! 1289:
! 1290: #define ONE_STEP2 if ( zzz = *pk ) { DMAR(zzz,a,*tk,md,*tk) } pk++; tk++;
! 1291:
1.1 noro 1292: int generic_gauss_elim_mod(mat,row,col,md,colstat)
1293: int **mat;
1294: int row,col,md;
1295: int *colstat;
1296: {
1.4 ! noro 1297: int i,j,k,l,inv,a,rank,zzz;
! 1298: int *t,*pivot,*pk,*tk;
1.1 noro 1299:
1300: for ( rank = 0, j = 0; j < col; j++ ) {
1301: for ( i = rank; i < row && !mat[i][j]; i++ );
1302: if ( i == row ) {
1303: colstat[j] = 0;
1304: continue;
1305: } else
1306: colstat[j] = 1;
1307: if ( i != rank ) {
1308: t = mat[i]; mat[i] = mat[rank]; mat[rank] = t;
1309: }
1310: pivot = mat[rank];
1311: inv = invm(pivot[j],md);
1.4 ! noro 1312: for ( k = j, pk = pivot+k; k < col; k++, pk++ )
! 1313: if ( *pk ) {
! 1314: DMAR(*pk,inv,0,md,*pk)
1.1 noro 1315: }
1316: for ( i = rank+1; i < row; i++ ) {
1317: t = mat[i];
1.4 ! noro 1318: if ( a = t[j] ) {
! 1319: a = md - a; pk = pivot+j; tk = t+j;
! 1320: k = col-j;
! 1321: for ( ; k >= 64; k -= 64 ) {
! 1322: ONE_STEP2 ONE_STEP2 ONE_STEP2 ONE_STEP2
! 1323: ONE_STEP2 ONE_STEP2 ONE_STEP2 ONE_STEP2
! 1324: ONE_STEP2 ONE_STEP2 ONE_STEP2 ONE_STEP2
! 1325: ONE_STEP2 ONE_STEP2 ONE_STEP2 ONE_STEP2
! 1326: ONE_STEP2 ONE_STEP2 ONE_STEP2 ONE_STEP2
! 1327: ONE_STEP2 ONE_STEP2 ONE_STEP2 ONE_STEP2
! 1328: ONE_STEP2 ONE_STEP2 ONE_STEP2 ONE_STEP2
! 1329: ONE_STEP2 ONE_STEP2 ONE_STEP2 ONE_STEP2
! 1330: ONE_STEP2 ONE_STEP2 ONE_STEP2 ONE_STEP2
! 1331: ONE_STEP2 ONE_STEP2 ONE_STEP2 ONE_STEP2
! 1332: ONE_STEP2 ONE_STEP2 ONE_STEP2 ONE_STEP2
! 1333: ONE_STEP2 ONE_STEP2 ONE_STEP2 ONE_STEP2
! 1334: ONE_STEP2 ONE_STEP2 ONE_STEP2 ONE_STEP2
! 1335: ONE_STEP2 ONE_STEP2 ONE_STEP2 ONE_STEP2
! 1336: ONE_STEP2 ONE_STEP2 ONE_STEP2 ONE_STEP2
! 1337: ONE_STEP2 ONE_STEP2 ONE_STEP2 ONE_STEP2
! 1338: }
! 1339: for ( ; k >= 0; k -- ) {
! 1340: if ( zzz = *pk ) { DMAR(zzz,a,*tk,md,*tk) } pk++; tk++;
! 1341: }
! 1342: }
1.1 noro 1343: }
1344: rank++;
1345: }
1346: for ( j = col-1, l = rank-1; j >= 0; j-- )
1347: if ( colstat[j] ) {
1348: pivot = mat[l];
1349: for ( i = 0; i < l; i++ ) {
1350: t = mat[i];
1.4 ! noro 1351: if ( a = t[j] ) {
! 1352: a = md-a; pk = pivot+j; tk = t+j;
! 1353: k = col-j;
! 1354: for ( ; k >= 64; k -= 64 ) {
! 1355: ONE_STEP2 ONE_STEP2 ONE_STEP2 ONE_STEP2
! 1356: ONE_STEP2 ONE_STEP2 ONE_STEP2 ONE_STEP2
! 1357: ONE_STEP2 ONE_STEP2 ONE_STEP2 ONE_STEP2
! 1358: ONE_STEP2 ONE_STEP2 ONE_STEP2 ONE_STEP2
! 1359: ONE_STEP2 ONE_STEP2 ONE_STEP2 ONE_STEP2
! 1360: ONE_STEP2 ONE_STEP2 ONE_STEP2 ONE_STEP2
! 1361: ONE_STEP2 ONE_STEP2 ONE_STEP2 ONE_STEP2
! 1362: ONE_STEP2 ONE_STEP2 ONE_STEP2 ONE_STEP2
! 1363: ONE_STEP2 ONE_STEP2 ONE_STEP2 ONE_STEP2
! 1364: ONE_STEP2 ONE_STEP2 ONE_STEP2 ONE_STEP2
! 1365: ONE_STEP2 ONE_STEP2 ONE_STEP2 ONE_STEP2
! 1366: ONE_STEP2 ONE_STEP2 ONE_STEP2 ONE_STEP2
! 1367: ONE_STEP2 ONE_STEP2 ONE_STEP2 ONE_STEP2
! 1368: ONE_STEP2 ONE_STEP2 ONE_STEP2 ONE_STEP2
! 1369: ONE_STEP2 ONE_STEP2 ONE_STEP2 ONE_STEP2
! 1370: ONE_STEP2 ONE_STEP2 ONE_STEP2 ONE_STEP2
! 1371: }
! 1372: for ( ; k >= 0; k -- ) {
! 1373: if ( zzz = *pk ) { DMAR(zzz,a,*tk,md,*tk) } pk++; tk++;
! 1374: }
! 1375: }
1.1 noro 1376: }
1377: l--;
1378: }
1379: return rank;
1380: }
1381:
1382: /* LU decomposition; a[i][i] = 1/U[i][i] */
1383:
1384: int lu_gfmmat(mat,md,perm)
1385: GFMMAT mat;
1386: unsigned int md;
1387: int *perm;
1388: {
1389: int row,col;
1390: int i,j,k,l;
1391: unsigned int *t,*pivot;
1392: unsigned int **a;
1393: unsigned int inv,m;
1394:
1395: row = mat->row; col = mat->col;
1396: a = mat->body;
1397: bzero(perm,row*sizeof(int));
1398:
1399: for ( i = 0; i < row; i++ )
1400: perm[i] = i;
1401: for ( k = 0; k < col; k++ ) {
1402: for ( i = k; i < row && !a[i][k]; i++ );
1403: if ( i == row )
1404: return 0;
1405: if ( i != k ) {
1406: j = perm[i]; perm[i] = perm[k]; perm[k] = j;
1407: t = a[i]; a[i] = a[k]; a[k] = t;
1408: }
1409: pivot = a[k];
1410: pivot[k] = inv = invm(pivot[k],md);
1411: for ( i = k+1; i < row; i++ ) {
1412: t = a[i];
1413: if ( m = t[k] ) {
1414: DMAR(inv,m,0,md,t[k])
1415: for ( j = k+1, m = md - t[k]; j < col; j++ )
1416: if ( pivot[j] ) {
1417: DMAR(m,pivot[j],t[j],md,t[j])
1418: }
1419: }
1420: }
1421: }
1422: return 1;
1423: }
1424:
1.3 noro 1425: /*
1426: Input
1427: a: a row x col matrix
1428: md : a modulus
1429:
1430: Output:
1431: return : d = the rank of mat
1432: a[0..(d-1)][0..(d-1)] : LU decomposition (a[i][i] = 1/U[i][i])
1433: rinfo: array of length row
1434: cinfo: array of length col
1435: i-th row in new a <-> rinfo[i]-th row in old a
1436: cinfo[j]=1 <=> j-th column is contained in the LU decomp.
1437: */
1438:
1439: int find_lhs_and_lu_mod(a,row,col,md,rinfo,cinfo)
1440: unsigned int **a;
1441: unsigned int md;
1442: int **rinfo,**cinfo;
1443: {
1444: int i,j,k,l,d;
1445: int *rp,*cp;
1446: unsigned int *t,*pivot;
1447: unsigned int inv,m;
1448:
1449: *rinfo = rp = (int *)MALLOC_ATOMIC(row*sizeof(int));
1450: *cinfo = cp = (int *)MALLOC_ATOMIC(col*sizeof(int));
1451: for ( i = 0; i < row; i++ )
1452: rp[i] = i;
1453: for ( k = 0, d = 0; k < col; k++ ) {
1454: for ( i = d; i < row && !a[i][k]; i++ );
1455: if ( i == row ) {
1456: cp[k] = 0;
1457: continue;
1458: } else
1459: cp[k] = 1;
1460: if ( i != d ) {
1461: j = rp[i]; rp[i] = rp[d]; rp[d] = j;
1462: t = a[i]; a[i] = a[d]; a[d] = t;
1463: }
1464: pivot = a[d];
1465: pivot[k] = inv = invm(pivot[k],md);
1466: for ( i = d+1; i < row; i++ ) {
1467: t = a[i];
1468: if ( m = t[k] ) {
1469: DMAR(inv,m,0,md,t[k])
1470: for ( j = k+1, m = md - t[k]; j < col; j++ )
1471: if ( pivot[j] ) {
1472: DMAR(m,pivot[j],t[j],md,t[j])
1473: }
1474: }
1475: }
1476: d++;
1477: }
1478: return d;
1479: }
1480:
1481: /*
1482: Input
1483: a : n x n matrix; a result of LU-decomposition
1484: md : modulus
1485: b : n x l matrix
1486: Output
1487: b = a^(-1)b
1488: */
1489:
1490: void solve_by_lu_mod(a,n,md,b,l)
1491: int **a;
1492: int n;
1493: int md;
1494: int **b;
1495: int l;
1496: {
1497: unsigned int *y,*c;
1498: int i,j,k;
1499: unsigned int t,m,m2;
1500:
1501: y = (int *)MALLOC_ATOMIC(n*sizeof(int));
1502: c = (int *)MALLOC_ATOMIC(n*sizeof(int));
1503: m2 = md>>1;
1504: for ( k = 0; k < l; k++ ) {
1505: /* copy b[.][k] to c */
1506: for ( i = 0; i < n; i++ )
1507: c[i] = (unsigned int)b[i][k];
1508: /* solve Ly=c */
1509: for ( i = 0; i < n; i++ ) {
1510: for ( t = c[i], j = 0; j < i; j++ )
1511: if ( a[i][j] ) {
1512: m = md - a[i][j];
1513: DMAR(m,y[j],t,md,t)
1514: }
1515: y[i] = t;
1516: }
1517: /* solve Uc=y */
1518: for ( i = n-1; i >= 0; i-- ) {
1519: for ( t = y[i], j =i+1; j < n; j++ )
1520: if ( a[i][j] ) {
1521: m = md - a[i][j];
1522: DMAR(m,c[j],t,md,t)
1523: }
1524: /* a[i][i] = 1/U[i][i] */
1525: DMAR(t,a[i][i],0,md,c[i])
1526: }
1527: /* copy c to b[.][k] with normalization */
1528: for ( i = 0; i < n; i++ )
1529: b[i][k] = (int)(c[i]>m2 ? c[i]-md : c[i]);
1530: }
1531: }
1532:
1.1 noro 1533: void Pleqm1(arg,rp)
1534: NODE arg;
1535: VECT *rp;
1536: {
1537: MAT m;
1538: VECT vect;
1539: pointer **mat;
1540: Q *v;
1541: Q q;
1542: int **wmat;
1543: int md,i,j,row,col,t,n,status;
1544:
1545: asir_assert(ARG0(arg),O_MAT,"leqm1");
1546: asir_assert(ARG1(arg),O_N,"leqm1");
1547: m = (MAT)ARG0(arg); md = QTOS((Q)ARG1(arg));
1548: row = m->row; col = m->col; mat = m->body;
1549: wmat = (int **)almat(row,col);
1550: for ( i = 0; i < row; i++ )
1551: for ( j = 0; j < col; j++ )
1552: if ( q = (Q)mat[i][j] ) {
1553: t = rem(NM(q),md);
1554: if ( SGN(q) < 0 )
1555: t = (md - t) % md;
1556: wmat[i][j] = t;
1557: } else
1558: wmat[i][j] = 0;
1559: status = gauss_elim_mod1(wmat,row,col,md);
1560: if ( status < 0 )
1561: *rp = 0;
1562: else if ( status > 0 )
1563: *rp = (VECT)ONE;
1564: else {
1565: n = col - 1;
1566: MKVECT(vect,n);
1567: for ( i = 0, v = (Q *)vect->body; i < n; i++ ) {
1568: t = (md-wmat[i][n])%md; STOQ(t,v[i]);
1569: }
1570: *rp = vect;
1571: }
1572: }
1573:
1574: gauss_elim_mod1(mat,row,col,md)
1575: int **mat;
1576: int row,col,md;
1577: {
1578: int i,j,k,inv,a,n;
1579: int *t,*pivot;
1580:
1581: n = col - 1;
1582: for ( j = 0; j < n; j++ ) {
1583: for ( i = j; i < row && !mat[i][j]; i++ );
1584: if ( i == row )
1585: return 1;
1586: if ( i != j ) {
1587: t = mat[i]; mat[i] = mat[j]; mat[j] = t;
1588: }
1589: pivot = mat[j];
1590: inv = invm(pivot[j],md);
1591: for ( k = j; k <= n; k++ )
1592: pivot[k] = dmar(pivot[k],inv,0,md);
1593: for ( i = j+1; i < row; i++ ) {
1594: t = mat[i];
1595: if ( i != j && (a = t[j]) )
1596: for ( k = j, a = md - a; k <= n; k++ )
1597: t[k] = dmar(pivot[k],a,t[k],md);
1598: }
1599: }
1600: for ( i = n; i < row && !mat[i][n]; i++ );
1601: if ( i == row ) {
1602: for ( j = n-1; j >= 0; j-- ) {
1603: for ( i = j-1, a = (md-mat[j][n])%md; i >= 0; i-- ) {
1604: mat[i][n] = dmar(mat[i][j],a,mat[i][n],md);
1605: mat[i][j] = 0;
1606: }
1607: }
1608: return 0;
1609: } else
1610: return -1;
1611: }
1612:
1613: void Pgeninvm(arg,rp)
1614: NODE arg;
1615: LIST *rp;
1616: {
1617: MAT m;
1618: pointer **mat;
1619: Q **tmat;
1620: Q q;
1621: unsigned int **wmat;
1622: int md,i,j,row,col,t,status;
1623: MAT mat1,mat2;
1624: NODE node1,node2;
1625:
1626: asir_assert(ARG0(arg),O_MAT,"leqm1");
1627: asir_assert(ARG1(arg),O_N,"leqm1");
1628: m = (MAT)ARG0(arg); md = QTOS((Q)ARG1(arg));
1629: row = m->row; col = m->col; mat = m->body;
1630: wmat = (unsigned int **)almat(row,col+row);
1631: for ( i = 0; i < row; i++ ) {
1632: bzero((char *)wmat[i],(col+row)*sizeof(int));
1633: for ( j = 0; j < col; j++ )
1634: if ( q = (Q)mat[i][j] ) {
1635: t = rem(NM(q),md);
1636: if ( SGN(q) < 0 )
1637: t = (md - t) % md;
1638: wmat[i][j] = t;
1639: }
1640: wmat[i][col+i] = 1;
1641: }
1642: status = gauss_elim_geninv_mod(wmat,row,col,md);
1643: if ( status > 0 )
1644: *rp = 0;
1645: else {
1646: MKMAT(mat1,col,row); MKMAT(mat2,row-col,row);
1647: for ( i = 0, tmat = (Q **)mat1->body; i < col; i++ )
1648: for ( j = 0; j < row; j++ )
1649: STOQ(wmat[i][j+col],tmat[i][j]);
1650: for ( tmat = (Q **)mat2->body; i < row; i++ )
1651: for ( j = 0; j < row; j++ )
1652: STOQ(wmat[i][j+col],tmat[i-col][j]);
1653: MKNODE(node2,mat2,0); MKNODE(node1,mat1,node2); MKLIST(*rp,node1);
1654: }
1655: }
1656:
1657: int gauss_elim_geninv_mod(mat,row,col,md)
1658: unsigned int **mat;
1659: int row,col,md;
1660: {
1661: int i,j,k,inv,a,n,m;
1662: unsigned int *t,*pivot;
1663:
1664: n = col; m = row+col;
1665: for ( j = 0; j < n; j++ ) {
1666: for ( i = j; i < row && !mat[i][j]; i++ );
1667: if ( i == row )
1668: return 1;
1669: if ( i != j ) {
1670: t = mat[i]; mat[i] = mat[j]; mat[j] = t;
1671: }
1672: pivot = mat[j];
1673: inv = invm(pivot[j],md);
1674: for ( k = j; k < m; k++ )
1675: pivot[k] = dmar(pivot[k],inv,0,md);
1676: for ( i = j+1; i < row; i++ ) {
1677: t = mat[i];
1678: if ( a = t[j] )
1679: for ( k = j, a = md - a; k < m; k++ )
1680: t[k] = dmar(pivot[k],a,t[k],md);
1681: }
1682: }
1683: for ( j = n-1; j >= 0; j-- ) {
1684: pivot = mat[j];
1685: for ( i = j-1; i >= 0; i-- ) {
1686: t = mat[i];
1687: if ( a = t[j] )
1688: for ( k = j, a = md - a; k < m; k++ )
1689: t[k] = dmar(pivot[k],a,t[k],md);
1690: }
1691: }
1692: return 0;
1693: }
1694:
1695: void Psolve_by_lu_gfmmat(arg,rp)
1696: NODE arg;
1697: VECT *rp;
1698: {
1699: GFMMAT lu;
1700: Q *perm,*rhs,*v;
1701: int n,i;
1702: unsigned int md;
1703: unsigned int *b,*sol;
1704: VECT r;
1705:
1706: lu = (GFMMAT)ARG0(arg);
1707: perm = (Q *)BDY((VECT)ARG1(arg));
1708: rhs = (Q *)BDY((VECT)ARG2(arg));
1709: md = (unsigned int)QTOS((Q)ARG3(arg));
1710: n = lu->col;
1711: b = (unsigned int *)MALLOC_ATOMIC(n*sizeof(int));
1712: sol = (unsigned int *)MALLOC_ATOMIC(n*sizeof(int));
1713: for ( i = 0; i < n; i++ )
1714: b[i] = QTOS(rhs[QTOS(perm[i])]);
1715: solve_by_lu_gfmmat(lu,md,b,sol);
1716: MKVECT(r,n);
1717: for ( i = 0, v = (Q *)r->body; i < n; i++ )
1718: STOQ(sol[i],v[i]);
1719: *rp = r;
1720: }
1721:
1722: void solve_by_lu_gfmmat(lu,md,b,x)
1723: GFMMAT lu;
1724: unsigned int md;
1725: unsigned int *b;
1726: unsigned int *x;
1727: {
1728: int n;
1729: unsigned int **a;
1730: unsigned int *y;
1731: int i,j;
1732: unsigned int t,m;
1733:
1734: n = lu->col;
1735: a = lu->body;
1736: y = (unsigned int *)MALLOC_ATOMIC(n*sizeof(int));
1737: /* solve Ly=b */
1738: for ( i = 0; i < n; i++ ) {
1739: for ( t = b[i], j = 0; j < i; j++ )
1740: if ( a[i][j] ) {
1741: m = md - a[i][j];
1742: DMAR(m,y[j],t,md,t)
1743: }
1744: y[i] = t;
1745: }
1746: /* solve Ux=y */
1747: for ( i = n-1; i >= 0; i-- ) {
1748: for ( t = y[i], j =i+1; j < n; j++ )
1749: if ( a[i][j] ) {
1750: m = md - a[i][j];
1751: DMAR(m,x[j],t,md,t)
1752: }
1753: /* a[i][i] = 1/U[i][i] */
1754: DMAR(t,a[i][i],0,md,x[i])
1755: }
1756: }
1757:
1758: void Plu_gfmmat(arg,rp)
1759: NODE arg;
1760: LIST *rp;
1761: {
1762: MAT m;
1763: GFMMAT mm;
1764: unsigned int md;
1765: int i,row,col,status;
1766: int *iperm;
1767: Q *v;
1768: VECT perm;
1769: NODE n0;
1770:
1771: asir_assert(ARG0(arg),O_MAT,"mat_to_gfmmat");
1772: asir_assert(ARG1(arg),O_N,"mat_to_gfmmat");
1773: m = (MAT)ARG0(arg); md = (unsigned int)QTOS((Q)ARG1(arg));
1774: mat_to_gfmmat(m,md,&mm);
1775: row = m->row;
1776: col = m->col;
1777: iperm = (int *)MALLOC_ATOMIC(row*sizeof(int));
1778: status = lu_gfmmat(mm,md,iperm);
1779: if ( !status )
1780: n0 = 0;
1781: else {
1782: MKVECT(perm,row);
1783: for ( i = 0, v = (Q *)perm->body; i < row; i++ )
1784: STOQ(iperm[i],v[i]);
1785: n0 = mknode(2,mm,perm);
1786: }
1787: MKLIST(*rp,n0);
1788: }
1789:
1790: void Pmat_to_gfmmat(arg,rp)
1791: NODE arg;
1792: GFMMAT *rp;
1793: {
1794: MAT m;
1795: unsigned int md;
1796:
1797: asir_assert(ARG0(arg),O_MAT,"mat_to_gfmmat");
1798: asir_assert(ARG1(arg),O_N,"mat_to_gfmmat");
1799: m = (MAT)ARG0(arg); md = (unsigned int)QTOS((Q)ARG1(arg));
1800: mat_to_gfmmat(m,md,rp);
1801: }
1802:
1803: void mat_to_gfmmat(m,md,rp)
1804: MAT m;
1805: unsigned int md;
1806: GFMMAT *rp;
1807: {
1808: unsigned int **wmat;
1809: unsigned int t;
1810: Q **mat;
1811: Q q;
1812: int i,j,row,col;
1813:
1814: row = m->row; col = m->col; mat = (Q **)m->body;
1815: wmat = (unsigned int **)almat(row,col);
1816: for ( i = 0; i < row; i++ ) {
1817: bzero((char *)wmat[i],col*sizeof(unsigned int));
1818: for ( j = 0; j < col; j++ )
1819: if ( q = mat[i][j] ) {
1820: t = (unsigned int)rem(NM(q),md);
1821: if ( SGN(q) < 0 )
1822: t = (md - t) % md;
1823: wmat[i][j] = t;
1824: }
1825: }
1826: TOGFMMAT(row,col,wmat,*rp);
1827: }
1828:
1829: void Pgeninvm_swap(arg,rp)
1830: NODE arg;
1831: LIST *rp;
1832: {
1833: MAT m;
1834: pointer **mat;
1835: Q **tmat;
1836: Q *tvect;
1837: Q q;
1838: unsigned int **wmat,**invmat;
1839: int *index;
1840: unsigned int t,md;
1841: int i,j,row,col,status;
1842: MAT mat1;
1843: VECT vect1;
1844: NODE node1,node2;
1845:
1846: asir_assert(ARG0(arg),O_MAT,"geninvm_swap");
1847: asir_assert(ARG1(arg),O_N,"geninvm_swap");
1848: m = (MAT)ARG0(arg); md = QTOS((Q)ARG1(arg));
1849: row = m->row; col = m->col; mat = m->body;
1850: wmat = (unsigned int **)almat(row,col+row);
1851: for ( i = 0; i < row; i++ ) {
1852: bzero((char *)wmat[i],(col+row)*sizeof(int));
1853: for ( j = 0; j < col; j++ )
1854: if ( q = (Q)mat[i][j] ) {
1855: t = (unsigned int)rem(NM(q),md);
1856: if ( SGN(q) < 0 )
1857: t = (md - t) % md;
1858: wmat[i][j] = t;
1859: }
1860: wmat[i][col+i] = 1;
1861: }
1862: status = gauss_elim_geninv_mod_swap(wmat,row,col,md,&invmat,&index);
1863: if ( status > 0 )
1864: *rp = 0;
1865: else {
1866: MKMAT(mat1,col,col);
1867: for ( i = 0, tmat = (Q **)mat1->body; i < col; i++ )
1868: for ( j = 0; j < col; j++ )
1869: UTOQ(invmat[i][j],tmat[i][j]);
1870: MKVECT(vect1,row);
1871: for ( i = 0, tvect = (Q *)vect1->body; i < row; i++ )
1872: STOQ(index[i],tvect[i]);
1873: MKNODE(node2,vect1,0); MKNODE(node1,mat1,node2); MKLIST(*rp,node1);
1874: }
1875: }
1876:
1877: gauss_elim_geninv_mod_swap(mat,row,col,md,invmatp,indexp)
1878: unsigned int **mat;
1879: int row,col;
1880: unsigned int md;
1881: unsigned int ***invmatp;
1882: int **indexp;
1883: {
1884: int i,j,k,inv,a,n,m;
1885: unsigned int *t,*pivot,*s;
1886: int *index;
1887: unsigned int **invmat;
1888:
1889: n = col; m = row+col;
1890: *indexp = index = (int *)MALLOC_ATOMIC(row*sizeof(int));
1891: for ( i = 0; i < row; i++ )
1892: index[i] = i;
1893: for ( j = 0; j < n; j++ ) {
1894: for ( i = j; i < row && !mat[i][j]; i++ );
1895: if ( i == row ) {
1896: *indexp = 0; *invmatp = 0; return 1;
1897: }
1898: if ( i != j ) {
1899: t = mat[i]; mat[i] = mat[j]; mat[j] = t;
1900: k = index[i]; index[i] = index[j]; index[j] = k;
1901: }
1902: pivot = mat[j];
1903: inv = (unsigned int)invm(pivot[j],md);
1904: for ( k = j; k < m; k++ )
1905: if ( pivot[k] )
1906: pivot[k] = (unsigned int)dmar(pivot[k],inv,0,md);
1907: for ( i = j+1; i < row; i++ ) {
1908: t = mat[i];
1909: if ( a = t[j] )
1910: for ( k = j, a = md - a; k < m; k++ )
1911: if ( pivot[k] )
1912: t[k] = dmar(pivot[k],a,t[k],md);
1913: }
1914: }
1915: for ( j = n-1; j >= 0; j-- ) {
1916: pivot = mat[j];
1917: for ( i = j-1; i >= 0; i-- ) {
1918: t = mat[i];
1919: if ( a = t[j] )
1920: for ( k = j, a = md - a; k < m; k++ )
1921: if ( pivot[k] )
1922: t[k] = dmar(pivot[k],a,t[k],md);
1923: }
1924: }
1925: *invmatp = invmat = (unsigned int **)almat(col,col);
1926: for ( i = 0; i < col; i++ )
1927: for ( j = 0, s = invmat[i], t = mat[i]; j < col; j++ )
1928: s[j] = t[col+index[j]];
1929: return 0;
1930: }
1931:
1932: void _addn(N,N,N);
1933: int _subn(N,N,N);
1934: void _muln(N,N,N);
1935:
1936: void inner_product_int(a,b,n,r)
1937: Q *a,*b;
1938: int n;
1939: Q *r;
1940: {
1941: int la,lb,i;
1942: int sgn,sgn1;
1943: N wm,wma,sum,t;
1944:
1945: for ( la = lb = 0, i = 0; i < n; i++ ) {
1946: if ( a[i] )
1947: if ( DN(a[i]) )
1948: error("inner_product_int : invalid argument");
1949: else
1950: la = MAX(PL(NM(a[i])),la);
1951: if ( b[i] )
1952: if ( DN(b[i]) )
1953: error("inner_product_int : invalid argument");
1954: else
1955: lb = MAX(PL(NM(b[i])),lb);
1956: }
1957: sgn = 0;
1958: sum= NALLOC(la+lb+2);
1959: bzero((char *)sum,(la+lb+3)*sizeof(unsigned int));
1960: wm = NALLOC(la+lb+2);
1961: wma = NALLOC(la+lb+2);
1962: for ( i = 0; i < n; i++ ) {
1963: if ( !a[i] || !b[i] )
1964: continue;
1965: _muln(NM(a[i]),NM(b[i]),wm);
1966: sgn1 = SGN(a[i])*SGN(b[i]);
1967: if ( !sgn ) {
1968: sgn = sgn1;
1969: t = wm; wm = sum; sum = t;
1970: } else if ( sgn == sgn1 ) {
1971: _addn(sum,wm,wma);
1972: if ( !PL(wma) )
1973: sgn = 0;
1974: t = wma; wma = sum; sum = t;
1975: } else {
1976: /* sgn*sum+sgn1*wm = sgn*(sum-wm) */
1977: sgn *= _subn(sum,wm,wma);
1978: t = wma; wma = sum; sum = t;
1979: }
1980: }
1981: GC_free(wm);
1982: GC_free(wma);
1983: if ( !sgn ) {
1984: GC_free(sum);
1985: *r = 0;
1986: } else
1987: NTOQ(sum,sgn,*r);
1988: }
1989:
1.3 noro 1990: /* (k,l) element of a*b where a: .x n matrix, b: n x . integer matrix */
1991:
1992: void inner_product_mat_int_mod(a,b,n,k,l,r)
1993: Q **a;
1994: int **b;
1995: int n,k,l;
1996: Q *r;
1997: {
1998: int la,lb,i;
1999: int sgn,sgn1;
2000: N wm,wma,sum,t;
2001: Q aki;
2002: int bil,bilsgn;
2003: struct oN tn;
2004:
2005: for ( la = 0, i = 0; i < n; i++ ) {
2006: if ( aki = a[k][i] )
2007: if ( DN(aki) )
2008: error("inner_product_int : invalid argument");
2009: else
2010: la = MAX(PL(NM(aki)),la);
2011: }
2012: lb = 1;
2013: sgn = 0;
2014: sum= NALLOC(la+lb+2);
2015: bzero((char *)sum,(la+lb+3)*sizeof(unsigned int));
2016: wm = NALLOC(la+lb+2);
2017: wma = NALLOC(la+lb+2);
2018: for ( i = 0; i < n; i++ ) {
2019: if ( !(aki = a[k][i]) || !(bil = b[i][l]) )
2020: continue;
2021: tn.p = 1;
2022: if ( bil > 0 ) {
2023: tn.b[0] = bil; bilsgn = 1;
2024: } else {
2025: tn.b[0] = -bil; bilsgn = -1;
2026: }
2027: _muln(NM(aki),&tn,wm);
2028: sgn1 = SGN(aki)*bilsgn;
2029: if ( !sgn ) {
2030: sgn = sgn1;
2031: t = wm; wm = sum; sum = t;
2032: } else if ( sgn == sgn1 ) {
2033: _addn(sum,wm,wma);
2034: if ( !PL(wma) )
2035: sgn = 0;
2036: t = wma; wma = sum; sum = t;
2037: } else {
2038: /* sgn*sum+sgn1*wm = sgn*(sum-wm) */
2039: sgn *= _subn(sum,wm,wma);
2040: t = wma; wma = sum; sum = t;
2041: }
2042: }
2043: GC_free(wm);
2044: GC_free(wma);
2045: if ( !sgn ) {
2046: GC_free(sum);
2047: *r = 0;
2048: } else
2049: NTOQ(sum,sgn,*r);
2050: }
2051:
1.1 noro 2052: void Pmul_mat_vect_int(arg,rp)
2053: NODE arg;
2054: VECT *rp;
2055: {
2056: MAT mat;
2057: VECT vect,r;
2058: int row,col,i;
2059:
2060: mat = (MAT)ARG0(arg);
2061: vect = (VECT)ARG1(arg);
2062: row = mat->row;
2063: col = mat->col;
2064: MKVECT(r,row);
2065: for ( i = 0; i < row; i++ )
2066: inner_product_int(mat->body[i],vect->body,col,&r->body[i]);
2067: *rp = r;
2068: }
2069:
2070: void Pnbpoly_up2(arg,rp)
2071: NODE arg;
2072: GF2N *rp;
2073: {
2074: int m,type,ret;
2075: UP2 r;
2076:
2077: m = QTOS((Q)ARG0(arg));
2078: type = QTOS((Q)ARG1(arg));
2079: ret = generate_ONB_polynomial(&r,m,type);
2080: if ( ret == 0 )
2081: MKGF2N(r,*rp);
2082: else
2083: *rp = 0;
2084: }
2085:
2086: void Px962_irredpoly_up2(arg,rp)
2087: NODE arg;
2088: GF2N *rp;
2089: {
2090: int m,type,ret,w;
2091: GF2N prev;
2092: UP2 r;
2093:
2094: m = QTOS((Q)ARG0(arg));
2095: prev = (GF2N)ARG1(arg);
2096: if ( !prev ) {
2097: w = (m>>5)+1; NEWUP2(r,w); r->w = 0;
2098: bzero((char *)r->b,w*sizeof(unsigned int));
2099: } else {
2100: r = prev->body;
2101: if ( degup2(r) != m ) {
2102: w = (m>>5)+1; NEWUP2(r,w); r->w = 0;
2103: bzero((char *)r->b,w*sizeof(unsigned int));
2104: }
2105: }
2106: ret = _generate_irreducible_polynomial(r,m,type);
2107: if ( ret == 0 )
2108: MKGF2N(r,*rp);
2109: else
2110: *rp = 0;
2111: }
2112:
2113: void Pirredpoly_up2(arg,rp)
2114: NODE arg;
2115: GF2N *rp;
2116: {
2117: int m,type,ret,w;
2118: GF2N prev;
2119: UP2 r;
2120:
2121: m = QTOS((Q)ARG0(arg));
2122: prev = (GF2N)ARG1(arg);
2123: if ( !prev ) {
2124: w = (m>>5)+1; NEWUP2(r,w); r->w = 0;
2125: bzero((char *)r->b,w*sizeof(unsigned int));
2126: } else {
2127: r = prev->body;
2128: if ( degup2(r) != m ) {
2129: w = (m>>5)+1; NEWUP2(r,w); r->w = 0;
2130: bzero((char *)r->b,w*sizeof(unsigned int));
2131: }
2132: }
2133: ret = _generate_good_irreducible_polynomial(r,m,type);
2134: if ( ret == 0 )
2135: MKGF2N(r,*rp);
2136: else
2137: *rp = 0;
2138: }
2139:
2140: /*
2141: * f = type 'type' normal polynomial of degree m if exists
2142: * IEEE P1363 A.7.2
2143: *
2144: * return value : 0 --- exists
2145: * 1 --- does not exist
2146: * -1 --- failure (memory allocation error)
2147: */
2148:
2149: int generate_ONB_polynomial(UP2 *rp,int m,int type)
2150: {
2151: int i,r;
2152: int w;
2153: UP2 f,f0,f1,f2,t;
2154:
2155: w = (m>>5)+1;
2156: switch ( type ) {
2157: case 1:
2158: if ( !TypeT_NB_check(m,1) ) return 1;
2159: NEWUP2(f,w); *rp = f; f->w = w;
2160: /* set all the bits */
2161: for ( i = 0; i < w; i++ )
2162: f->b[i] = 0xffffffff;
2163: /* mask the top word if necessary */
2164: if ( r = (m+1)&31 )
2165: f->b[w-1] &= (1<<r)-1;
2166: return 0;
2167: break;
2168: case 2:
2169: if ( !TypeT_NB_check(m,2) ) return 1;
2170: NEWUP2(f,w); *rp = f;
2171: W_NEWUP2(f0,w);
2172: W_NEWUP2(f1,w);
2173: W_NEWUP2(f2,w);
2174:
2175: /* recursion for genrating Type II normal polynomial */
2176:
2177: /* f0 = 1, f1 = t+1 */
2178: f0->w = 1; f0->b[0] = 1;
2179: f1->w = 1; f1->b[0] = 3;
2180: for ( i = 2; i <= m; i++ ) {
2181: /* f2 = t*f1+f0 */
2182: _bshiftup2(f1,-1,f2);
2183: _addup2_destructive(f2,f0);
2184: /* cyclic change of the variables */
2185: t = f0; f0 = f1; f1 = f2; f2 = t;
2186: }
2187: _copyup2(f1,f);
2188: return 0;
2189: break;
2190: default:
2191: return -1;
2192: break;
2193: }
2194: }
2195:
2196: /*
2197: * f = an irreducible trinomial or pentanomial of degree d 'after' f
2198: * return value : 0 --- exists
2199: * 1 --- does not exist (exhaustion)
2200: */
2201:
2202: int _generate_irreducible_polynomial(UP2 f,int d)
2203: {
2204: int ret,i,j,k,nz,i0,j0,k0;
2205: int w;
2206: unsigned int *fd;
2207:
2208: /*
2209: * if f = x^d+x^i+1 then i0 <- i, j0 <- 0, k0 <-0.
2210: * if f = x^d+x^k+x^j+x^i+1 (k>j>i) then i0 <- i, j0 <- j, k0 <-k.
2211: * otherwise i0,j0,k0 is set to 0.
2212: */
2213:
2214: fd = f->b;
2215: w = (d>>5)+1;
2216: if ( f->w && (d==degup2(f)) ) {
2217: for ( nz = 0, i = d; i >= 0; i-- )
2218: if ( fd[i>>5]&(1<<(i&31)) ) nz++;
2219: switch ( nz ) {
2220: case 3:
2221: for ( i0 = 1; !(fd[i0>>5]&(1<<(i0&31))) ; i0++ );
2222: /* reset i0-th bit */
2223: fd[i0>>5] &= ~(1<<(i0&31));
2224: j0 = k0 = 0;
2225: break;
2226: case 5:
2227: for ( i0 = 1; !(fd[i0>>5]&(1<<(i0&31))) ; i0++ );
2228: /* reset i0-th bit */
2229: fd[i0>>5] &= ~(1<<(i0&31));
2230: for ( j0 = i0+1; !(fd[j0>>5]&(1<<(j0&31))) ; j0++ );
2231: /* reset j0-th bit */
2232: fd[j0>>5] &= ~(1<<(j0&31));
2233: for ( k0 = j0+1; !(fd[k0>>5]&(1<<(k0&31))) ; k0++ );
2234: /* reset k0-th bit */
2235: fd[k0>>5] &= ~(1<<(k0&31));
2236: break;
2237: default:
2238: f->w = 0; break;
2239: }
2240: } else
2241: f->w = 0;
2242:
2243: if ( !f->w ) {
2244: fd = f->b;
2245: f->w = w; fd[0] |= 1; fd[d>>5] |= (1<<(d&31));
2246: i0 = j0 = k0 = 0;
2247: }
2248: /* if j0 > 0 then f is already a pentanomial */
2249: if ( j0 > 0 ) goto PENTA;
2250:
2251: /* searching for an irreducible trinomial */
2252:
2253: for ( i = 1; 2*i <= d; i++ ) {
2254: /* skip the polynomials 'before' f */
2255: if ( i < i0 ) continue;
2256: if ( i == i0 ) { i0 = 0; continue; }
2257: /* set i-th bit */
2258: fd[i>>5] |= (1<<(i&31));
2259: ret = irredcheck_dddup2(f);
2260: if ( ret == 1 ) return 0;
2261: /* reset i-th bit */
2262: fd[i>>5] &= ~(1<<(i&31));
2263: }
2264:
2265: /* searching for an irreducible pentanomial */
2266: PENTA:
2267: for ( i = 1; i < d; i++ ) {
2268: /* skip the polynomials 'before' f */
2269: if ( i < i0 ) continue;
2270: if ( i == i0 ) i0 = 0;
2271: /* set i-th bit */
2272: fd[i>>5] |= (1<<(i&31));
2273: for ( j = i+1; j < d; j++ ) {
2274: /* skip the polynomials 'before' f */
2275: if ( j < j0 ) continue;
2276: if ( j == j0 ) j0 = 0;
2277: /* set j-th bit */
2278: fd[j>>5] |= (1<<(j&31));
2279: for ( k = j+1; k < d; k++ ) {
2280: /* skip the polynomials 'before' f */
2281: if ( k < k0 ) continue;
2282: else if ( k == k0 ) { k0 = 0; continue; }
2283: /* set k-th bit */
2284: fd[k>>5] |= (1<<(k&31));
2285: ret = irredcheck_dddup2(f);
2286: if ( ret == 1 ) return 0;
2287: /* reset k-th bit */
2288: fd[k>>5] &= ~(1<<(k&31));
2289: }
2290: /* reset j-th bit */
2291: fd[j>>5] &= ~(1<<(j&31));
2292: }
2293: /* reset i-th bit */
2294: fd[i>>5] &= ~(1<<(i&31));
2295: }
2296: /* exhausted */
2297: return 1;
2298: }
2299:
2300: /*
2301: * f = an irreducible trinomial or pentanomial of degree d 'after' f
2302: *
2303: * searching strategy:
2304: * trinomial x^d+x^i+1:
2305: * i is as small as possible.
2306: * trinomial x^d+x^i+x^j+x^k+1:
2307: * i is as small as possible.
2308: * For such i, j is as small as possible.
2309: * For such i and j, 'k' is as small as possible.
2310: *
2311: * return value : 0 --- exists
2312: * 1 --- does not exist (exhaustion)
2313: */
2314:
2315: int _generate_good_irreducible_polynomial(UP2 f,int d)
2316: {
2317: int ret,i,j,k,nz,i0,j0,k0;
2318: int w;
2319: unsigned int *fd;
2320:
2321: /*
2322: * if f = x^d+x^i+1 then i0 <- i, j0 <- 0, k0 <-0.
2323: * if f = x^d+x^k+x^j+x^i+1 (k>j>i) then i0 <- i, j0 <- j, k0 <-k.
2324: * otherwise i0,j0,k0 is set to 0.
2325: */
2326:
2327: fd = f->b;
2328: w = (d>>5)+1;
2329: if ( f->w && (d==degup2(f)) ) {
2330: for ( nz = 0, i = d; i >= 0; i-- )
2331: if ( fd[i>>5]&(1<<(i&31)) ) nz++;
2332: switch ( nz ) {
2333: case 3:
2334: for ( i0 = 1; !(fd[i0>>5]&(1<<(i0&31))) ; i0++ );
2335: /* reset i0-th bit */
2336: fd[i0>>5] &= ~(1<<(i0&31));
2337: j0 = k0 = 0;
2338: break;
2339: case 5:
2340: for ( i0 = 1; !(fd[i0>>5]&(1<<(i0&31))) ; i0++ );
2341: /* reset i0-th bit */
2342: fd[i0>>5] &= ~(1<<(i0&31));
2343: for ( j0 = i0+1; !(fd[j0>>5]&(1<<(j0&31))) ; j0++ );
2344: /* reset j0-th bit */
2345: fd[j0>>5] &= ~(1<<(j0&31));
2346: for ( k0 = j0+1; !(fd[k0>>5]&(1<<(k0&31))) ; k0++ );
2347: /* reset k0-th bit */
2348: fd[k0>>5] &= ~(1<<(k0&31));
2349: break;
2350: default:
2351: f->w = 0; break;
2352: }
2353: } else
2354: f->w = 0;
2355:
2356: if ( !f->w ) {
2357: fd = f->b;
2358: f->w = w; fd[0] |= 1; fd[d>>5] |= (1<<(d&31));
2359: i0 = j0 = k0 = 0;
2360: }
2361: /* if j0 > 0 then f is already a pentanomial */
2362: if ( j0 > 0 ) goto PENTA;
2363:
2364: /* searching for an irreducible trinomial */
2365:
2366: for ( i = 1; 2*i <= d; i++ ) {
2367: /* skip the polynomials 'before' f */
2368: if ( i < i0 ) continue;
2369: if ( i == i0 ) { i0 = 0; continue; }
2370: /* set i-th bit */
2371: fd[i>>5] |= (1<<(i&31));
2372: ret = irredcheck_dddup2(f);
2373: if ( ret == 1 ) return 0;
2374: /* reset i-th bit */
2375: fd[i>>5] &= ~(1<<(i&31));
2376: }
2377:
2378: /* searching for an irreducible pentanomial */
2379: PENTA:
2380: for ( i = 3; i < d; i++ ) {
2381: /* skip the polynomials 'before' f */
2382: if ( i < i0 ) continue;
2383: if ( i == i0 ) i0 = 0;
2384: /* set i-th bit */
2385: fd[i>>5] |= (1<<(i&31));
2386: for ( j = 2; j < i; j++ ) {
2387: /* skip the polynomials 'before' f */
2388: if ( j < j0 ) continue;
2389: if ( j == j0 ) j0 = 0;
2390: /* set j-th bit */
2391: fd[j>>5] |= (1<<(j&31));
2392: for ( k = 1; k < j; k++ ) {
2393: /* skip the polynomials 'before' f */
2394: if ( k < k0 ) continue;
2395: else if ( k == k0 ) { k0 = 0; continue; }
2396: /* set k-th bit */
2397: fd[k>>5] |= (1<<(k&31));
2398: ret = irredcheck_dddup2(f);
2399: if ( ret == 1 ) return 0;
2400: /* reset k-th bit */
2401: fd[k>>5] &= ~(1<<(k&31));
2402: }
2403: /* reset j-th bit */
2404: fd[j>>5] &= ~(1<<(j&31));
2405: }
2406: /* reset i-th bit */
2407: fd[i>>5] &= ~(1<<(i&31));
2408: }
2409: /* exhausted */
2410: return 1;
1.3 noro 2411: }
2412:
2413: printqmat(mat,row,col)
2414: Q **mat;
2415: int row,col;
2416: {
2417: int i,j;
2418:
2419: for ( i = 0; i < row; i++ ) {
2420: for ( j = 0; j < col; j++ ) {
2421: printnum(mat[i][j]); printf(" ");
2422: }
2423: printf("\n");
2424: }
2425: }
2426:
2427: printimat(mat,row,col)
2428: int **mat;
2429: int row,col;
2430: {
2431: int i,j;
2432:
2433: for ( i = 0; i < row; i++ ) {
2434: for ( j = 0; j < col; j++ ) {
2435: printf("%d ",mat[i][j]);
2436: }
2437: printf("\n");
2438: }
1.1 noro 2439: }