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