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