![[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.1 noro 1: /* $OpenXM: OpenXM/src/asir99/builtin/array.c,v 1.2 1999/11/23 07:14:14 noro Exp $ */
2: #include "ca.h"
3: #include "base.h"
4: #include "parse.h"
5: #include "inline.h"
6: /*
7: #undef DMAR
8: #define DMAR(a1,a2,a3,d,r) (r)=dmar(a1,a2,a3,d);
9: */
10:
11: extern int Print; /* XXX */
12:
13: void solve_by_lu_gfmmat(GFMMAT,unsigned int,unsigned int *,unsigned int *);
14: int lu_gfmmat(GFMMAT,unsigned int,int *);
15: void mat_to_gfmmat(MAT,unsigned int,GFMMAT *);
16:
17: int generic_gauss_elim_mod(int **,int,int,int,int *);
18: int generic_gauss_elim(MAT ,MAT *,Q *,int **,int **);
19:
20: int gauss_elim_mod(int **,int,int,int);
21: int gauss_elim_mod1(int **,int,int,int);
22: int gauss_elim_geninv_mod(unsigned int **,int,int,int);
23: int gauss_elim_geninv_mod_swap(unsigned int **,int,int,unsigned int,unsigned int ***,int **);
24: void Pnewvect(), Pnewmat(), Psepvect(), Psize(), Pdet(), Pleqm(), Pleqm1(), Pgeninvm();
25:
26: void Pgeneric_gauss_elim_mod();
27:
28: void Pmat_to_gfmmat(),Plu_gfmmat(),Psolve_by_lu_gfmmat();
29: void Pgeninvm_swap(), Premainder(), Psremainder(), Pvtol();
30: void sepvect();
31: void Pmulmat_gf2n();
32: void Pbconvmat_gf2n();
33: void Pmul_vect_mat_gf2n();
34: void PNBmul_gf2n();
35: void Pmul_mat_vect_int();
36: void Psepmat_destructive();
37: void Px962_irredpoly_up2();
38: void Pirredpoly_up2();
39: void Pnbpoly_up2();
40: void Pqsort();
41:
42: struct ftab array_tab[] = {
43: {"solve_by_lu_gfmmat",Psolve_by_lu_gfmmat,4},
44: {"lu_gfmmat",Plu_gfmmat,2},
45: {"mat_to_gfmmat",Pmat_to_gfmmat,2},
46: {"generic_gauss_elim_mod",Pgeneric_gauss_elim_mod,2},
47: {"newvect",Pnewvect,-2},
48: {"newmat",Pnewmat,-3},
49: {"sepmat_destructive",Psepmat_destructive,2},
50: {"sepvect",Psepvect,2},
51: {"qsort",Pqsort,-2},
52: {"vtol",Pvtol,1},
53: {"size",Psize,1},
54: {"det",Pdet,-2},
55: {"leqm",Pleqm,2},
56: {"leqm1",Pleqm1,2},
57: {"geninvm",Pgeninvm,2},
58: {"geninvm_swap",Pgeninvm_swap,2},
59: {"remainder",Premainder,2},
60: {"sremainder",Psremainder,2},
61: {"mulmat_gf2n",Pmulmat_gf2n,1},
62: {"bconvmat_gf2n",Pbconvmat_gf2n,-4},
63: {"mul_vect_mat_gf2n",Pmul_vect_mat_gf2n,2},
64: {"mul_mat_vect_int",Pmul_mat_vect_int,2},
65: {"nbmul_gf2n",PNBmul_gf2n,3},
66: {"x962_irredpoly_up2",Px962_irredpoly_up2,2},
67: {"irredpoly_up2",Pirredpoly_up2,2},
68: {"nbpoly_up2",Pnbpoly_up2,2},
69: {0,0,0},
70: };
71:
72: int comp_obj(a,b)
73: Obj *a,*b;
74: {
75: return arf_comp(CO,*a,*b);
76: }
77:
78: static FUNC generic_comp_obj_func;
79: static NODE generic_comp_obj_arg;
80:
81: int generic_comp_obj(a,b)
82: Obj *a,*b;
83: {
84: Q r;
85:
86: BDY(generic_comp_obj_arg)=(pointer)(*a);
87: BDY(NEXT(generic_comp_obj_arg))=(pointer)(*b);
88: r = (Q)bevalf(generic_comp_obj_func,generic_comp_obj_arg);
89: if ( !r )
90: return 0;
91: else
92: return SGN(r)>0?1:-1;
93: }
94:
95:
96: void Pqsort(arg,rp)
97: NODE arg;
98: VECT *rp;
99: {
100: VECT vect;
101: char buf[BUFSIZ];
102: char *fname;
103: NODE n;
104: P p;
105: V v;
106:
107: asir_assert(ARG0(arg),O_VECT,"qsort");
108: vect = (VECT)ARG0(arg);
109: if ( argc(arg) == 1 )
110: qsort(BDY(vect),vect->len,sizeof(Obj),(int (*)(const void *,const void *))comp_obj);
111: else {
112: p = (P)ARG1(arg);
113: if ( !p || OID(p)!=2 )
114: error("qsort : invalid argument");
115: v = VR(p);
116: if ( (int)v->attr != V_SR )
117: error("qsort : no such function");
118: generic_comp_obj_func = (FUNC)v->priv;
119: MKNODE(n,0,0); MKNODE(generic_comp_obj_arg,0,n);
120: qsort(BDY(vect),vect->len,sizeof(Obj),(int (*)(const void *,const void *))generic_comp_obj);
121: }
122: *rp = vect;
123: }
124:
125: void PNBmul_gf2n(arg,rp)
126: NODE arg;
127: GF2N *rp;
128: {
129: GF2N a,b;
130: GF2MAT mat;
131: int n,w;
132: unsigned int *ab,*bb;
133: UP2 r;
134:
135: a = (GF2N)ARG0(arg);
136: b = (GF2N)ARG1(arg);
137: mat = (GF2MAT)ARG2(arg);
138: if ( !a || !b )
139: *rp = 0;
140: else {
141: n = mat->row;
142: w = (n+BSH-1)/BSH;
143:
144: ab = (unsigned int *)ALLOCA(w*sizeof(unsigned int));
145: bzero((char *)ab,w*sizeof(unsigned int));
146: bcopy(a->body->b,ab,(a->body->w)*sizeof(unsigned int));
147:
148: bb = (unsigned int *)ALLOCA(w*sizeof(unsigned int));
149: bzero((char *)bb,w*sizeof(unsigned int));
150: bcopy(b->body->b,bb,(b->body->w)*sizeof(unsigned int));
151:
152: NEWUP2(r,w);
153: bzero((char *)r->b,w*sizeof(unsigned int));
154: mul_nb(mat,ab,bb,r->b);
155: r->w = w;
156: _adjup2(r);
157: if ( !r->w )
158: *rp = 0;
159: else
160: MKGF2N(r,*rp);
161: }
162: }
163:
164: void Pmul_vect_mat_gf2n(arg,rp)
165: NODE arg;
166: GF2N *rp;
167: {
168: GF2N a;
169: GF2MAT mat;
170: int n,w;
171: unsigned int *b;
172: UP2 r;
173:
174: a = (GF2N)ARG0(arg);
175: mat = (GF2MAT)ARG1(arg);
176: if ( !a )
177: *rp = 0;
178: else {
179: n = mat->row;
180: w = (n+BSH-1)/BSH;
181: b = (unsigned int *)ALLOCA(w*sizeof(unsigned int));
182: bzero((char *)b,w*sizeof(unsigned int));
183: bcopy(a->body->b,b,(a->body->w)*sizeof(unsigned int));
184: NEWUP2(r,w);
185: bzero((char *)r->b,w*sizeof(unsigned int));
186: mulgf2vectmat(mat->row,b,mat->body,r->b);
187: r->w = w;
188: _adjup2(r);
189: if ( !r->w )
190: *rp = 0;
191: else {
192: MKGF2N(r,*rp);
193: }
194: }
195: }
196:
197: void Pbconvmat_gf2n(arg,rp)
198: NODE arg;
199: LIST *rp;
200: {
201: P p0,p1;
202: int to;
203: GF2MAT p01,p10;
204: GF2N root;
205: NODE n0,n1;
206:
207: p0 = (P)ARG0(arg);
208: p1 = (P)ARG1(arg);
209: to = ARG2(arg)?1:0;
210: if ( argc(arg) == 4 ) {
211: root = (GF2N)ARG3(arg);
212: compute_change_of_basis_matrix_with_root(p0,p1,to,root,&p01,&p10);
213: } else
214: compute_change_of_basis_matrix(p0,p1,to,&p01,&p10);
215: MKNODE(n1,p10,0); MKNODE(n0,p01,n1);
216: MKLIST(*rp,n0);
217: }
218:
219: void Pmulmat_gf2n(arg,rp)
220: NODE arg;
221: GF2MAT *rp;
222: {
223: GF2MAT m;
224:
225: if ( !compute_multiplication_matrix((P)ARG0(arg),&m) )
226: error("mulmat_gf2n : input is not a normal polynomial");
227: *rp = m;
228: }
229:
230: void Psepmat_destructive(arg,rp)
231: NODE arg;
232: LIST *rp;
233: {
234: MAT mat,mat1;
235: int i,j,row,col;
236: Q **a,**a1;
237: Q ent;
238: N nm,mod,rem,quo;
239: int sgn;
240: NODE n0,n1;
241:
242: mat = (MAT)ARG0(arg); mod = NM((Q)ARG1(arg));
243: row = mat->row; col = mat->col;
244: MKMAT(mat1,row,col);
245: a = (Q **)mat->body; a1 = (Q **)mat1->body;
246: for ( i = 0; i < row; i++ )
247: for ( j = 0; j < col; j++ ) {
248: ent = a[i][j];
249: if ( !ent )
250: continue;
251: nm = NM(ent);
252: sgn = SGN(ent);
253: divn(nm,mod,&quo,&rem);
254: /* if ( quo != nm && rem != nm ) */
255: /* GC_free(nm); */
256: /* GC_free(ent); */
257: NTOQ(rem,sgn,a[i][j]); NTOQ(quo,sgn,a1[i][j]);
258: }
259: MKNODE(n1,mat1,0); MKNODE(n0,mat,n1);
260: MKLIST(*rp,n0);
261: }
262:
263: void Psepvect(arg,rp)
264: NODE arg;
265: VECT *rp;
266: {
267: sepvect((VECT)ARG0(arg),QTOS((Q)ARG1(arg)),rp);
268: }
269:
270: void sepvect(v,d,rp)
271: VECT v;
272: int d;
273: VECT *rp;
274: {
275: int i,j,k,n,q,q1,r;
276: pointer *pv,*pw,*pu;
277: VECT w,u;
278:
279: n = v->len;
280: if ( d > n )
281: d = n;
282: q = n/d; r = n%d; q1 = q+1;
283: MKVECT(w,d); *rp = w;
284: pv = BDY(v); pw = BDY(w); k = 0;
285: for ( i = 0; i < r; i++ ) {
286: MKVECT(u,q1); pw[i] = (pointer)u;
287: for ( pu = BDY(u), j = 0; j < q1; j++, k++ )
288: pu[j] = pv[k];
289: }
290: for ( ; i < d; i++ ) {
291: MKVECT(u,q); pw[i] = (pointer)u;
292: for ( pu = BDY(u), j = 0; j < q; j++, k++ )
293: pu[j] = pv[k];
294: }
295: }
296:
297: void Pnewvect(arg,rp)
298: NODE arg;
299: VECT *rp;
300: {
301: int len,i,r;
302: VECT vect;
303: pointer *vb;
304: LIST list;
305: NODE tn;
306:
307: asir_assert(ARG0(arg),O_N,"newvect");
308: len = QTOS((Q)ARG0(arg));
309: if ( len <= 0 )
310: error("newvect : invalid size");
311: MKVECT(vect,len);
312: if ( argc(arg) == 2 ) {
313: list = (LIST)ARG1(arg);
314: asir_assert(list,O_LIST,"newvect");
315: for ( r = 0, tn = BDY(list); tn; r++, tn = NEXT(tn) );
316: if ( r > len ) {
317: *rp = vect;
318: return;
319: }
320: for ( i = 0, tn = BDY(list), vb = BDY(vect); tn; i++, tn = NEXT(tn) )
321: vb[i] = (pointer)BDY(tn);
322: }
323: *rp = vect;
324: }
325:
326: void Pnewmat(arg,rp)
327: NODE arg;
328: MAT *rp;
329: {
330: int row,col;
331: int i,j,r,c;
332: NODE tn,sn;
333: MAT m;
334: pointer **mb;
335: LIST list;
336:
337: asir_assert(ARG0(arg),O_N,"newmat");
338: asir_assert(ARG1(arg),O_N,"newmat");
339: row = QTOS((Q)ARG0(arg)); col = QTOS((Q)ARG1(arg));
340: if ( row <= 0 || col <= 0 )
341: error("newmat : invalid size");
342: MKMAT(m,row,col);
343: if ( argc(arg) == 3 ) {
344: list = (LIST)ARG2(arg);
345: asir_assert(list,O_LIST,"newmat");
346: for ( r = 0, c = 0, tn = BDY(list); tn; r++, tn = NEXT(tn) ) {
347: for ( j = 0, sn = BDY((LIST)BDY(tn)); sn; j++, sn = NEXT(sn) );
348: c = MAX(c,j);
349: }
350: if ( (r > row) || (c > col) ) {
351: *rp = m;
352: return;
353: }
354: for ( i = 0, tn = BDY(list), mb = BDY(m); tn; i++, tn = NEXT(tn) ) {
355: asir_assert(BDY(tn),O_LIST,"newmat");
356: for ( j = 0, sn = BDY((LIST)BDY(tn)); sn; j++, sn = NEXT(sn) )
357: mb[i][j] = (pointer)BDY(sn);
358: }
359: }
360: *rp = m;
361: }
362:
363: void Pvtol(arg,rp)
364: NODE arg;
365: LIST *rp;
366: {
367: NODE n,n1;
368: VECT v;
369: pointer *a;
370: int len,i;
371:
372: asir_assert(ARG0(arg),O_VECT,"vtol");
373: v = (VECT)ARG0(arg); len = v->len; a = BDY(v);
374: for ( i = len - 1, n = 0; i >= 0; i-- ) {
375: MKNODE(n1,a[i],n); n = n1;
376: }
377: MKLIST(*rp,n);
378: }
379:
380: void Premainder(arg,rp)
381: NODE arg;
382: Obj *rp;
383: {
384: Obj a;
385: VECT v,w;
386: MAT m,l;
387: pointer *vb,*wb;
388: pointer **mb,**lb;
389: int id,i,j,n,row,col,t,smd,sgn;
390: Q md,q;
391:
392: a = (Obj)ARG0(arg); md = (Q)ARG1(arg);
393: if ( !a )
394: *rp = 0;
395: else {
396: id = OID(a);
397: switch ( id ) {
398: case O_N:
399: case O_P:
400: cmp(md,(P)a,(P *)rp); break;
401: case O_VECT:
402: smd = QTOS(md);
403: v = (VECT)a; n = v->len; vb = v->body;
404: MKVECT(w,n); wb = w->body;
405: for ( i = 0; i < n; i++ ) {
406: if ( q = (Q)vb[i] ) {
407: sgn = SGN(q); t = rem(NM(q),smd);
408: STOQ(t,q);
409: if ( q )
410: SGN(q) = sgn;
411: }
412: wb[i] = (pointer)q;
413: }
414: *rp = (Obj)w;
415: break;
416: case O_MAT:
417: m = (MAT)a; row = m->row; col = m->col; mb = m->body;
418: MKMAT(l,row,col); lb = l->body;
419: for ( i = 0; i < row; i++ )
420: for ( j = 0, vb = mb[i], wb = lb[i]; j < col; j++ )
421: cmp(md,(P)vb[j],(P *)&wb[j]);
422: *rp = (Obj)l;
423: break;
424: default:
425: error("remainder : invalid argument");
426: }
427: }
428: }
429:
430: void Psremainder(arg,rp)
431: NODE arg;
432: Obj *rp;
433: {
434: Obj a;
435: VECT v,w;
436: MAT m,l;
437: pointer *vb,*wb;
438: pointer **mb,**lb;
439: unsigned int t,smd;
440: int id,i,j,n,row,col;
441: Q md,q;
442:
443: a = (Obj)ARG0(arg); md = (Q)ARG1(arg);
444: if ( !a )
445: *rp = 0;
446: else {
447: id = OID(a);
448: switch ( id ) {
449: case O_N:
450: case O_P:
451: cmp(md,(P)a,(P *)rp); break;
452: case O_VECT:
453: smd = QTOS(md);
454: v = (VECT)a; n = v->len; vb = v->body;
455: MKVECT(w,n); wb = w->body;
456: for ( i = 0; i < n; i++ ) {
457: if ( q = (Q)vb[i] ) {
458: t = (unsigned int)rem(NM(q),smd);
459: if ( SGN(q) < 0 )
460: t = (smd - t) % smd;
461: UTOQ(t,q);
462: }
463: wb[i] = (pointer)q;
464: }
465: *rp = (Obj)w;
466: break;
467: case O_MAT:
468: m = (MAT)a; row = m->row; col = m->col; mb = m->body;
469: MKMAT(l,row,col); lb = l->body;
470: for ( i = 0; i < row; i++ )
471: for ( j = 0, vb = mb[i], wb = lb[i]; j < col; j++ )
472: cmp(md,(P)vb[j],(P *)&wb[j]);
473: *rp = (Obj)l;
474: break;
475: default:
476: error("remainder : invalid argument");
477: }
478: }
479: }
480:
481: void Psize(arg,rp)
482: NODE arg;
483: LIST *rp;
484: {
485:
486: int n,m;
487: Q q;
488: NODE t,s;
489:
490: if ( !ARG0(arg) )
491: t = 0;
492: else {
493: switch (OID(ARG0(arg))) {
494: case O_VECT:
495: n = ((VECT)ARG0(arg))->len;
496: STOQ(n,q); MKNODE(t,q,0);
497: break;
498: case O_MAT:
499: n = ((MAT)ARG0(arg))->row; m = ((MAT)ARG0(arg))->col;
500: STOQ(m,q); MKNODE(s,q,0); STOQ(n,q); MKNODE(t,q,s);
501: break;
502: default:
503: error("size : invalid argument"); break;
504: }
505: }
506: MKLIST(*rp,t);
507: }
508:
509: void Pdet(arg,rp)
510: NODE arg;
511: P *rp;
512: {
513: MAT m;
514: int n,i,j,mod;
515: P d;
516: P **mat,**w;
517:
518: m = (MAT)ARG0(arg);
519: asir_assert(m,O_MAT,"det");
520: if ( m->row != m->col )
521: error("det : non-square matrix");
522: else if ( argc(arg) == 1 )
523: detp(CO,(P **)BDY(m),m->row,rp);
524: else {
525: n = m->row; mod = QTOS((Q)ARG1(arg)); mat = (P **)BDY(m);
526: w = (P **)almat_pointer(n,n);
527: for ( i = 0; i < n; i++ )
528: for ( j = 0; j < n; j++ )
529: ptomp(mod,mat[i][j],&w[i][j]);
530: detmp(CO,mod,w,n,&d);
531: mptop(d,rp);
532: }
533: }
534:
535: /*
536: input : a row x col matrix A
537: A[I] <-> A[I][0]*x_0+A[I][1]*x_1+...
538:
539: output : [B,R,C]
540: B : a rank(A) x col-rank(A) matrix
541: R : a vector of length rank(A)
542: C : a vector of length col-rank(A)
543: B[I] <-> x_{R[I]}+B[I][0]x_{C[0]}+B[I][1]x_{C[1]}+...
544: */
545:
546: void Pgeneric_gauss_elim_mod(arg,rp)
547: NODE arg;
548: LIST *rp;
549: {
550: NODE n0;
551: MAT m,mat;
552: VECT rind,cind;
553: Q **tmat;
554: int **wmat;
555: Q *rib,*cib;
556: int *colstat;
557: Q q;
558: int md,i,j,k,l,row,col,t,n,rank;
559:
560: asir_assert(ARG0(arg),O_MAT,"generic_gauss_elim_mod");
561: asir_assert(ARG1(arg),O_N,"generic_gauss_elim_mod");
562: m = (MAT)ARG0(arg); md = QTOS((Q)ARG1(arg));
563: row = m->row; col = m->col; tmat = (Q **)m->body;
564: wmat = (int **)almat(row,col);
565: colstat = (int *)MALLOC_ATOMIC(col*sizeof(int));
566: for ( i = 0; i < row; i++ )
567: for ( j = 0; j < col; j++ )
568: if ( q = (Q)tmat[i][j] ) {
569: t = rem(NM(q),md);
570: if ( t && SGN(q) < 0 )
571: t = (md - t) % md;
572: wmat[i][j] = t;
573: } else
574: wmat[i][j] = 0;
575: rank = generic_gauss_elim_mod(wmat,row,col,md,colstat);
576:
577: MKMAT(mat,rank,col-rank);
578: tmat = (Q **)mat->body;
579: for ( i = 0; i < rank; i++ )
580: for ( j = k = 0; j < col; j++ )
581: if ( !colstat[j] ) {
582: UTOQ(wmat[i][j],tmat[i][k]); k++;
583: }
584:
585: MKVECT(rind,rank);
586: MKVECT(cind,col-rank);
587: rib = (Q *)rind->body; cib = (Q *)cind->body;
588: for ( j = k = l = 0; j < col; j++ )
589: if ( colstat[j] ) {
590: STOQ(j,rib[k]); k++;
591: } else {
592: STOQ(j,cib[l]); l++;
593: }
594: n0 = mknode(3,mat,rind,cind);
595: MKLIST(*rp,n0);
596: }
597:
598: void Pleqm(arg,rp)
599: NODE arg;
600: VECT *rp;
601: {
602: MAT m;
603: VECT vect;
604: pointer **mat;
605: Q *v;
606: Q q;
607: int **wmat;
608: int md,i,j,row,col,t,n,status;
609:
610: asir_assert(ARG0(arg),O_MAT,"leqm");
611: asir_assert(ARG1(arg),O_N,"leqm");
612: m = (MAT)ARG0(arg); md = QTOS((Q)ARG1(arg));
613: row = m->row; col = m->col; mat = m->body;
614: wmat = (int **)almat(row,col);
615: for ( i = 0; i < row; i++ )
616: for ( j = 0; j < col; j++ )
617: if ( q = (Q)mat[i][j] ) {
618: t = rem(NM(q),md);
619: if ( SGN(q) < 0 )
620: t = (md - t) % md;
621: wmat[i][j] = t;
622: } else
623: wmat[i][j] = 0;
624: status = gauss_elim_mod(wmat,row,col,md);
625: if ( status < 0 )
626: *rp = 0;
627: else if ( status > 0 )
628: *rp = (VECT)ONE;
629: else {
630: n = col - 1;
631: MKVECT(vect,n);
632: for ( i = 0, v = (Q *)vect->body; i < n; i++ ) {
633: t = (md-wmat[i][n])%md; STOQ(t,v[i]);
634: }
635: *rp = vect;
636: }
637: }
638:
639: int gauss_elim_mod(mat,row,col,md)
640: int **mat;
641: int row,col,md;
642: {
643: int i,j,k,inv,a,n;
644: int *t,*pivot;
645:
646: n = col - 1;
647: for ( j = 0; j < n; j++ ) {
648: for ( i = j; i < row && !mat[i][j]; i++ );
649: if ( i == row )
650: return 1;
651: if ( i != j ) {
652: t = mat[i]; mat[i] = mat[j]; mat[j] = t;
653: }
654: pivot = mat[j];
655: inv = invm(pivot[j],md);
656: for ( k = j; k <= n; k++ ) {
657: /* pivot[k] = dmar(pivot[k],inv,0,md); */
658: DMAR(pivot[k],inv,0,md,pivot[k])
659: }
660: for ( i = 0; i < row; i++ ) {
661: t = mat[i];
662: if ( i != j && (a = t[j]) )
663: for ( k = j, a = md - a; k <= n; k++ ) {
664: /* t[k] = dmar(pivot[k],a,t[k],md); */
665: DMAR(pivot[k],a,t[k],md,t[k])
666: }
667: }
668: }
669: for ( i = n; i < row && !mat[i][n]; i++ );
670: if ( i == row )
671: return 0;
672: else
673: return -1;
674: }
675:
676: struct oEGT eg_mod,eg_elim,eg_chrem,eg_gschk,eg_intrat,eg_symb;
677:
678: int generic_gauss_elim(mat,nm,dn,rindp,cindp)
679: MAT mat;
680: MAT *nm;
681: Q *dn;
682: int **rindp,**cindp;
683: {
684: int **wmat;
685: Q **bmat;
686: N **tmat;
687: Q *bmi;
688: N *tmi;
689: Q q;
690: int *wmi;
691: int *colstat,*wcolstat,*rind,*cind;
692: int row,col,ind,md,i,j,k,l,t,t1,rank,rank0,inv;
693: N m1,m2,m3,s,u;
694: MAT r,crmat;
695: struct oEGT tmp0,tmp1;
696: struct oEGT eg_mod_split,eg_elim_split,eg_chrem_split;
697: struct oEGT eg_intrat_split,eg_gschk_split;
698: int ret;
699:
700: init_eg(&eg_mod_split); init_eg(&eg_chrem_split);
701: init_eg(&eg_elim_split); init_eg(&eg_intrat_split);
702: init_eg(&eg_gschk_split);
703: bmat = (Q **)mat->body;
704: row = mat->row; col = mat->col;
705: wmat = (int **)almat(row,col);
706: colstat = (int *)MALLOC_ATOMIC(col*sizeof(int));
707: wcolstat = (int *)MALLOC_ATOMIC(col*sizeof(int));
708: for ( ind = 0; ; ind++ ) {
709: if ( Print )
710: fprintf(asir_out,".");
711: md = lprime[ind];
712: get_eg(&tmp0);
713: for ( i = 0; i < row; i++ )
714: for ( j = 0, bmi = bmat[i], wmi = wmat[i]; j < col; j++ )
715: if ( q = (Q)bmi[j] ) {
716: t = rem(NM(q),md);
717: if ( t && SGN(q) < 0 )
718: t = (md - t) % md;
719: wmi[j] = t;
720: } else
721: wmi[j] = 0;
722: get_eg(&tmp1);
723: add_eg(&eg_mod,&tmp0,&tmp1);
724: add_eg(&eg_mod_split,&tmp0,&tmp1);
725: get_eg(&tmp0);
726: rank = generic_gauss_elim_mod(wmat,row,col,md,wcolstat);
727: get_eg(&tmp1);
728: add_eg(&eg_elim,&tmp0,&tmp1);
729: add_eg(&eg_elim_split,&tmp0,&tmp1);
730: if ( !ind ) {
731: RESET:
732: UTON(md,m1);
733: rank0 = rank;
734: bcopy(wcolstat,colstat,col*sizeof(int));
735: MKMAT(crmat,rank,col-rank);
736: MKMAT(r,rank,col-rank); *nm = r;
737: tmat = (N **)crmat->body;
738: for ( i = 0; i < rank; i++ )
739: for ( j = k = 0, tmi = tmat[i], wmi = wmat[i]; j < col; j++ )
740: if ( !colstat[j] ) {
741: UTON(wmi[j],tmi[k]); k++;
742: }
743: } else {
744: if ( rank < rank0 ) {
745: if ( Print )
746: fprintf(asir_out,"lower rank matrix; continuing...\n");
747: continue;
748: } else if ( rank > rank0 ) {
749: if ( Print )
750: fprintf(asir_out,"higher rank matrix; resetting...\n");
751: goto RESET;
752: } else {
753: for ( j = 0; (j<col) && (colstat[j]==wcolstat[j]); j++ );
754: if ( j < col ) {
755: if ( Print )
756: fprintf(asir_out,"inconsitent colstat; resetting...\n");
757: goto RESET;
758: }
759: }
760:
761: get_eg(&tmp0);
762: inv = invm(rem(m1,md),md);
763: UTON(md,m2); muln(m1,m2,&m3);
764: for ( i = 0; i < rank; i++ )
765: for ( j = k = 0, tmi = tmat[i], wmi = wmat[i]; j < col; j++ )
766: if ( !colstat[j] ) {
767: if ( tmi[k] ) {
768: /* f3 = f1+m1*(m1 mod m2)^(-1)*(f2 - f1 mod m2) */
769: t = rem(tmi[k],md);
770: if ( wmi[j] >= t )
771: t = wmi[j]-t;
772: else
773: t = md-(t-wmi[j]);
774: DMAR(t,inv,0,md,t1)
775: UTON(t1,u);
776: muln(m1,u,&s);
777: addn(tmi[k],s,&u); tmi[k] = u;
778: } else if ( wmi[j] ) {
779: /* f3 = m1*(m1 mod m2)^(-1)*f2 */
780: DMAR(wmi[j],inv,0,md,t)
781: UTON(t,u);
782: muln(m1,u,&s); tmi[k] = s;
783: }
784: k++;
785: }
786: m1 = m3;
787: get_eg(&tmp1);
788: add_eg(&eg_chrem,&tmp0,&tmp1);
789: add_eg(&eg_chrem_split,&tmp0,&tmp1);
790:
791: get_eg(&tmp0);
792: ret = intmtoratm(crmat,m1,*nm,dn);
793: get_eg(&tmp1);
794: add_eg(&eg_intrat,&tmp0,&tmp1);
795: add_eg(&eg_intrat_split,&tmp0,&tmp1);
796: if ( ret ) {
797: *rindp = rind = (int *)MALLOC_ATOMIC(rank*sizeof(int));
798: *cindp = cind = (int *)MALLOC_ATOMIC((col-rank)*sizeof(int));
799: for ( j = k = l = 0; j < col; j++ )
800: if ( colstat[j] )
801: rind[k++] = j;
802: else
803: cind[l++] = j;
804: get_eg(&tmp0);
805: if ( gensolve_check(mat,*nm,*dn,rind,cind) )
806: get_eg(&tmp1);
807: add_eg(&eg_gschk,&tmp0,&tmp1);
808: add_eg(&eg_gschk_split,&tmp0,&tmp1);
809: if ( Print ) {
810: print_eg("Mod",&eg_mod_split);
811: print_eg("Elim",&eg_elim_split);
812: print_eg("ChRem",&eg_chrem_split);
813: print_eg("IntRat",&eg_intrat_split);
814: print_eg("Check",&eg_gschk_split);
815: fflush(asir_out);
816: }
817: return rank;
818: }
819: }
820: }
821: }
822:
823: int f4_nocheck;
824:
825: int gensolve_check(mat,nm,dn,rind,cind)
826: MAT mat,nm;
827: Q dn;
828: int *rind,*cind;
829: {
830: int row,col,rank,clen,i,j,k,l;
831: Q s,t,u;
832: Q *w;
833: Q *mati,*nmk;
834:
835: if ( f4_nocheck )
836: return 1;
837: row = mat->row; col = mat->col;
838: rank = nm->row; clen = nm->col;
839: w = (Q *)MALLOC(clen*sizeof(Q));
840: for ( i = 0; i < row; i++ ) {
841: mati = (Q *)mat->body[i];
842: #if 1
843: bzero(w,clen*sizeof(Q));
844: for ( k = 0; k < rank; k++ )
845: for ( l = 0, nmk = (Q *)nm->body[k]; l < clen; l++ ) {
846: mulq(mati[rind[k]],nmk[l],&t);
847: addq(w[l],t,&s); w[l] = s;
848: }
849: for ( j = 0; j < clen; j++ ) {
850: mulq(dn,mati[cind[j]],&t);
851: if ( cmpq(w[j],t) )
852: break;
853: }
854: #else
855: for ( j = 0; j < clen; j++ ) {
856: for ( k = 0, s = 0; k < rank; k++ ) {
857: mulq(mati[rind[k]],nm->body[k][j],&t);
858: addq(s,t,&u); s = u;
859: }
860: mulq(dn,mati[cind[j]],&t);
861: if ( cmpq(s,t) )
862: break;
863: }
864: #endif
865: if ( j != clen )
866: break;
867: }
868: if ( i != row )
869: return 0;
870: else
871: return 1;
872: }
873:
874: /* assuming 0 < c < m */
875:
876: int inttorat(c,m,b,sgnp,nmp,dnp)
877: N c,m,b;
878: int *sgnp;
879: N *nmp,*dnp;
880: {
881: Q qq,t,u1,v1,r1,nm;
882: N q,r,u2,v2,r2;
883:
884: u1 = 0; v1 = ONE; u2 = m; v2 = c;
885: while ( cmpn(v2,b) >= 0 ) {
886: divn(u2,v2,&q,&r2); u2 = v2; v2 = r2;
887: NTOQ(q,1,qq); mulq(qq,v1,&t); subq(u1,t,&r1); u1 = v1; v1 = r1;
888: }
889: if ( cmpn(NM(v1),b) >= 0 )
890: return 0;
891: else {
892: *nmp = v2;
893: *dnp = NM(v1);
894: *sgnp = SGN(v1);
895: return 1;
896: }
897: }
898:
899: /* mat->body = N ** */
900:
901: int intmtoratm(mat,md,nm,dn)
902: MAT mat;
903: N md;
904: MAT nm;
905: Q *dn;
906: {
907: N t,s,b;
908: Q bound,dn0,dn1,nm1,q,tq;
909: int i,j,k,l,row,col;
910: Q **rmat;
911: N **tmat;
912: N *tmi;
913: Q *nmk;
914: N u,unm,udn;
915: int sgn,ret;
916:
917: row = mat->row; col = mat->col;
918: bshiftn(md,1,&t);
919: isqrt(t,&s);
920: bshiftn(s,64,&b);
921: if ( !b )
922: b = ONEN;
923: dn0 = ONE;
924: tmat = (N **)mat->body;
925: rmat = (Q **)nm->body;
926: for ( i = 0; i < row; i++ )
927: for ( j = 0, tmi = tmat[i]; j < col; j++ )
928: if ( tmi[j] ) {
929: muln(tmi[j],NM(dn0),&s);
930: remn(s,md,&u);
931: ret = inttorat(u,md,b,&sgn,&unm,&udn);
932: if ( !ret )
933: return 0;
934: else {
935: NTOQ(unm,sgn,nm1);
936: NTOQ(udn,1,dn1);
937: if ( !UNIQ(dn1) ) {
938: for ( k = 0; k < i; k++ )
939: for ( l = 0, nmk = rmat[k]; l < col; l++ ) {
940: mulq(nmk[l],dn1,&q); nmk[l] = q;
941: }
942: for ( l = 0, nmk = rmat[i]; l < j; l++ ) {
943: mulq(nmk[l],dn1,&q); nmk[l] = q;
944: }
945: }
946: rmat[i][j] = nm1;
947: mulq(dn0,dn1,&q); dn0 = q;
948: }
949: }
950: *dn = dn0;
951: return 1;
952: }
953:
954: int generic_gauss_elim_mod(mat,row,col,md,colstat)
955: int **mat;
956: int row,col,md;
957: int *colstat;
958: {
959: int i,j,k,l,inv,a,rank;
960: int *t,*pivot;
961:
962: for ( rank = 0, j = 0; j < col; j++ ) {
963: for ( i = rank; i < row && !mat[i][j]; i++ );
964: if ( i == row ) {
965: colstat[j] = 0;
966: continue;
967: } else
968: colstat[j] = 1;
969: if ( i != rank ) {
970: t = mat[i]; mat[i] = mat[rank]; mat[rank] = t;
971: }
972: pivot = mat[rank];
973: inv = invm(pivot[j],md);
974: for ( k = j; k < col; k++ )
975: if ( pivot[k] ) {
976: DMAR(pivot[k],inv,0,md,pivot[k])
977: }
978: for ( i = rank+1; i < row; i++ ) {
979: t = mat[i];
980: if ( a = t[j] )
981: for ( k = j, a = md - a; k < col; k++ )
982: if ( pivot[k] ) {
983: DMAR(pivot[k],a,t[k],md,t[k])
984: }
985: }
986: rank++;
987: }
988: for ( j = col-1, l = rank-1; j >= 0; j-- )
989: if ( colstat[j] ) {
990: pivot = mat[l];
991: for ( i = 0; i < l; i++ ) {
992: t = mat[i];
993: if ( a = t[j] )
994: for ( k = j, a = md-a; k < col; k++ )
995: if ( pivot[k] ) {
996: DMAR(pivot[k],a,t[k],md,t[k])
997: }
998: }
999: l--;
1000: }
1001: return rank;
1002: }
1003:
1004: /* LU decomposition; a[i][i] = 1/U[i][i] */
1005:
1006: int lu_gfmmat(mat,md,perm)
1007: GFMMAT mat;
1008: unsigned int md;
1009: int *perm;
1010: {
1011: int row,col;
1012: int i,j,k,l;
1013: unsigned int *t,*pivot;
1014: unsigned int **a;
1015: unsigned int inv,m;
1016:
1017: row = mat->row; col = mat->col;
1018: a = mat->body;
1019: bzero(perm,row*sizeof(int));
1020:
1021: for ( i = 0; i < row; i++ )
1022: perm[i] = i;
1023: for ( k = 0; k < col; k++ ) {
1024: for ( i = k; i < row && !a[i][k]; i++ );
1025: if ( i == row )
1026: return 0;
1027: if ( i != k ) {
1028: j = perm[i]; perm[i] = perm[k]; perm[k] = j;
1029: t = a[i]; a[i] = a[k]; a[k] = t;
1030: }
1031: pivot = a[k];
1032: pivot[k] = inv = invm(pivot[k],md);
1033: for ( i = k+1; i < row; i++ ) {
1034: t = a[i];
1035: if ( m = t[k] ) {
1036: DMAR(inv,m,0,md,t[k])
1037: for ( j = k+1, m = md - t[k]; j < col; j++ )
1038: if ( pivot[j] ) {
1039: DMAR(m,pivot[j],t[j],md,t[j])
1040: }
1041: }
1042: }
1043: }
1044: return 1;
1045: }
1046:
1047: void Pleqm1(arg,rp)
1048: NODE arg;
1049: VECT *rp;
1050: {
1051: MAT m;
1052: VECT vect;
1053: pointer **mat;
1054: Q *v;
1055: Q q;
1056: int **wmat;
1057: int md,i,j,row,col,t,n,status;
1058:
1059: asir_assert(ARG0(arg),O_MAT,"leqm1");
1060: asir_assert(ARG1(arg),O_N,"leqm1");
1061: m = (MAT)ARG0(arg); md = QTOS((Q)ARG1(arg));
1062: row = m->row; col = m->col; mat = m->body;
1063: wmat = (int **)almat(row,col);
1064: for ( i = 0; i < row; i++ )
1065: for ( j = 0; j < col; j++ )
1066: if ( q = (Q)mat[i][j] ) {
1067: t = rem(NM(q),md);
1068: if ( SGN(q) < 0 )
1069: t = (md - t) % md;
1070: wmat[i][j] = t;
1071: } else
1072: wmat[i][j] = 0;
1073: status = gauss_elim_mod1(wmat,row,col,md);
1074: if ( status < 0 )
1075: *rp = 0;
1076: else if ( status > 0 )
1077: *rp = (VECT)ONE;
1078: else {
1079: n = col - 1;
1080: MKVECT(vect,n);
1081: for ( i = 0, v = (Q *)vect->body; i < n; i++ ) {
1082: t = (md-wmat[i][n])%md; STOQ(t,v[i]);
1083: }
1084: *rp = vect;
1085: }
1086: }
1087:
1088: gauss_elim_mod1(mat,row,col,md)
1089: int **mat;
1090: int row,col,md;
1091: {
1092: int i,j,k,inv,a,n;
1093: int *t,*pivot;
1094:
1095: n = col - 1;
1096: for ( j = 0; j < n; j++ ) {
1097: for ( i = j; i < row && !mat[i][j]; i++ );
1098: if ( i == row )
1099: return 1;
1100: if ( i != j ) {
1101: t = mat[i]; mat[i] = mat[j]; mat[j] = t;
1102: }
1103: pivot = mat[j];
1104: inv = invm(pivot[j],md);
1105: for ( k = j; k <= n; k++ )
1106: pivot[k] = dmar(pivot[k],inv,0,md);
1107: for ( i = j+1; i < row; i++ ) {
1108: t = mat[i];
1109: if ( i != j && (a = t[j]) )
1110: for ( k = j, a = md - a; k <= n; k++ )
1111: t[k] = dmar(pivot[k],a,t[k],md);
1112: }
1113: }
1114: for ( i = n; i < row && !mat[i][n]; i++ );
1115: if ( i == row ) {
1116: for ( j = n-1; j >= 0; j-- ) {
1117: for ( i = j-1, a = (md-mat[j][n])%md; i >= 0; i-- ) {
1118: mat[i][n] = dmar(mat[i][j],a,mat[i][n],md);
1119: mat[i][j] = 0;
1120: }
1121: }
1122: return 0;
1123: } else
1124: return -1;
1125: }
1126:
1127: void Pgeninvm(arg,rp)
1128: NODE arg;
1129: LIST *rp;
1130: {
1131: MAT m;
1132: pointer **mat;
1133: Q **tmat;
1134: Q q;
1135: unsigned int **wmat;
1136: int md,i,j,row,col,t,status;
1137: MAT mat1,mat2;
1138: NODE node1,node2;
1139:
1140: asir_assert(ARG0(arg),O_MAT,"leqm1");
1141: asir_assert(ARG1(arg),O_N,"leqm1");
1142: m = (MAT)ARG0(arg); md = QTOS((Q)ARG1(arg));
1143: row = m->row; col = m->col; mat = m->body;
1144: wmat = (unsigned int **)almat(row,col+row);
1145: for ( i = 0; i < row; i++ ) {
1146: bzero((char *)wmat[i],(col+row)*sizeof(int));
1147: for ( j = 0; j < col; j++ )
1148: if ( q = (Q)mat[i][j] ) {
1149: t = rem(NM(q),md);
1150: if ( SGN(q) < 0 )
1151: t = (md - t) % md;
1152: wmat[i][j] = t;
1153: }
1154: wmat[i][col+i] = 1;
1155: }
1156: status = gauss_elim_geninv_mod(wmat,row,col,md);
1157: if ( status > 0 )
1158: *rp = 0;
1159: else {
1160: MKMAT(mat1,col,row); MKMAT(mat2,row-col,row);
1161: for ( i = 0, tmat = (Q **)mat1->body; i < col; i++ )
1162: for ( j = 0; j < row; j++ )
1163: STOQ(wmat[i][j+col],tmat[i][j]);
1164: for ( tmat = (Q **)mat2->body; i < row; i++ )
1165: for ( j = 0; j < row; j++ )
1166: STOQ(wmat[i][j+col],tmat[i-col][j]);
1167: MKNODE(node2,mat2,0); MKNODE(node1,mat1,node2); MKLIST(*rp,node1);
1168: }
1169: }
1170:
1171: int gauss_elim_geninv_mod(mat,row,col,md)
1172: unsigned int **mat;
1173: int row,col,md;
1174: {
1175: int i,j,k,inv,a,n,m;
1176: unsigned int *t,*pivot;
1177:
1178: n = col; m = row+col;
1179: for ( j = 0; j < n; j++ ) {
1180: for ( i = j; i < row && !mat[i][j]; i++ );
1181: if ( i == row )
1182: return 1;
1183: if ( i != j ) {
1184: t = mat[i]; mat[i] = mat[j]; mat[j] = t;
1185: }
1186: pivot = mat[j];
1187: inv = invm(pivot[j],md);
1188: for ( k = j; k < m; k++ )
1189: pivot[k] = dmar(pivot[k],inv,0,md);
1190: for ( i = j+1; i < row; i++ ) {
1191: t = mat[i];
1192: if ( a = t[j] )
1193: for ( k = j, a = md - a; k < m; k++ )
1194: t[k] = dmar(pivot[k],a,t[k],md);
1195: }
1196: }
1197: for ( j = n-1; j >= 0; j-- ) {
1198: pivot = mat[j];
1199: for ( i = j-1; i >= 0; i-- ) {
1200: t = mat[i];
1201: if ( a = t[j] )
1202: for ( k = j, a = md - a; k < m; k++ )
1203: t[k] = dmar(pivot[k],a,t[k],md);
1204: }
1205: }
1206: return 0;
1207: }
1208:
1209: void Psolve_by_lu_gfmmat(arg,rp)
1210: NODE arg;
1211: VECT *rp;
1212: {
1213: GFMMAT lu;
1214: Q *perm,*rhs,*v;
1215: int n,i;
1216: unsigned int md;
1217: unsigned int *b,*sol;
1218: VECT r;
1219:
1220: lu = (GFMMAT)ARG0(arg);
1221: perm = (Q *)BDY((VECT)ARG1(arg));
1222: rhs = (Q *)BDY((VECT)ARG2(arg));
1223: md = (unsigned int)QTOS((Q)ARG3(arg));
1224: n = lu->col;
1225: b = (unsigned int *)MALLOC_ATOMIC(n*sizeof(int));
1226: sol = (unsigned int *)MALLOC_ATOMIC(n*sizeof(int));
1227: for ( i = 0; i < n; i++ )
1228: b[i] = QTOS(rhs[QTOS(perm[i])]);
1229: solve_by_lu_gfmmat(lu,md,b,sol);
1230: MKVECT(r,n);
1231: for ( i = 0, v = (Q *)r->body; i < n; i++ )
1232: STOQ(sol[i],v[i]);
1233: *rp = r;
1234: }
1235:
1236: void solve_by_lu_gfmmat(lu,md,b,x)
1237: GFMMAT lu;
1238: unsigned int md;
1239: unsigned int *b;
1240: unsigned int *x;
1241: {
1242: int n;
1243: unsigned int **a;
1244: unsigned int *y;
1245: int i,j;
1246: unsigned int t,m;
1247:
1248: n = lu->col;
1249: a = lu->body;
1250: y = (unsigned int *)MALLOC_ATOMIC(n*sizeof(int));
1251: /* solve Ly=b */
1252: for ( i = 0; i < n; i++ ) {
1253: for ( t = b[i], j = 0; j < i; j++ )
1254: if ( a[i][j] ) {
1255: m = md - a[i][j];
1256: DMAR(m,y[j],t,md,t)
1257: }
1258: y[i] = t;
1259: }
1260: /* solve Ux=y */
1261: for ( i = n-1; i >= 0; i-- ) {
1262: for ( t = y[i], j =i+1; j < n; j++ )
1263: if ( a[i][j] ) {
1264: m = md - a[i][j];
1265: DMAR(m,x[j],t,md,t)
1266: }
1267: /* a[i][i] = 1/U[i][i] */
1268: DMAR(t,a[i][i],0,md,x[i])
1269: }
1270: }
1271:
1272: void Plu_gfmmat(arg,rp)
1273: NODE arg;
1274: LIST *rp;
1275: {
1276: MAT m;
1277: GFMMAT mm;
1278: unsigned int md;
1279: int i,row,col,status;
1280: int *iperm;
1281: Q *v;
1282: VECT perm;
1283: NODE n0;
1284:
1285: asir_assert(ARG0(arg),O_MAT,"mat_to_gfmmat");
1286: asir_assert(ARG1(arg),O_N,"mat_to_gfmmat");
1287: m = (MAT)ARG0(arg); md = (unsigned int)QTOS((Q)ARG1(arg));
1288: mat_to_gfmmat(m,md,&mm);
1289: row = m->row;
1290: col = m->col;
1291: iperm = (int *)MALLOC_ATOMIC(row*sizeof(int));
1292: status = lu_gfmmat(mm,md,iperm);
1293: if ( !status )
1294: n0 = 0;
1295: else {
1296: MKVECT(perm,row);
1297: for ( i = 0, v = (Q *)perm->body; i < row; i++ )
1298: STOQ(iperm[i],v[i]);
1299: n0 = mknode(2,mm,perm);
1300: }
1301: MKLIST(*rp,n0);
1302: }
1303:
1304: void Pmat_to_gfmmat(arg,rp)
1305: NODE arg;
1306: GFMMAT *rp;
1307: {
1308: MAT m;
1309: unsigned int md;
1310:
1311: asir_assert(ARG0(arg),O_MAT,"mat_to_gfmmat");
1312: asir_assert(ARG1(arg),O_N,"mat_to_gfmmat");
1313: m = (MAT)ARG0(arg); md = (unsigned int)QTOS((Q)ARG1(arg));
1314: mat_to_gfmmat(m,md,rp);
1315: }
1316:
1317: void mat_to_gfmmat(m,md,rp)
1318: MAT m;
1319: unsigned int md;
1320: GFMMAT *rp;
1321: {
1322: unsigned int **wmat;
1323: unsigned int t;
1324: Q **mat;
1325: Q q;
1326: int i,j,row,col;
1327:
1328: row = m->row; col = m->col; mat = (Q **)m->body;
1329: wmat = (unsigned int **)almat(row,col);
1330: for ( i = 0; i < row; i++ ) {
1331: bzero((char *)wmat[i],col*sizeof(unsigned int));
1332: for ( j = 0; j < col; j++ )
1333: if ( q = mat[i][j] ) {
1334: t = (unsigned int)rem(NM(q),md);
1335: if ( SGN(q) < 0 )
1336: t = (md - t) % md;
1337: wmat[i][j] = t;
1338: }
1339: }
1340: TOGFMMAT(row,col,wmat,*rp);
1341: }
1342:
1343: void Pgeninvm_swap(arg,rp)
1344: NODE arg;
1345: LIST *rp;
1346: {
1347: MAT m;
1348: pointer **mat;
1349: Q **tmat;
1350: Q *tvect;
1351: Q q;
1352: unsigned int **wmat,**invmat;
1353: int *index;
1354: unsigned int t,md;
1355: int i,j,row,col,status;
1356: MAT mat1;
1357: VECT vect1;
1358: NODE node1,node2;
1359:
1360: asir_assert(ARG0(arg),O_MAT,"geninvm_swap");
1361: asir_assert(ARG1(arg),O_N,"geninvm_swap");
1362: m = (MAT)ARG0(arg); md = QTOS((Q)ARG1(arg));
1363: row = m->row; col = m->col; mat = m->body;
1364: wmat = (unsigned int **)almat(row,col+row);
1365: for ( i = 0; i < row; i++ ) {
1366: bzero((char *)wmat[i],(col+row)*sizeof(int));
1367: for ( j = 0; j < col; j++ )
1368: if ( q = (Q)mat[i][j] ) {
1369: t = (unsigned int)rem(NM(q),md);
1370: if ( SGN(q) < 0 )
1371: t = (md - t) % md;
1372: wmat[i][j] = t;
1373: }
1374: wmat[i][col+i] = 1;
1375: }
1376: status = gauss_elim_geninv_mod_swap(wmat,row,col,md,&invmat,&index);
1377: if ( status > 0 )
1378: *rp = 0;
1379: else {
1380: MKMAT(mat1,col,col);
1381: for ( i = 0, tmat = (Q **)mat1->body; i < col; i++ )
1382: for ( j = 0; j < col; j++ )
1383: UTOQ(invmat[i][j],tmat[i][j]);
1384: MKVECT(vect1,row);
1385: for ( i = 0, tvect = (Q *)vect1->body; i < row; i++ )
1386: STOQ(index[i],tvect[i]);
1387: MKNODE(node2,vect1,0); MKNODE(node1,mat1,node2); MKLIST(*rp,node1);
1388: }
1389: }
1390:
1391: gauss_elim_geninv_mod_swap(mat,row,col,md,invmatp,indexp)
1392: unsigned int **mat;
1393: int row,col;
1394: unsigned int md;
1395: unsigned int ***invmatp;
1396: int **indexp;
1397: {
1398: int i,j,k,inv,a,n,m;
1399: unsigned int *t,*pivot,*s;
1400: int *index;
1401: unsigned int **invmat;
1402:
1403: n = col; m = row+col;
1404: *indexp = index = (int *)MALLOC_ATOMIC(row*sizeof(int));
1405: for ( i = 0; i < row; i++ )
1406: index[i] = i;
1407: for ( j = 0; j < n; j++ ) {
1408: for ( i = j; i < row && !mat[i][j]; i++ );
1409: if ( i == row ) {
1410: *indexp = 0; *invmatp = 0; return 1;
1411: }
1412: if ( i != j ) {
1413: t = mat[i]; mat[i] = mat[j]; mat[j] = t;
1414: k = index[i]; index[i] = index[j]; index[j] = k;
1415: }
1416: pivot = mat[j];
1417: inv = (unsigned int)invm(pivot[j],md);
1418: for ( k = j; k < m; k++ )
1419: if ( pivot[k] )
1420: pivot[k] = (unsigned int)dmar(pivot[k],inv,0,md);
1421: for ( i = j+1; i < row; i++ ) {
1422: t = mat[i];
1423: if ( a = t[j] )
1424: for ( k = j, a = md - a; k < m; k++ )
1425: if ( pivot[k] )
1426: t[k] = dmar(pivot[k],a,t[k],md);
1427: }
1428: }
1429: for ( j = n-1; j >= 0; j-- ) {
1430: pivot = mat[j];
1431: for ( i = j-1; i >= 0; i-- ) {
1432: t = mat[i];
1433: if ( a = t[j] )
1434: for ( k = j, a = md - a; k < m; k++ )
1435: if ( pivot[k] )
1436: t[k] = dmar(pivot[k],a,t[k],md);
1437: }
1438: }
1439: *invmatp = invmat = (unsigned int **)almat(col,col);
1440: for ( i = 0; i < col; i++ )
1441: for ( j = 0, s = invmat[i], t = mat[i]; j < col; j++ )
1442: s[j] = t[col+index[j]];
1443: return 0;
1444: }
1445:
1446: void _addn(N,N,N);
1447: int _subn(N,N,N);
1448: void _muln(N,N,N);
1449:
1450: void inner_product_int(a,b,n,r)
1451: Q *a,*b;
1452: int n;
1453: Q *r;
1454: {
1455: int la,lb,i;
1456: int sgn,sgn1;
1457: N wm,wma,sum,t;
1458:
1459: for ( la = lb = 0, i = 0; i < n; i++ ) {
1460: if ( a[i] )
1461: if ( DN(a[i]) )
1462: error("inner_product_int : invalid argument");
1463: else
1464: la = MAX(PL(NM(a[i])),la);
1465: if ( b[i] )
1466: if ( DN(b[i]) )
1467: error("inner_product_int : invalid argument");
1468: else
1469: lb = MAX(PL(NM(b[i])),lb);
1470: }
1471: sgn = 0;
1472: sum= NALLOC(la+lb+2);
1473: bzero((char *)sum,(la+lb+3)*sizeof(unsigned int));
1474: wm = NALLOC(la+lb+2);
1475: wma = NALLOC(la+lb+2);
1476: for ( i = 0; i < n; i++ ) {
1477: if ( !a[i] || !b[i] )
1478: continue;
1479: _muln(NM(a[i]),NM(b[i]),wm);
1480: sgn1 = SGN(a[i])*SGN(b[i]);
1481: if ( !sgn ) {
1482: sgn = sgn1;
1483: t = wm; wm = sum; sum = t;
1484: } else if ( sgn == sgn1 ) {
1485: _addn(sum,wm,wma);
1486: if ( !PL(wma) )
1487: sgn = 0;
1488: t = wma; wma = sum; sum = t;
1489: } else {
1490: /* sgn*sum+sgn1*wm = sgn*(sum-wm) */
1491: sgn *= _subn(sum,wm,wma);
1492: t = wma; wma = sum; sum = t;
1493: }
1494: }
1495: GC_free(wm);
1496: GC_free(wma);
1497: if ( !sgn ) {
1498: GC_free(sum);
1499: *r = 0;
1500: } else
1501: NTOQ(sum,sgn,*r);
1502: }
1503:
1504: void Pmul_mat_vect_int(arg,rp)
1505: NODE arg;
1506: VECT *rp;
1507: {
1508: MAT mat;
1509: VECT vect,r;
1510: int row,col,i;
1511:
1512: mat = (MAT)ARG0(arg);
1513: vect = (VECT)ARG1(arg);
1514: row = mat->row;
1515: col = mat->col;
1516: MKVECT(r,row);
1517: for ( i = 0; i < row; i++ )
1518: inner_product_int(mat->body[i],vect->body,col,&r->body[i]);
1519: *rp = r;
1520: }
1521:
1522: void Pnbpoly_up2(arg,rp)
1523: NODE arg;
1524: GF2N *rp;
1525: {
1526: int m,type,ret;
1527: UP2 r;
1528:
1529: m = QTOS((Q)ARG0(arg));
1530: type = QTOS((Q)ARG1(arg));
1531: ret = generate_ONB_polynomial(&r,m,type);
1532: if ( ret == 0 )
1533: MKGF2N(r,*rp);
1534: else
1535: *rp = 0;
1536: }
1537:
1538: void Px962_irredpoly_up2(arg,rp)
1539: NODE arg;
1540: GF2N *rp;
1541: {
1542: int m,type,ret,w;
1543: GF2N prev;
1544: UP2 r;
1545:
1546: m = QTOS((Q)ARG0(arg));
1547: prev = (GF2N)ARG1(arg);
1548: if ( !prev ) {
1549: w = (m>>5)+1; NEWUP2(r,w); r->w = 0;
1550: bzero((char *)r->b,w*sizeof(unsigned int));
1551: } else {
1552: r = prev->body;
1553: if ( degup2(r) != m ) {
1554: w = (m>>5)+1; NEWUP2(r,w); r->w = 0;
1555: bzero((char *)r->b,w*sizeof(unsigned int));
1556: }
1557: }
1558: ret = _generate_irreducible_polynomial(r,m,type);
1559: if ( ret == 0 )
1560: MKGF2N(r,*rp);
1561: else
1562: *rp = 0;
1563: }
1564:
1565: void Pirredpoly_up2(arg,rp)
1566: NODE arg;
1567: GF2N *rp;
1568: {
1569: int m,type,ret,w;
1570: GF2N prev;
1571: UP2 r;
1572:
1573: m = QTOS((Q)ARG0(arg));
1574: prev = (GF2N)ARG1(arg);
1575: if ( !prev ) {
1576: w = (m>>5)+1; NEWUP2(r,w); r->w = 0;
1577: bzero((char *)r->b,w*sizeof(unsigned int));
1578: } else {
1579: r = prev->body;
1580: if ( degup2(r) != m ) {
1581: w = (m>>5)+1; NEWUP2(r,w); r->w = 0;
1582: bzero((char *)r->b,w*sizeof(unsigned int));
1583: }
1584: }
1585: ret = _generate_good_irreducible_polynomial(r,m,type);
1586: if ( ret == 0 )
1587: MKGF2N(r,*rp);
1588: else
1589: *rp = 0;
1590: }
1591:
1592: /*
1593: * f = type 'type' normal polynomial of degree m if exists
1594: * IEEE P1363 A.7.2
1595: *
1596: * return value : 0 --- exists
1597: * 1 --- does not exist
1598: * -1 --- failure (memory allocation error)
1599: */
1600:
1601: int generate_ONB_polynomial(UP2 *rp,int m,int type)
1602: {
1603: int i,r;
1604: int w;
1605: UP2 f,f0,f1,f2,t;
1606:
1607: w = (m>>5)+1;
1608: switch ( type ) {
1609: case 1:
1610: if ( !TypeT_NB_check(m,1) ) return 1;
1611: NEWUP2(f,w); *rp = f; f->w = w;
1612: /* set all the bits */
1613: for ( i = 0; i < w; i++ )
1614: f->b[i] = 0xffffffff;
1615: /* mask the top word if necessary */
1616: if ( r = (m+1)&31 )
1617: f->b[w-1] &= (1<<r)-1;
1618: return 0;
1619: break;
1620: case 2:
1621: if ( !TypeT_NB_check(m,2) ) return 1;
1622: NEWUP2(f,w); *rp = f;
1623: W_NEWUP2(f0,w);
1624: W_NEWUP2(f1,w);
1625: W_NEWUP2(f2,w);
1626:
1627: /* recursion for genrating Type II normal polynomial */
1628:
1629: /* f0 = 1, f1 = t+1 */
1630: f0->w = 1; f0->b[0] = 1;
1631: f1->w = 1; f1->b[0] = 3;
1632: for ( i = 2; i <= m; i++ ) {
1633: /* f2 = t*f1+f0 */
1634: _bshiftup2(f1,-1,f2);
1635: _addup2_destructive(f2,f0);
1636: /* cyclic change of the variables */
1637: t = f0; f0 = f1; f1 = f2; f2 = t;
1638: }
1639: _copyup2(f1,f);
1640: return 0;
1641: break;
1642: default:
1643: return -1;
1644: break;
1645: }
1646: }
1647:
1648: /*
1649: * f = an irreducible trinomial or pentanomial of degree d 'after' f
1650: * return value : 0 --- exists
1651: * 1 --- does not exist (exhaustion)
1652: */
1653:
1654: int _generate_irreducible_polynomial(UP2 f,int d)
1655: {
1656: int ret,i,j,k,nz,i0,j0,k0;
1657: int w;
1658: unsigned int *fd;
1659:
1660: /*
1661: * if f = x^d+x^i+1 then i0 <- i, j0 <- 0, k0 <-0.
1662: * if f = x^d+x^k+x^j+x^i+1 (k>j>i) then i0 <- i, j0 <- j, k0 <-k.
1663: * otherwise i0,j0,k0 is set to 0.
1664: */
1665:
1666: fd = f->b;
1667: w = (d>>5)+1;
1668: if ( f->w && (d==degup2(f)) ) {
1669: for ( nz = 0, i = d; i >= 0; i-- )
1670: if ( fd[i>>5]&(1<<(i&31)) ) nz++;
1671: switch ( nz ) {
1672: case 3:
1673: for ( i0 = 1; !(fd[i0>>5]&(1<<(i0&31))) ; i0++ );
1674: /* reset i0-th bit */
1675: fd[i0>>5] &= ~(1<<(i0&31));
1676: j0 = k0 = 0;
1677: break;
1678: case 5:
1679: for ( i0 = 1; !(fd[i0>>5]&(1<<(i0&31))) ; i0++ );
1680: /* reset i0-th bit */
1681: fd[i0>>5] &= ~(1<<(i0&31));
1682: for ( j0 = i0+1; !(fd[j0>>5]&(1<<(j0&31))) ; j0++ );
1683: /* reset j0-th bit */
1684: fd[j0>>5] &= ~(1<<(j0&31));
1685: for ( k0 = j0+1; !(fd[k0>>5]&(1<<(k0&31))) ; k0++ );
1686: /* reset k0-th bit */
1687: fd[k0>>5] &= ~(1<<(k0&31));
1688: break;
1689: default:
1690: f->w = 0; break;
1691: }
1692: } else
1693: f->w = 0;
1694:
1695: if ( !f->w ) {
1696: fd = f->b;
1697: f->w = w; fd[0] |= 1; fd[d>>5] |= (1<<(d&31));
1698: i0 = j0 = k0 = 0;
1699: }
1700: /* if j0 > 0 then f is already a pentanomial */
1701: if ( j0 > 0 ) goto PENTA;
1702:
1703: /* searching for an irreducible trinomial */
1704:
1705: for ( i = 1; 2*i <= d; i++ ) {
1706: /* skip the polynomials 'before' f */
1707: if ( i < i0 ) continue;
1708: if ( i == i0 ) { i0 = 0; continue; }
1709: /* set i-th bit */
1710: fd[i>>5] |= (1<<(i&31));
1711: ret = irredcheck_dddup2(f);
1712: if ( ret == 1 ) return 0;
1713: /* reset i-th bit */
1714: fd[i>>5] &= ~(1<<(i&31));
1715: }
1716:
1717: /* searching for an irreducible pentanomial */
1718: PENTA:
1719: for ( i = 1; i < d; i++ ) {
1720: /* skip the polynomials 'before' f */
1721: if ( i < i0 ) continue;
1722: if ( i == i0 ) i0 = 0;
1723: /* set i-th bit */
1724: fd[i>>5] |= (1<<(i&31));
1725: for ( j = i+1; j < d; j++ ) {
1726: /* skip the polynomials 'before' f */
1727: if ( j < j0 ) continue;
1728: if ( j == j0 ) j0 = 0;
1729: /* set j-th bit */
1730: fd[j>>5] |= (1<<(j&31));
1731: for ( k = j+1; k < d; k++ ) {
1732: /* skip the polynomials 'before' f */
1733: if ( k < k0 ) continue;
1734: else if ( k == k0 ) { k0 = 0; continue; }
1735: /* set k-th bit */
1736: fd[k>>5] |= (1<<(k&31));
1737: ret = irredcheck_dddup2(f);
1738: if ( ret == 1 ) return 0;
1739: /* reset k-th bit */
1740: fd[k>>5] &= ~(1<<(k&31));
1741: }
1742: /* reset j-th bit */
1743: fd[j>>5] &= ~(1<<(j&31));
1744: }
1745: /* reset i-th bit */
1746: fd[i>>5] &= ~(1<<(i&31));
1747: }
1748: /* exhausted */
1749: return 1;
1750: }
1751:
1752: /*
1753: * f = an irreducible trinomial or pentanomial of degree d 'after' f
1754: *
1755: * searching strategy:
1756: * trinomial x^d+x^i+1:
1757: * i is as small as possible.
1758: * trinomial x^d+x^i+x^j+x^k+1:
1759: * i is as small as possible.
1760: * For such i, j is as small as possible.
1761: * For such i and j, 'k' is as small as possible.
1762: *
1763: * return value : 0 --- exists
1764: * 1 --- does not exist (exhaustion)
1765: */
1766:
1767: int _generate_good_irreducible_polynomial(UP2 f,int d)
1768: {
1769: int ret,i,j,k,nz,i0,j0,k0;
1770: int w;
1771: unsigned int *fd;
1772:
1773: /*
1774: * if f = x^d+x^i+1 then i0 <- i, j0 <- 0, k0 <-0.
1775: * if f = x^d+x^k+x^j+x^i+1 (k>j>i) then i0 <- i, j0 <- j, k0 <-k.
1776: * otherwise i0,j0,k0 is set to 0.
1777: */
1778:
1779: fd = f->b;
1780: w = (d>>5)+1;
1781: if ( f->w && (d==degup2(f)) ) {
1782: for ( nz = 0, i = d; i >= 0; i-- )
1783: if ( fd[i>>5]&(1<<(i&31)) ) nz++;
1784: switch ( nz ) {
1785: case 3:
1786: for ( i0 = 1; !(fd[i0>>5]&(1<<(i0&31))) ; i0++ );
1787: /* reset i0-th bit */
1788: fd[i0>>5] &= ~(1<<(i0&31));
1789: j0 = k0 = 0;
1790: break;
1791: case 5:
1792: for ( i0 = 1; !(fd[i0>>5]&(1<<(i0&31))) ; i0++ );
1793: /* reset i0-th bit */
1794: fd[i0>>5] &= ~(1<<(i0&31));
1795: for ( j0 = i0+1; !(fd[j0>>5]&(1<<(j0&31))) ; j0++ );
1796: /* reset j0-th bit */
1797: fd[j0>>5] &= ~(1<<(j0&31));
1798: for ( k0 = j0+1; !(fd[k0>>5]&(1<<(k0&31))) ; k0++ );
1799: /* reset k0-th bit */
1800: fd[k0>>5] &= ~(1<<(k0&31));
1801: break;
1802: default:
1803: f->w = 0; break;
1804: }
1805: } else
1806: f->w = 0;
1807:
1808: if ( !f->w ) {
1809: fd = f->b;
1810: f->w = w; fd[0] |= 1; fd[d>>5] |= (1<<(d&31));
1811: i0 = j0 = k0 = 0;
1812: }
1813: /* if j0 > 0 then f is already a pentanomial */
1814: if ( j0 > 0 ) goto PENTA;
1815:
1816: /* searching for an irreducible trinomial */
1817:
1818: for ( i = 1; 2*i <= d; i++ ) {
1819: /* skip the polynomials 'before' f */
1820: if ( i < i0 ) continue;
1821: if ( i == i0 ) { i0 = 0; continue; }
1822: /* set i-th bit */
1823: fd[i>>5] |= (1<<(i&31));
1824: ret = irredcheck_dddup2(f);
1825: if ( ret == 1 ) return 0;
1826: /* reset i-th bit */
1827: fd[i>>5] &= ~(1<<(i&31));
1828: }
1829:
1830: /* searching for an irreducible pentanomial */
1831: PENTA:
1832: for ( i = 3; i < d; i++ ) {
1833: /* skip the polynomials 'before' f */
1834: if ( i < i0 ) continue;
1835: if ( i == i0 ) i0 = 0;
1836: /* set i-th bit */
1837: fd[i>>5] |= (1<<(i&31));
1838: for ( j = 2; j < i; j++ ) {
1839: /* skip the polynomials 'before' f */
1840: if ( j < j0 ) continue;
1841: if ( j == j0 ) j0 = 0;
1842: /* set j-th bit */
1843: fd[j>>5] |= (1<<(j&31));
1844: for ( k = 1; k < j; k++ ) {
1845: /* skip the polynomials 'before' f */
1846: if ( k < k0 ) continue;
1847: else if ( k == k0 ) { k0 = 0; continue; }
1848: /* set k-th bit */
1849: fd[k>>5] |= (1<<(k&31));
1850: ret = irredcheck_dddup2(f);
1851: if ( ret == 1 ) return 0;
1852: /* reset k-th bit */
1853: fd[k>>5] &= ~(1<<(k&31));
1854: }
1855: /* reset j-th bit */
1856: fd[j>>5] &= ~(1<<(j&31));
1857: }
1858: /* reset i-th bit */
1859: fd[i>>5] &= ~(1<<(i&31));
1860: }
1861: /* exhausted */
1862: return 1;
1863: }