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