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