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