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