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