![[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.2 ! noro 1: /* $OpenXM: OpenXM_contrib2/asir2000/builtin/array.c,v 1.1.1.1 1999/12/03 07:39:07 noro Exp $ */
1.1 noro 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++ ) {
1.2 ! noro 709: if ( Print ) {
! 710: fprintf(asir_out,"."); fflush(asir_out);
! 711: }
1.1 noro 712: md = lprime[ind];
713: get_eg(&tmp0);
714: for ( i = 0; i < row; i++ )
715: for ( j = 0, bmi = bmat[i], wmi = wmat[i]; j < col; j++ )
716: if ( q = (Q)bmi[j] ) {
717: t = rem(NM(q),md);
718: if ( t && SGN(q) < 0 )
719: t = (md - t) % md;
720: wmi[j] = t;
721: } else
722: wmi[j] = 0;
723: get_eg(&tmp1);
724: add_eg(&eg_mod,&tmp0,&tmp1);
725: add_eg(&eg_mod_split,&tmp0,&tmp1);
726: get_eg(&tmp0);
727: rank = generic_gauss_elim_mod(wmat,row,col,md,wcolstat);
728: get_eg(&tmp1);
729: add_eg(&eg_elim,&tmp0,&tmp1);
730: add_eg(&eg_elim_split,&tmp0,&tmp1);
731: if ( !ind ) {
732: RESET:
733: UTON(md,m1);
734: rank0 = rank;
735: bcopy(wcolstat,colstat,col*sizeof(int));
736: MKMAT(crmat,rank,col-rank);
737: MKMAT(r,rank,col-rank); *nm = r;
738: tmat = (N **)crmat->body;
739: for ( i = 0; i < rank; i++ )
740: for ( j = k = 0, tmi = tmat[i], wmi = wmat[i]; j < col; j++ )
741: if ( !colstat[j] ) {
742: UTON(wmi[j],tmi[k]); k++;
743: }
744: } else {
745: if ( rank < rank0 ) {
1.2 ! noro 746: if ( Print ) {
1.1 noro 747: fprintf(asir_out,"lower rank matrix; continuing...\n");
1.2 ! noro 748: fflush(asir_out);
! 749: }
1.1 noro 750: continue;
751: } else if ( rank > rank0 ) {
1.2 ! noro 752: if ( Print ) {
1.1 noro 753: fprintf(asir_out,"higher rank matrix; resetting...\n");
1.2 ! noro 754: fflush(asir_out);
! 755: }
1.1 noro 756: goto RESET;
757: } else {
758: for ( j = 0; (j<col) && (colstat[j]==wcolstat[j]); j++ );
759: if ( j < col ) {
1.2 ! noro 760: if ( Print ) {
1.1 noro 761: fprintf(asir_out,"inconsitent colstat; resetting...\n");
1.2 ! noro 762: fflush(asir_out);
! 763: }
1.1 noro 764: goto RESET;
765: }
766: }
767:
768: get_eg(&tmp0);
769: inv = invm(rem(m1,md),md);
770: UTON(md,m2); muln(m1,m2,&m3);
771: for ( i = 0; i < rank; i++ )
772: for ( j = k = 0, tmi = tmat[i], wmi = wmat[i]; j < col; j++ )
773: if ( !colstat[j] ) {
774: if ( tmi[k] ) {
775: /* f3 = f1+m1*(m1 mod m2)^(-1)*(f2 - f1 mod m2) */
776: t = rem(tmi[k],md);
777: if ( wmi[j] >= t )
778: t = wmi[j]-t;
779: else
780: t = md-(t-wmi[j]);
781: DMAR(t,inv,0,md,t1)
782: UTON(t1,u);
783: muln(m1,u,&s);
784: addn(tmi[k],s,&u); tmi[k] = u;
785: } else if ( wmi[j] ) {
786: /* f3 = m1*(m1 mod m2)^(-1)*f2 */
787: DMAR(wmi[j],inv,0,md,t)
788: UTON(t,u);
789: muln(m1,u,&s); tmi[k] = s;
790: }
791: k++;
792: }
793: m1 = m3;
794: get_eg(&tmp1);
795: add_eg(&eg_chrem,&tmp0,&tmp1);
796: add_eg(&eg_chrem_split,&tmp0,&tmp1);
797:
798: get_eg(&tmp0);
799: ret = intmtoratm(crmat,m1,*nm,dn);
800: get_eg(&tmp1);
801: add_eg(&eg_intrat,&tmp0,&tmp1);
802: add_eg(&eg_intrat_split,&tmp0,&tmp1);
803: if ( ret ) {
804: *rindp = rind = (int *)MALLOC_ATOMIC(rank*sizeof(int));
805: *cindp = cind = (int *)MALLOC_ATOMIC((col-rank)*sizeof(int));
806: for ( j = k = l = 0; j < col; j++ )
807: if ( colstat[j] )
808: rind[k++] = j;
809: else
810: cind[l++] = j;
811: get_eg(&tmp0);
812: if ( gensolve_check(mat,*nm,*dn,rind,cind) )
813: get_eg(&tmp1);
814: add_eg(&eg_gschk,&tmp0,&tmp1);
815: add_eg(&eg_gschk_split,&tmp0,&tmp1);
816: if ( Print ) {
817: print_eg("Mod",&eg_mod_split);
818: print_eg("Elim",&eg_elim_split);
819: print_eg("ChRem",&eg_chrem_split);
820: print_eg("IntRat",&eg_intrat_split);
821: print_eg("Check",&eg_gschk_split);
822: fflush(asir_out);
823: }
824: return rank;
825: }
826: }
827: }
828: }
829:
830: int f4_nocheck;
831:
832: int gensolve_check(mat,nm,dn,rind,cind)
833: MAT mat,nm;
834: Q dn;
835: int *rind,*cind;
836: {
837: int row,col,rank,clen,i,j,k,l;
838: Q s,t,u;
839: Q *w;
840: Q *mati,*nmk;
841:
842: if ( f4_nocheck )
843: return 1;
844: row = mat->row; col = mat->col;
845: rank = nm->row; clen = nm->col;
846: w = (Q *)MALLOC(clen*sizeof(Q));
847: for ( i = 0; i < row; i++ ) {
848: mati = (Q *)mat->body[i];
849: #if 1
850: bzero(w,clen*sizeof(Q));
851: for ( k = 0; k < rank; k++ )
852: for ( l = 0, nmk = (Q *)nm->body[k]; l < clen; l++ ) {
853: mulq(mati[rind[k]],nmk[l],&t);
854: addq(w[l],t,&s); w[l] = s;
855: }
856: for ( j = 0; j < clen; j++ ) {
857: mulq(dn,mati[cind[j]],&t);
858: if ( cmpq(w[j],t) )
859: break;
860: }
861: #else
862: for ( j = 0; j < clen; j++ ) {
863: for ( k = 0, s = 0; k < rank; k++ ) {
864: mulq(mati[rind[k]],nm->body[k][j],&t);
865: addq(s,t,&u); s = u;
866: }
867: mulq(dn,mati[cind[j]],&t);
868: if ( cmpq(s,t) )
869: break;
870: }
871: #endif
872: if ( j != clen )
873: break;
874: }
875: if ( i != row )
876: return 0;
877: else
878: return 1;
879: }
880:
881: /* assuming 0 < c < m */
882:
883: int inttorat(c,m,b,sgnp,nmp,dnp)
884: N c,m,b;
885: int *sgnp;
886: N *nmp,*dnp;
887: {
888: Q qq,t,u1,v1,r1,nm;
889: N q,r,u2,v2,r2;
890:
891: u1 = 0; v1 = ONE; u2 = m; v2 = c;
892: while ( cmpn(v2,b) >= 0 ) {
893: divn(u2,v2,&q,&r2); u2 = v2; v2 = r2;
894: NTOQ(q,1,qq); mulq(qq,v1,&t); subq(u1,t,&r1); u1 = v1; v1 = r1;
895: }
896: if ( cmpn(NM(v1),b) >= 0 )
897: return 0;
898: else {
899: *nmp = v2;
900: *dnp = NM(v1);
901: *sgnp = SGN(v1);
902: return 1;
903: }
904: }
905:
906: /* mat->body = N ** */
907:
908: int intmtoratm(mat,md,nm,dn)
909: MAT mat;
910: N md;
911: MAT nm;
912: Q *dn;
913: {
914: N t,s,b;
915: Q bound,dn0,dn1,nm1,q,tq;
916: int i,j,k,l,row,col;
917: Q **rmat;
918: N **tmat;
919: N *tmi;
920: Q *nmk;
921: N u,unm,udn;
922: int sgn,ret;
923:
924: row = mat->row; col = mat->col;
925: bshiftn(md,1,&t);
926: isqrt(t,&s);
927: bshiftn(s,64,&b);
928: if ( !b )
929: b = ONEN;
930: dn0 = ONE;
931: tmat = (N **)mat->body;
932: rmat = (Q **)nm->body;
933: for ( i = 0; i < row; i++ )
934: for ( j = 0, tmi = tmat[i]; j < col; j++ )
935: if ( tmi[j] ) {
936: muln(tmi[j],NM(dn0),&s);
937: remn(s,md,&u);
938: ret = inttorat(u,md,b,&sgn,&unm,&udn);
939: if ( !ret )
940: return 0;
941: else {
942: NTOQ(unm,sgn,nm1);
943: NTOQ(udn,1,dn1);
944: if ( !UNIQ(dn1) ) {
945: for ( k = 0; k < i; k++ )
946: for ( l = 0, nmk = rmat[k]; l < col; l++ ) {
947: mulq(nmk[l],dn1,&q); nmk[l] = q;
948: }
949: for ( l = 0, nmk = rmat[i]; l < j; l++ ) {
950: mulq(nmk[l],dn1,&q); nmk[l] = q;
951: }
952: }
953: rmat[i][j] = nm1;
954: mulq(dn0,dn1,&q); dn0 = q;
955: }
956: }
957: *dn = dn0;
958: return 1;
959: }
960:
961: int generic_gauss_elim_mod(mat,row,col,md,colstat)
962: int **mat;
963: int row,col,md;
964: int *colstat;
965: {
966: int i,j,k,l,inv,a,rank;
967: int *t,*pivot;
968:
969: for ( rank = 0, j = 0; j < col; j++ ) {
970: for ( i = rank; i < row && !mat[i][j]; i++ );
971: if ( i == row ) {
972: colstat[j] = 0;
973: continue;
974: } else
975: colstat[j] = 1;
976: if ( i != rank ) {
977: t = mat[i]; mat[i] = mat[rank]; mat[rank] = t;
978: }
979: pivot = mat[rank];
980: inv = invm(pivot[j],md);
981: for ( k = j; k < col; k++ )
982: if ( pivot[k] ) {
983: DMAR(pivot[k],inv,0,md,pivot[k])
984: }
985: for ( i = rank+1; i < row; i++ ) {
986: t = mat[i];
987: if ( a = t[j] )
988: for ( k = j, a = md - a; k < col; k++ )
989: if ( pivot[k] ) {
990: DMAR(pivot[k],a,t[k],md,t[k])
991: }
992: }
993: rank++;
994: }
995: for ( j = col-1, l = rank-1; j >= 0; j-- )
996: if ( colstat[j] ) {
997: pivot = mat[l];
998: for ( i = 0; i < l; i++ ) {
999: t = mat[i];
1000: if ( a = t[j] )
1001: for ( k = j, a = md-a; k < col; k++ )
1002: if ( pivot[k] ) {
1003: DMAR(pivot[k],a,t[k],md,t[k])
1004: }
1005: }
1006: l--;
1007: }
1008: return rank;
1009: }
1010:
1011: /* LU decomposition; a[i][i] = 1/U[i][i] */
1012:
1013: int lu_gfmmat(mat,md,perm)
1014: GFMMAT mat;
1015: unsigned int md;
1016: int *perm;
1017: {
1018: int row,col;
1019: int i,j,k,l;
1020: unsigned int *t,*pivot;
1021: unsigned int **a;
1022: unsigned int inv,m;
1023:
1024: row = mat->row; col = mat->col;
1025: a = mat->body;
1026: bzero(perm,row*sizeof(int));
1027:
1028: for ( i = 0; i < row; i++ )
1029: perm[i] = i;
1030: for ( k = 0; k < col; k++ ) {
1031: for ( i = k; i < row && !a[i][k]; i++ );
1032: if ( i == row )
1033: return 0;
1034: if ( i != k ) {
1035: j = perm[i]; perm[i] = perm[k]; perm[k] = j;
1036: t = a[i]; a[i] = a[k]; a[k] = t;
1037: }
1038: pivot = a[k];
1039: pivot[k] = inv = invm(pivot[k],md);
1040: for ( i = k+1; i < row; i++ ) {
1041: t = a[i];
1042: if ( m = t[k] ) {
1043: DMAR(inv,m,0,md,t[k])
1044: for ( j = k+1, m = md - t[k]; j < col; j++ )
1045: if ( pivot[j] ) {
1046: DMAR(m,pivot[j],t[j],md,t[j])
1047: }
1048: }
1049: }
1050: }
1051: return 1;
1052: }
1053:
1054: void Pleqm1(arg,rp)
1055: NODE arg;
1056: VECT *rp;
1057: {
1058: MAT m;
1059: VECT vect;
1060: pointer **mat;
1061: Q *v;
1062: Q q;
1063: int **wmat;
1064: int md,i,j,row,col,t,n,status;
1065:
1066: asir_assert(ARG0(arg),O_MAT,"leqm1");
1067: asir_assert(ARG1(arg),O_N,"leqm1");
1068: m = (MAT)ARG0(arg); md = QTOS((Q)ARG1(arg));
1069: row = m->row; col = m->col; mat = m->body;
1070: wmat = (int **)almat(row,col);
1071: for ( i = 0; i < row; i++ )
1072: for ( j = 0; j < col; j++ )
1073: if ( q = (Q)mat[i][j] ) {
1074: t = rem(NM(q),md);
1075: if ( SGN(q) < 0 )
1076: t = (md - t) % md;
1077: wmat[i][j] = t;
1078: } else
1079: wmat[i][j] = 0;
1080: status = gauss_elim_mod1(wmat,row,col,md);
1081: if ( status < 0 )
1082: *rp = 0;
1083: else if ( status > 0 )
1084: *rp = (VECT)ONE;
1085: else {
1086: n = col - 1;
1087: MKVECT(vect,n);
1088: for ( i = 0, v = (Q *)vect->body; i < n; i++ ) {
1089: t = (md-wmat[i][n])%md; STOQ(t,v[i]);
1090: }
1091: *rp = vect;
1092: }
1093: }
1094:
1095: gauss_elim_mod1(mat,row,col,md)
1096: int **mat;
1097: int row,col,md;
1098: {
1099: int i,j,k,inv,a,n;
1100: int *t,*pivot;
1101:
1102: n = col - 1;
1103: for ( j = 0; j < n; j++ ) {
1104: for ( i = j; i < row && !mat[i][j]; i++ );
1105: if ( i == row )
1106: return 1;
1107: if ( i != j ) {
1108: t = mat[i]; mat[i] = mat[j]; mat[j] = t;
1109: }
1110: pivot = mat[j];
1111: inv = invm(pivot[j],md);
1112: for ( k = j; k <= n; k++ )
1113: pivot[k] = dmar(pivot[k],inv,0,md);
1114: for ( i = j+1; i < row; i++ ) {
1115: t = mat[i];
1116: if ( i != j && (a = t[j]) )
1117: for ( k = j, a = md - a; k <= n; k++ )
1118: t[k] = dmar(pivot[k],a,t[k],md);
1119: }
1120: }
1121: for ( i = n; i < row && !mat[i][n]; i++ );
1122: if ( i == row ) {
1123: for ( j = n-1; j >= 0; j-- ) {
1124: for ( i = j-1, a = (md-mat[j][n])%md; i >= 0; i-- ) {
1125: mat[i][n] = dmar(mat[i][j],a,mat[i][n],md);
1126: mat[i][j] = 0;
1127: }
1128: }
1129: return 0;
1130: } else
1131: return -1;
1132: }
1133:
1134: void Pgeninvm(arg,rp)
1135: NODE arg;
1136: LIST *rp;
1137: {
1138: MAT m;
1139: pointer **mat;
1140: Q **tmat;
1141: Q q;
1142: unsigned int **wmat;
1143: int md,i,j,row,col,t,status;
1144: MAT mat1,mat2;
1145: NODE node1,node2;
1146:
1147: asir_assert(ARG0(arg),O_MAT,"leqm1");
1148: asir_assert(ARG1(arg),O_N,"leqm1");
1149: m = (MAT)ARG0(arg); md = QTOS((Q)ARG1(arg));
1150: row = m->row; col = m->col; mat = m->body;
1151: wmat = (unsigned int **)almat(row,col+row);
1152: for ( i = 0; i < row; i++ ) {
1153: bzero((char *)wmat[i],(col+row)*sizeof(int));
1154: for ( j = 0; j < col; j++ )
1155: if ( q = (Q)mat[i][j] ) {
1156: t = rem(NM(q),md);
1157: if ( SGN(q) < 0 )
1158: t = (md - t) % md;
1159: wmat[i][j] = t;
1160: }
1161: wmat[i][col+i] = 1;
1162: }
1163: status = gauss_elim_geninv_mod(wmat,row,col,md);
1164: if ( status > 0 )
1165: *rp = 0;
1166: else {
1167: MKMAT(mat1,col,row); MKMAT(mat2,row-col,row);
1168: for ( i = 0, tmat = (Q **)mat1->body; i < col; i++ )
1169: for ( j = 0; j < row; j++ )
1170: STOQ(wmat[i][j+col],tmat[i][j]);
1171: for ( tmat = (Q **)mat2->body; i < row; i++ )
1172: for ( j = 0; j < row; j++ )
1173: STOQ(wmat[i][j+col],tmat[i-col][j]);
1174: MKNODE(node2,mat2,0); MKNODE(node1,mat1,node2); MKLIST(*rp,node1);
1175: }
1176: }
1177:
1178: int gauss_elim_geninv_mod(mat,row,col,md)
1179: unsigned int **mat;
1180: int row,col,md;
1181: {
1182: int i,j,k,inv,a,n,m;
1183: unsigned int *t,*pivot;
1184:
1185: n = col; m = row+col;
1186: for ( j = 0; j < n; j++ ) {
1187: for ( i = j; i < row && !mat[i][j]; i++ );
1188: if ( i == row )
1189: return 1;
1190: if ( i != j ) {
1191: t = mat[i]; mat[i] = mat[j]; mat[j] = t;
1192: }
1193: pivot = mat[j];
1194: inv = invm(pivot[j],md);
1195: for ( k = j; k < m; k++ )
1196: pivot[k] = dmar(pivot[k],inv,0,md);
1197: for ( i = j+1; i < row; i++ ) {
1198: t = mat[i];
1199: if ( a = t[j] )
1200: for ( k = j, a = md - a; k < m; k++ )
1201: t[k] = dmar(pivot[k],a,t[k],md);
1202: }
1203: }
1204: for ( j = n-1; j >= 0; j-- ) {
1205: pivot = mat[j];
1206: for ( i = j-1; i >= 0; i-- ) {
1207: t = mat[i];
1208: if ( a = t[j] )
1209: for ( k = j, a = md - a; k < m; k++ )
1210: t[k] = dmar(pivot[k],a,t[k],md);
1211: }
1212: }
1213: return 0;
1214: }
1215:
1216: void Psolve_by_lu_gfmmat(arg,rp)
1217: NODE arg;
1218: VECT *rp;
1219: {
1220: GFMMAT lu;
1221: Q *perm,*rhs,*v;
1222: int n,i;
1223: unsigned int md;
1224: unsigned int *b,*sol;
1225: VECT r;
1226:
1227: lu = (GFMMAT)ARG0(arg);
1228: perm = (Q *)BDY((VECT)ARG1(arg));
1229: rhs = (Q *)BDY((VECT)ARG2(arg));
1230: md = (unsigned int)QTOS((Q)ARG3(arg));
1231: n = lu->col;
1232: b = (unsigned int *)MALLOC_ATOMIC(n*sizeof(int));
1233: sol = (unsigned int *)MALLOC_ATOMIC(n*sizeof(int));
1234: for ( i = 0; i < n; i++ )
1235: b[i] = QTOS(rhs[QTOS(perm[i])]);
1236: solve_by_lu_gfmmat(lu,md,b,sol);
1237: MKVECT(r,n);
1238: for ( i = 0, v = (Q *)r->body; i < n; i++ )
1239: STOQ(sol[i],v[i]);
1240: *rp = r;
1241: }
1242:
1243: void solve_by_lu_gfmmat(lu,md,b,x)
1244: GFMMAT lu;
1245: unsigned int md;
1246: unsigned int *b;
1247: unsigned int *x;
1248: {
1249: int n;
1250: unsigned int **a;
1251: unsigned int *y;
1252: int i,j;
1253: unsigned int t,m;
1254:
1255: n = lu->col;
1256: a = lu->body;
1257: y = (unsigned int *)MALLOC_ATOMIC(n*sizeof(int));
1258: /* solve Ly=b */
1259: for ( i = 0; i < n; i++ ) {
1260: for ( t = b[i], j = 0; j < i; j++ )
1261: if ( a[i][j] ) {
1262: m = md - a[i][j];
1263: DMAR(m,y[j],t,md,t)
1264: }
1265: y[i] = t;
1266: }
1267: /* solve Ux=y */
1268: for ( i = n-1; i >= 0; i-- ) {
1269: for ( t = y[i], j =i+1; j < n; j++ )
1270: if ( a[i][j] ) {
1271: m = md - a[i][j];
1272: DMAR(m,x[j],t,md,t)
1273: }
1274: /* a[i][i] = 1/U[i][i] */
1275: DMAR(t,a[i][i],0,md,x[i])
1276: }
1277: }
1278:
1279: void Plu_gfmmat(arg,rp)
1280: NODE arg;
1281: LIST *rp;
1282: {
1283: MAT m;
1284: GFMMAT mm;
1285: unsigned int md;
1286: int i,row,col,status;
1287: int *iperm;
1288: Q *v;
1289: VECT perm;
1290: NODE n0;
1291:
1292: asir_assert(ARG0(arg),O_MAT,"mat_to_gfmmat");
1293: asir_assert(ARG1(arg),O_N,"mat_to_gfmmat");
1294: m = (MAT)ARG0(arg); md = (unsigned int)QTOS((Q)ARG1(arg));
1295: mat_to_gfmmat(m,md,&mm);
1296: row = m->row;
1297: col = m->col;
1298: iperm = (int *)MALLOC_ATOMIC(row*sizeof(int));
1299: status = lu_gfmmat(mm,md,iperm);
1300: if ( !status )
1301: n0 = 0;
1302: else {
1303: MKVECT(perm,row);
1304: for ( i = 0, v = (Q *)perm->body; i < row; i++ )
1305: STOQ(iperm[i],v[i]);
1306: n0 = mknode(2,mm,perm);
1307: }
1308: MKLIST(*rp,n0);
1309: }
1310:
1311: void Pmat_to_gfmmat(arg,rp)
1312: NODE arg;
1313: GFMMAT *rp;
1314: {
1315: MAT m;
1316: unsigned int md;
1317:
1318: asir_assert(ARG0(arg),O_MAT,"mat_to_gfmmat");
1319: asir_assert(ARG1(arg),O_N,"mat_to_gfmmat");
1320: m = (MAT)ARG0(arg); md = (unsigned int)QTOS((Q)ARG1(arg));
1321: mat_to_gfmmat(m,md,rp);
1322: }
1323:
1324: void mat_to_gfmmat(m,md,rp)
1325: MAT m;
1326: unsigned int md;
1327: GFMMAT *rp;
1328: {
1329: unsigned int **wmat;
1330: unsigned int t;
1331: Q **mat;
1332: Q q;
1333: int i,j,row,col;
1334:
1335: row = m->row; col = m->col; mat = (Q **)m->body;
1336: wmat = (unsigned int **)almat(row,col);
1337: for ( i = 0; i < row; i++ ) {
1338: bzero((char *)wmat[i],col*sizeof(unsigned int));
1339: for ( j = 0; j < col; j++ )
1340: if ( q = mat[i][j] ) {
1341: t = (unsigned int)rem(NM(q),md);
1342: if ( SGN(q) < 0 )
1343: t = (md - t) % md;
1344: wmat[i][j] = t;
1345: }
1346: }
1347: TOGFMMAT(row,col,wmat,*rp);
1348: }
1349:
1350: void Pgeninvm_swap(arg,rp)
1351: NODE arg;
1352: LIST *rp;
1353: {
1354: MAT m;
1355: pointer **mat;
1356: Q **tmat;
1357: Q *tvect;
1358: Q q;
1359: unsigned int **wmat,**invmat;
1360: int *index;
1361: unsigned int t,md;
1362: int i,j,row,col,status;
1363: MAT mat1;
1364: VECT vect1;
1365: NODE node1,node2;
1366:
1367: asir_assert(ARG0(arg),O_MAT,"geninvm_swap");
1368: asir_assert(ARG1(arg),O_N,"geninvm_swap");
1369: m = (MAT)ARG0(arg); md = QTOS((Q)ARG1(arg));
1370: row = m->row; col = m->col; mat = m->body;
1371: wmat = (unsigned int **)almat(row,col+row);
1372: for ( i = 0; i < row; i++ ) {
1373: bzero((char *)wmat[i],(col+row)*sizeof(int));
1374: for ( j = 0; j < col; j++ )
1375: if ( q = (Q)mat[i][j] ) {
1376: t = (unsigned int)rem(NM(q),md);
1377: if ( SGN(q) < 0 )
1378: t = (md - t) % md;
1379: wmat[i][j] = t;
1380: }
1381: wmat[i][col+i] = 1;
1382: }
1383: status = gauss_elim_geninv_mod_swap(wmat,row,col,md,&invmat,&index);
1384: if ( status > 0 )
1385: *rp = 0;
1386: else {
1387: MKMAT(mat1,col,col);
1388: for ( i = 0, tmat = (Q **)mat1->body; i < col; i++ )
1389: for ( j = 0; j < col; j++ )
1390: UTOQ(invmat[i][j],tmat[i][j]);
1391: MKVECT(vect1,row);
1392: for ( i = 0, tvect = (Q *)vect1->body; i < row; i++ )
1393: STOQ(index[i],tvect[i]);
1394: MKNODE(node2,vect1,0); MKNODE(node1,mat1,node2); MKLIST(*rp,node1);
1395: }
1396: }
1397:
1398: gauss_elim_geninv_mod_swap(mat,row,col,md,invmatp,indexp)
1399: unsigned int **mat;
1400: int row,col;
1401: unsigned int md;
1402: unsigned int ***invmatp;
1403: int **indexp;
1404: {
1405: int i,j,k,inv,a,n,m;
1406: unsigned int *t,*pivot,*s;
1407: int *index;
1408: unsigned int **invmat;
1409:
1410: n = col; m = row+col;
1411: *indexp = index = (int *)MALLOC_ATOMIC(row*sizeof(int));
1412: for ( i = 0; i < row; i++ )
1413: index[i] = i;
1414: for ( j = 0; j < n; j++ ) {
1415: for ( i = j; i < row && !mat[i][j]; i++ );
1416: if ( i == row ) {
1417: *indexp = 0; *invmatp = 0; return 1;
1418: }
1419: if ( i != j ) {
1420: t = mat[i]; mat[i] = mat[j]; mat[j] = t;
1421: k = index[i]; index[i] = index[j]; index[j] = k;
1422: }
1423: pivot = mat[j];
1424: inv = (unsigned int)invm(pivot[j],md);
1425: for ( k = j; k < m; k++ )
1426: if ( pivot[k] )
1427: pivot[k] = (unsigned int)dmar(pivot[k],inv,0,md);
1428: for ( i = j+1; i < row; i++ ) {
1429: t = mat[i];
1430: if ( a = t[j] )
1431: for ( k = j, a = md - a; k < m; k++ )
1432: if ( pivot[k] )
1433: t[k] = dmar(pivot[k],a,t[k],md);
1434: }
1435: }
1436: for ( j = n-1; j >= 0; j-- ) {
1437: pivot = mat[j];
1438: for ( i = j-1; i >= 0; i-- ) {
1439: t = mat[i];
1440: if ( a = t[j] )
1441: for ( k = j, a = md - a; k < m; k++ )
1442: if ( pivot[k] )
1443: t[k] = dmar(pivot[k],a,t[k],md);
1444: }
1445: }
1446: *invmatp = invmat = (unsigned int **)almat(col,col);
1447: for ( i = 0; i < col; i++ )
1448: for ( j = 0, s = invmat[i], t = mat[i]; j < col; j++ )
1449: s[j] = t[col+index[j]];
1450: return 0;
1451: }
1452:
1453: void _addn(N,N,N);
1454: int _subn(N,N,N);
1455: void _muln(N,N,N);
1456:
1457: void inner_product_int(a,b,n,r)
1458: Q *a,*b;
1459: int n;
1460: Q *r;
1461: {
1462: int la,lb,i;
1463: int sgn,sgn1;
1464: N wm,wma,sum,t;
1465:
1466: for ( la = lb = 0, i = 0; i < n; i++ ) {
1467: if ( a[i] )
1468: if ( DN(a[i]) )
1469: error("inner_product_int : invalid argument");
1470: else
1471: la = MAX(PL(NM(a[i])),la);
1472: if ( b[i] )
1473: if ( DN(b[i]) )
1474: error("inner_product_int : invalid argument");
1475: else
1476: lb = MAX(PL(NM(b[i])),lb);
1477: }
1478: sgn = 0;
1479: sum= NALLOC(la+lb+2);
1480: bzero((char *)sum,(la+lb+3)*sizeof(unsigned int));
1481: wm = NALLOC(la+lb+2);
1482: wma = NALLOC(la+lb+2);
1483: for ( i = 0; i < n; i++ ) {
1484: if ( !a[i] || !b[i] )
1485: continue;
1486: _muln(NM(a[i]),NM(b[i]),wm);
1487: sgn1 = SGN(a[i])*SGN(b[i]);
1488: if ( !sgn ) {
1489: sgn = sgn1;
1490: t = wm; wm = sum; sum = t;
1491: } else if ( sgn == sgn1 ) {
1492: _addn(sum,wm,wma);
1493: if ( !PL(wma) )
1494: sgn = 0;
1495: t = wma; wma = sum; sum = t;
1496: } else {
1497: /* sgn*sum+sgn1*wm = sgn*(sum-wm) */
1498: sgn *= _subn(sum,wm,wma);
1499: t = wma; wma = sum; sum = t;
1500: }
1501: }
1502: GC_free(wm);
1503: GC_free(wma);
1504: if ( !sgn ) {
1505: GC_free(sum);
1506: *r = 0;
1507: } else
1508: NTOQ(sum,sgn,*r);
1509: }
1510:
1511: void Pmul_mat_vect_int(arg,rp)
1512: NODE arg;
1513: VECT *rp;
1514: {
1515: MAT mat;
1516: VECT vect,r;
1517: int row,col,i;
1518:
1519: mat = (MAT)ARG0(arg);
1520: vect = (VECT)ARG1(arg);
1521: row = mat->row;
1522: col = mat->col;
1523: MKVECT(r,row);
1524: for ( i = 0; i < row; i++ )
1525: inner_product_int(mat->body[i],vect->body,col,&r->body[i]);
1526: *rp = r;
1527: }
1528:
1529: void Pnbpoly_up2(arg,rp)
1530: NODE arg;
1531: GF2N *rp;
1532: {
1533: int m,type,ret;
1534: UP2 r;
1535:
1536: m = QTOS((Q)ARG0(arg));
1537: type = QTOS((Q)ARG1(arg));
1538: ret = generate_ONB_polynomial(&r,m,type);
1539: if ( ret == 0 )
1540: MKGF2N(r,*rp);
1541: else
1542: *rp = 0;
1543: }
1544:
1545: void Px962_irredpoly_up2(arg,rp)
1546: NODE arg;
1547: GF2N *rp;
1548: {
1549: int m,type,ret,w;
1550: GF2N prev;
1551: UP2 r;
1552:
1553: m = QTOS((Q)ARG0(arg));
1554: prev = (GF2N)ARG1(arg);
1555: if ( !prev ) {
1556: w = (m>>5)+1; NEWUP2(r,w); r->w = 0;
1557: bzero((char *)r->b,w*sizeof(unsigned int));
1558: } else {
1559: r = prev->body;
1560: if ( degup2(r) != m ) {
1561: w = (m>>5)+1; NEWUP2(r,w); r->w = 0;
1562: bzero((char *)r->b,w*sizeof(unsigned int));
1563: }
1564: }
1565: ret = _generate_irreducible_polynomial(r,m,type);
1566: if ( ret == 0 )
1567: MKGF2N(r,*rp);
1568: else
1569: *rp = 0;
1570: }
1571:
1572: void Pirredpoly_up2(arg,rp)
1573: NODE arg;
1574: GF2N *rp;
1575: {
1576: int m,type,ret,w;
1577: GF2N prev;
1578: UP2 r;
1579:
1580: m = QTOS((Q)ARG0(arg));
1581: prev = (GF2N)ARG1(arg);
1582: if ( !prev ) {
1583: w = (m>>5)+1; NEWUP2(r,w); r->w = 0;
1584: bzero((char *)r->b,w*sizeof(unsigned int));
1585: } else {
1586: r = prev->body;
1587: if ( degup2(r) != m ) {
1588: w = (m>>5)+1; NEWUP2(r,w); r->w = 0;
1589: bzero((char *)r->b,w*sizeof(unsigned int));
1590: }
1591: }
1592: ret = _generate_good_irreducible_polynomial(r,m,type);
1593: if ( ret == 0 )
1594: MKGF2N(r,*rp);
1595: else
1596: *rp = 0;
1597: }
1598:
1599: /*
1600: * f = type 'type' normal polynomial of degree m if exists
1601: * IEEE P1363 A.7.2
1602: *
1603: * return value : 0 --- exists
1604: * 1 --- does not exist
1605: * -1 --- failure (memory allocation error)
1606: */
1607:
1608: int generate_ONB_polynomial(UP2 *rp,int m,int type)
1609: {
1610: int i,r;
1611: int w;
1612: UP2 f,f0,f1,f2,t;
1613:
1614: w = (m>>5)+1;
1615: switch ( type ) {
1616: case 1:
1617: if ( !TypeT_NB_check(m,1) ) return 1;
1618: NEWUP2(f,w); *rp = f; f->w = w;
1619: /* set all the bits */
1620: for ( i = 0; i < w; i++ )
1621: f->b[i] = 0xffffffff;
1622: /* mask the top word if necessary */
1623: if ( r = (m+1)&31 )
1624: f->b[w-1] &= (1<<r)-1;
1625: return 0;
1626: break;
1627: case 2:
1628: if ( !TypeT_NB_check(m,2) ) return 1;
1629: NEWUP2(f,w); *rp = f;
1630: W_NEWUP2(f0,w);
1631: W_NEWUP2(f1,w);
1632: W_NEWUP2(f2,w);
1633:
1634: /* recursion for genrating Type II normal polynomial */
1635:
1636: /* f0 = 1, f1 = t+1 */
1637: f0->w = 1; f0->b[0] = 1;
1638: f1->w = 1; f1->b[0] = 3;
1639: for ( i = 2; i <= m; i++ ) {
1640: /* f2 = t*f1+f0 */
1641: _bshiftup2(f1,-1,f2);
1642: _addup2_destructive(f2,f0);
1643: /* cyclic change of the variables */
1644: t = f0; f0 = f1; f1 = f2; f2 = t;
1645: }
1646: _copyup2(f1,f);
1647: return 0;
1648: break;
1649: default:
1650: return -1;
1651: break;
1652: }
1653: }
1654:
1655: /*
1656: * f = an irreducible trinomial or pentanomial of degree d 'after' f
1657: * return value : 0 --- exists
1658: * 1 --- does not exist (exhaustion)
1659: */
1660:
1661: int _generate_irreducible_polynomial(UP2 f,int d)
1662: {
1663: int ret,i,j,k,nz,i0,j0,k0;
1664: int w;
1665: unsigned int *fd;
1666:
1667: /*
1668: * if f = x^d+x^i+1 then i0 <- i, j0 <- 0, k0 <-0.
1669: * if f = x^d+x^k+x^j+x^i+1 (k>j>i) then i0 <- i, j0 <- j, k0 <-k.
1670: * otherwise i0,j0,k0 is set to 0.
1671: */
1672:
1673: fd = f->b;
1674: w = (d>>5)+1;
1675: if ( f->w && (d==degup2(f)) ) {
1676: for ( nz = 0, i = d; i >= 0; i-- )
1677: if ( fd[i>>5]&(1<<(i&31)) ) nz++;
1678: switch ( nz ) {
1679: case 3:
1680: for ( i0 = 1; !(fd[i0>>5]&(1<<(i0&31))) ; i0++ );
1681: /* reset i0-th bit */
1682: fd[i0>>5] &= ~(1<<(i0&31));
1683: j0 = k0 = 0;
1684: break;
1685: case 5:
1686: for ( i0 = 1; !(fd[i0>>5]&(1<<(i0&31))) ; i0++ );
1687: /* reset i0-th bit */
1688: fd[i0>>5] &= ~(1<<(i0&31));
1689: for ( j0 = i0+1; !(fd[j0>>5]&(1<<(j0&31))) ; j0++ );
1690: /* reset j0-th bit */
1691: fd[j0>>5] &= ~(1<<(j0&31));
1692: for ( k0 = j0+1; !(fd[k0>>5]&(1<<(k0&31))) ; k0++ );
1693: /* reset k0-th bit */
1694: fd[k0>>5] &= ~(1<<(k0&31));
1695: break;
1696: default:
1697: f->w = 0; break;
1698: }
1699: } else
1700: f->w = 0;
1701:
1702: if ( !f->w ) {
1703: fd = f->b;
1704: f->w = w; fd[0] |= 1; fd[d>>5] |= (1<<(d&31));
1705: i0 = j0 = k0 = 0;
1706: }
1707: /* if j0 > 0 then f is already a pentanomial */
1708: if ( j0 > 0 ) goto PENTA;
1709:
1710: /* searching for an irreducible trinomial */
1711:
1712: for ( i = 1; 2*i <= d; i++ ) {
1713: /* skip the polynomials 'before' f */
1714: if ( i < i0 ) continue;
1715: if ( i == i0 ) { i0 = 0; continue; }
1716: /* set i-th bit */
1717: fd[i>>5] |= (1<<(i&31));
1718: ret = irredcheck_dddup2(f);
1719: if ( ret == 1 ) return 0;
1720: /* reset i-th bit */
1721: fd[i>>5] &= ~(1<<(i&31));
1722: }
1723:
1724: /* searching for an irreducible pentanomial */
1725: PENTA:
1726: for ( i = 1; i < d; i++ ) {
1727: /* skip the polynomials 'before' f */
1728: if ( i < i0 ) continue;
1729: if ( i == i0 ) i0 = 0;
1730: /* set i-th bit */
1731: fd[i>>5] |= (1<<(i&31));
1732: for ( j = i+1; j < d; j++ ) {
1733: /* skip the polynomials 'before' f */
1734: if ( j < j0 ) continue;
1735: if ( j == j0 ) j0 = 0;
1736: /* set j-th bit */
1737: fd[j>>5] |= (1<<(j&31));
1738: for ( k = j+1; k < d; k++ ) {
1739: /* skip the polynomials 'before' f */
1740: if ( k < k0 ) continue;
1741: else if ( k == k0 ) { k0 = 0; continue; }
1742: /* set k-th bit */
1743: fd[k>>5] |= (1<<(k&31));
1744: ret = irredcheck_dddup2(f);
1745: if ( ret == 1 ) return 0;
1746: /* reset k-th bit */
1747: fd[k>>5] &= ~(1<<(k&31));
1748: }
1749: /* reset j-th bit */
1750: fd[j>>5] &= ~(1<<(j&31));
1751: }
1752: /* reset i-th bit */
1753: fd[i>>5] &= ~(1<<(i&31));
1754: }
1755: /* exhausted */
1756: return 1;
1757: }
1758:
1759: /*
1760: * f = an irreducible trinomial or pentanomial of degree d 'after' f
1761: *
1762: * searching strategy:
1763: * trinomial x^d+x^i+1:
1764: * i is as small as possible.
1765: * trinomial x^d+x^i+x^j+x^k+1:
1766: * i is as small as possible.
1767: * For such i, j is as small as possible.
1768: * For such i and j, 'k' is as small as possible.
1769: *
1770: * return value : 0 --- exists
1771: * 1 --- does not exist (exhaustion)
1772: */
1773:
1774: int _generate_good_irreducible_polynomial(UP2 f,int d)
1775: {
1776: int ret,i,j,k,nz,i0,j0,k0;
1777: int w;
1778: unsigned int *fd;
1779:
1780: /*
1781: * if f = x^d+x^i+1 then i0 <- i, j0 <- 0, k0 <-0.
1782: * if f = x^d+x^k+x^j+x^i+1 (k>j>i) then i0 <- i, j0 <- j, k0 <-k.
1783: * otherwise i0,j0,k0 is set to 0.
1784: */
1785:
1786: fd = f->b;
1787: w = (d>>5)+1;
1788: if ( f->w && (d==degup2(f)) ) {
1789: for ( nz = 0, i = d; i >= 0; i-- )
1790: if ( fd[i>>5]&(1<<(i&31)) ) nz++;
1791: switch ( nz ) {
1792: case 3:
1793: for ( i0 = 1; !(fd[i0>>5]&(1<<(i0&31))) ; i0++ );
1794: /* reset i0-th bit */
1795: fd[i0>>5] &= ~(1<<(i0&31));
1796: j0 = k0 = 0;
1797: break;
1798: case 5:
1799: for ( i0 = 1; !(fd[i0>>5]&(1<<(i0&31))) ; i0++ );
1800: /* reset i0-th bit */
1801: fd[i0>>5] &= ~(1<<(i0&31));
1802: for ( j0 = i0+1; !(fd[j0>>5]&(1<<(j0&31))) ; j0++ );
1803: /* reset j0-th bit */
1804: fd[j0>>5] &= ~(1<<(j0&31));
1805: for ( k0 = j0+1; !(fd[k0>>5]&(1<<(k0&31))) ; k0++ );
1806: /* reset k0-th bit */
1807: fd[k0>>5] &= ~(1<<(k0&31));
1808: break;
1809: default:
1810: f->w = 0; break;
1811: }
1812: } else
1813: f->w = 0;
1814:
1815: if ( !f->w ) {
1816: fd = f->b;
1817: f->w = w; fd[0] |= 1; fd[d>>5] |= (1<<(d&31));
1818: i0 = j0 = k0 = 0;
1819: }
1820: /* if j0 > 0 then f is already a pentanomial */
1821: if ( j0 > 0 ) goto PENTA;
1822:
1823: /* searching for an irreducible trinomial */
1824:
1825: for ( i = 1; 2*i <= d; i++ ) {
1826: /* skip the polynomials 'before' f */
1827: if ( i < i0 ) continue;
1828: if ( i == i0 ) { i0 = 0; continue; }
1829: /* set i-th bit */
1830: fd[i>>5] |= (1<<(i&31));
1831: ret = irredcheck_dddup2(f);
1832: if ( ret == 1 ) return 0;
1833: /* reset i-th bit */
1834: fd[i>>5] &= ~(1<<(i&31));
1835: }
1836:
1837: /* searching for an irreducible pentanomial */
1838: PENTA:
1839: for ( i = 3; i < d; i++ ) {
1840: /* skip the polynomials 'before' f */
1841: if ( i < i0 ) continue;
1842: if ( i == i0 ) i0 = 0;
1843: /* set i-th bit */
1844: fd[i>>5] |= (1<<(i&31));
1845: for ( j = 2; j < i; j++ ) {
1846: /* skip the polynomials 'before' f */
1847: if ( j < j0 ) continue;
1848: if ( j == j0 ) j0 = 0;
1849: /* set j-th bit */
1850: fd[j>>5] |= (1<<(j&31));
1851: for ( k = 1; k < j; k++ ) {
1852: /* skip the polynomials 'before' f */
1853: if ( k < k0 ) continue;
1854: else if ( k == k0 ) { k0 = 0; continue; }
1855: /* set k-th bit */
1856: fd[k>>5] |= (1<<(k&31));
1857: ret = irredcheck_dddup2(f);
1858: if ( ret == 1 ) return 0;
1859: /* reset k-th bit */
1860: fd[k>>5] &= ~(1<<(k&31));
1861: }
1862: /* reset j-th bit */
1863: fd[j>>5] &= ~(1<<(j&31));
1864: }
1865: /* reset i-th bit */
1866: fd[i>>5] &= ~(1<<(i&31));
1867: }
1868: /* exhausted */
1869: return 1;
1870: }