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