![[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.19 ! noro 48: * $OpenXM: OpenXM_contrib2/asir2000/builtin/array.c,v 1.18 2001/09/17 01:18:34 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: int red_by_compress(m,p,r,hc,len)
1442: int m;
1443: unsigned int *p;
1.19 ! noro 1444: register struct oCM *r;
1.18 noro 1445: unsigned int hc;
1.19 ! noro 1446: register int len;
1.18 noro 1447: {
1.19 ! noro 1448: unsigned int up,lo;
1.18 noro 1449: unsigned int dmy;
1450: unsigned int *pj;
1451:
1452: p[r->index] = 0; r++;
1.19 ! noro 1453: for ( len--; len; len--, r++ ) {
1.18 noro 1454: pj = p+r->index;
1455: DMA0(r->c,hc,*pj,up,lo);
1456: if ( up ) {
1457: DSAB(m,up,lo,dmy,*pj);
1458: } else
1459: *pj = lo;
1460: }
1461: }
1462:
1463: /* p -= hc*r */
1464:
1465: int red_by_vect(m,p,r,hc,len)
1466: int m;
1467: unsigned int *p,*r;
1468: unsigned int hc;
1469: int len;
1470: {
1471: register unsigned int up,lo;
1472: unsigned int dmy;
1473:
1474: *p++ = 0; r++; len--;
1475: for ( ; len; len--, r++, p++ )
1476: if ( *r ) {
1477: DMA0(*r,hc,*p,up,lo);
1478: if ( up ) {
1479: DSAB(m,up,lo,dmy,*p);
1480: } else
1481: *p = lo;
1482: }
1483: }
1484:
1.15 noro 1485: void reduce_sp_by_red_mod_compress (sp,redmat,ind,nred,col,md)
1486: int *sp;
1487: CDP *redmat;
1488: int *ind;
1489: int nred,col;
1490: int md;
1491: {
1.18 noro 1492: int i,j,k,len;
1493: unsigned int *tj;
1.15 noro 1494: CDP ri;
1.18 noro 1495: unsigned int hc,up,lo,up1,lo1,c;
1496: unsigned int *usp;
1497: struct oCM *rib;
1.15 noro 1498:
1.18 noro 1499: usp = (unsigned int *)sp;
1.15 noro 1500: /* reduce the spolys by redmat */
1501: for ( i = nred-1; i >= 0; i-- ) {
1502: /* reduce sp by redmat[i] */
1.18 noro 1503: usp[ind[i]] %= md;
1504: if ( hc = usp[ind[i]] ) {
1.15 noro 1505: /* sp = sp-hc*redmat[i] */
1506: hc = md-hc;
1507: ri = redmat[i];
1508: len = ri->len;
1.18 noro 1509: red_by_compress(md,usp,ri->body,hc,len);
1.4 noro 1510: }
1511: }
1.18 noro 1512: for ( i = 0; i < col; i++ )
1513: if ( usp[i] >= md )
1514: usp[i] %= md;
1.4 noro 1515: }
1516:
1517: #define ONE_STEP2 if ( zzz = *pk ) { DMAR(zzz,a,*tk,md,*tk) } pk++; tk++;
1518:
1.18 noro 1519: int generic_gauss_elim_mod(mat0,row,col,md,colstat)
1520: int **mat0;
1.1 noro 1521: int row,col,md;
1522: int *colstat;
1523: {
1.4 noro 1524: int i,j,k,l,inv,a,rank,zzz;
1.18 noro 1525: unsigned int *t,*pivot,*pk,*tk;
1526: unsigned int **mat;
1.1 noro 1527:
1.18 noro 1528: mat = (unsigned int **)mat0;
1.1 noro 1529: for ( rank = 0, j = 0; j < col; j++ ) {
1.18 noro 1530: for ( i = rank; i < row; i++ )
1531: mat[i][j] %= md;
1532: for ( i = rank; i < row; i++ )
1533: if ( mat[i][j] )
1534: break;
1.1 noro 1535: if ( i == row ) {
1536: colstat[j] = 0;
1537: continue;
1538: } else
1539: colstat[j] = 1;
1540: if ( i != rank ) {
1541: t = mat[i]; mat[i] = mat[rank]; mat[rank] = t;
1542: }
1543: pivot = mat[rank];
1544: inv = invm(pivot[j],md);
1.4 noro 1545: for ( k = j, pk = pivot+k; k < col; k++, pk++ )
1546: if ( *pk ) {
1.18 noro 1547: if ( *pk >= md )
1548: *pk %= md;
1.4 noro 1549: DMAR(*pk,inv,0,md,*pk)
1.1 noro 1550: }
1551: for ( i = rank+1; i < row; i++ ) {
1552: t = mat[i];
1.18 noro 1553: if ( a = t[j] )
1554: red_by_vect(md,t+j,pivot+j,md-a,col-j);
1.1 noro 1555: }
1556: rank++;
1557: }
1558: for ( j = col-1, l = rank-1; j >= 0; j-- )
1559: if ( colstat[j] ) {
1560: pivot = mat[l];
1561: for ( i = 0; i < l; i++ ) {
1562: t = mat[i];
1.18 noro 1563: t[j] %= md;
1564: if ( a = t[j] )
1565: red_by_vect(md,t+j,pivot+j,md-a,col-j);
1.1 noro 1566: }
1567: l--;
1.18 noro 1568: }
1569: for ( j = 0, l = 0; l < rank; j++ )
1570: if ( colstat[j] ) {
1571: t = mat[l];
1572: for ( k = j; k < col; k++ )
1573: if ( t[k] >= md )
1574: t[k] %= md;
1575: l++;
1.1 noro 1576: }
1577: return rank;
1578: }
1579:
1580: /* LU decomposition; a[i][i] = 1/U[i][i] */
1581:
1582: int lu_gfmmat(mat,md,perm)
1583: GFMMAT mat;
1584: unsigned int md;
1585: int *perm;
1586: {
1587: int row,col;
1588: int i,j,k,l;
1589: unsigned int *t,*pivot;
1590: unsigned int **a;
1591: unsigned int inv,m;
1592:
1593: row = mat->row; col = mat->col;
1594: a = mat->body;
1595: bzero(perm,row*sizeof(int));
1596:
1597: for ( i = 0; i < row; i++ )
1598: perm[i] = i;
1599: for ( k = 0; k < col; k++ ) {
1600: for ( i = k; i < row && !a[i][k]; i++ );
1601: if ( i == row )
1602: return 0;
1603: if ( i != k ) {
1604: j = perm[i]; perm[i] = perm[k]; perm[k] = j;
1605: t = a[i]; a[i] = a[k]; a[k] = t;
1606: }
1607: pivot = a[k];
1608: pivot[k] = inv = invm(pivot[k],md);
1609: for ( i = k+1; i < row; i++ ) {
1610: t = a[i];
1611: if ( m = t[k] ) {
1612: DMAR(inv,m,0,md,t[k])
1613: for ( j = k+1, m = md - t[k]; j < col; j++ )
1614: if ( pivot[j] ) {
1.8 noro 1615: unsigned int tj;
1616:
1617: DMAR(m,pivot[j],t[j],md,tj)
1618: t[j] = tj;
1.1 noro 1619: }
1620: }
1621: }
1622: }
1623: return 1;
1624: }
1625:
1.3 noro 1626: /*
1627: Input
1628: a: a row x col matrix
1629: md : a modulus
1630:
1631: Output:
1632: return : d = the rank of mat
1633: a[0..(d-1)][0..(d-1)] : LU decomposition (a[i][i] = 1/U[i][i])
1634: rinfo: array of length row
1635: cinfo: array of length col
1636: i-th row in new a <-> rinfo[i]-th row in old a
1637: cinfo[j]=1 <=> j-th column is contained in the LU decomp.
1638: */
1639:
1640: int find_lhs_and_lu_mod(a,row,col,md,rinfo,cinfo)
1641: unsigned int **a;
1642: unsigned int md;
1643: int **rinfo,**cinfo;
1644: {
1645: int i,j,k,l,d;
1646: int *rp,*cp;
1647: unsigned int *t,*pivot;
1648: unsigned int inv,m;
1649:
1650: *rinfo = rp = (int *)MALLOC_ATOMIC(row*sizeof(int));
1651: *cinfo = cp = (int *)MALLOC_ATOMIC(col*sizeof(int));
1652: for ( i = 0; i < row; i++ )
1653: rp[i] = i;
1654: for ( k = 0, d = 0; k < col; k++ ) {
1655: for ( i = d; i < row && !a[i][k]; i++ );
1656: if ( i == row ) {
1657: cp[k] = 0;
1658: continue;
1659: } else
1660: cp[k] = 1;
1661: if ( i != d ) {
1662: j = rp[i]; rp[i] = rp[d]; rp[d] = j;
1663: t = a[i]; a[i] = a[d]; a[d] = t;
1664: }
1665: pivot = a[d];
1666: pivot[k] = inv = invm(pivot[k],md);
1667: for ( i = d+1; i < row; i++ ) {
1668: t = a[i];
1669: if ( m = t[k] ) {
1670: DMAR(inv,m,0,md,t[k])
1671: for ( j = k+1, m = md - t[k]; j < col; j++ )
1672: if ( pivot[j] ) {
1.8 noro 1673: unsigned int tj;
1674: DMAR(m,pivot[j],t[j],md,tj)
1675: t[j] = tj;
1.3 noro 1676: }
1677: }
1678: }
1679: d++;
1680: }
1681: return d;
1682: }
1683:
1684: /*
1685: Input
1686: a : n x n matrix; a result of LU-decomposition
1687: md : modulus
1688: b : n x l matrix
1689: Output
1690: b = a^(-1)b
1691: */
1692:
1693: void solve_by_lu_mod(a,n,md,b,l)
1694: int **a;
1695: int n;
1696: int md;
1697: int **b;
1698: int l;
1699: {
1700: unsigned int *y,*c;
1701: int i,j,k;
1702: unsigned int t,m,m2;
1703:
1704: y = (int *)MALLOC_ATOMIC(n*sizeof(int));
1705: c = (int *)MALLOC_ATOMIC(n*sizeof(int));
1706: m2 = md>>1;
1707: for ( k = 0; k < l; k++ ) {
1708: /* copy b[.][k] to c */
1709: for ( i = 0; i < n; i++ )
1710: c[i] = (unsigned int)b[i][k];
1711: /* solve Ly=c */
1712: for ( i = 0; i < n; i++ ) {
1713: for ( t = c[i], j = 0; j < i; j++ )
1714: if ( a[i][j] ) {
1715: m = md - a[i][j];
1716: DMAR(m,y[j],t,md,t)
1717: }
1718: y[i] = t;
1719: }
1720: /* solve Uc=y */
1721: for ( i = n-1; i >= 0; i-- ) {
1722: for ( t = y[i], j =i+1; j < n; j++ )
1723: if ( a[i][j] ) {
1724: m = md - a[i][j];
1725: DMAR(m,c[j],t,md,t)
1726: }
1727: /* a[i][i] = 1/U[i][i] */
1728: DMAR(t,a[i][i],0,md,c[i])
1729: }
1730: /* copy c to b[.][k] with normalization */
1731: for ( i = 0; i < n; i++ )
1732: b[i][k] = (int)(c[i]>m2 ? c[i]-md : c[i]);
1733: }
1734: }
1735:
1.1 noro 1736: void Pleqm1(arg,rp)
1737: NODE arg;
1738: VECT *rp;
1739: {
1740: MAT m;
1741: VECT vect;
1742: pointer **mat;
1743: Q *v;
1744: Q q;
1745: int **wmat;
1746: int md,i,j,row,col,t,n,status;
1747:
1748: asir_assert(ARG0(arg),O_MAT,"leqm1");
1749: asir_assert(ARG1(arg),O_N,"leqm1");
1750: m = (MAT)ARG0(arg); md = QTOS((Q)ARG1(arg));
1751: row = m->row; col = m->col; mat = m->body;
1752: wmat = (int **)almat(row,col);
1753: for ( i = 0; i < row; i++ )
1754: for ( j = 0; j < col; j++ )
1755: if ( q = (Q)mat[i][j] ) {
1756: t = rem(NM(q),md);
1757: if ( SGN(q) < 0 )
1758: t = (md - t) % md;
1759: wmat[i][j] = t;
1760: } else
1761: wmat[i][j] = 0;
1762: status = gauss_elim_mod1(wmat,row,col,md);
1763: if ( status < 0 )
1764: *rp = 0;
1765: else if ( status > 0 )
1766: *rp = (VECT)ONE;
1767: else {
1768: n = col - 1;
1769: MKVECT(vect,n);
1770: for ( i = 0, v = (Q *)vect->body; i < n; i++ ) {
1771: t = (md-wmat[i][n])%md; STOQ(t,v[i]);
1772: }
1773: *rp = vect;
1774: }
1775: }
1776:
1777: gauss_elim_mod1(mat,row,col,md)
1778: int **mat;
1779: int row,col,md;
1780: {
1781: int i,j,k,inv,a,n;
1782: int *t,*pivot;
1783:
1784: n = col - 1;
1785: for ( j = 0; j < n; j++ ) {
1786: for ( i = j; i < row && !mat[i][j]; i++ );
1787: if ( i == row )
1788: return 1;
1789: if ( i != j ) {
1790: t = mat[i]; mat[i] = mat[j]; mat[j] = t;
1791: }
1792: pivot = mat[j];
1793: inv = invm(pivot[j],md);
1794: for ( k = j; k <= n; k++ )
1795: pivot[k] = dmar(pivot[k],inv,0,md);
1796: for ( i = j+1; i < row; i++ ) {
1797: t = mat[i];
1798: if ( i != j && (a = t[j]) )
1799: for ( k = j, a = md - a; k <= n; k++ )
1800: t[k] = dmar(pivot[k],a,t[k],md);
1801: }
1802: }
1803: for ( i = n; i < row && !mat[i][n]; i++ );
1804: if ( i == row ) {
1805: for ( j = n-1; j >= 0; j-- ) {
1806: for ( i = j-1, a = (md-mat[j][n])%md; i >= 0; i-- ) {
1807: mat[i][n] = dmar(mat[i][j],a,mat[i][n],md);
1808: mat[i][j] = 0;
1809: }
1810: }
1811: return 0;
1812: } else
1813: return -1;
1814: }
1815:
1816: void Pgeninvm(arg,rp)
1817: NODE arg;
1818: LIST *rp;
1819: {
1820: MAT m;
1821: pointer **mat;
1822: Q **tmat;
1823: Q q;
1824: unsigned int **wmat;
1825: int md,i,j,row,col,t,status;
1826: MAT mat1,mat2;
1827: NODE node1,node2;
1828:
1829: asir_assert(ARG0(arg),O_MAT,"leqm1");
1830: asir_assert(ARG1(arg),O_N,"leqm1");
1831: m = (MAT)ARG0(arg); md = QTOS((Q)ARG1(arg));
1832: row = m->row; col = m->col; mat = m->body;
1833: wmat = (unsigned int **)almat(row,col+row);
1834: for ( i = 0; i < row; i++ ) {
1835: bzero((char *)wmat[i],(col+row)*sizeof(int));
1836: for ( j = 0; j < col; j++ )
1837: if ( q = (Q)mat[i][j] ) {
1838: t = rem(NM(q),md);
1839: if ( SGN(q) < 0 )
1840: t = (md - t) % md;
1841: wmat[i][j] = t;
1842: }
1843: wmat[i][col+i] = 1;
1844: }
1845: status = gauss_elim_geninv_mod(wmat,row,col,md);
1846: if ( status > 0 )
1847: *rp = 0;
1848: else {
1849: MKMAT(mat1,col,row); MKMAT(mat2,row-col,row);
1850: for ( i = 0, tmat = (Q **)mat1->body; i < col; i++ )
1851: for ( j = 0; j < row; j++ )
1852: STOQ(wmat[i][j+col],tmat[i][j]);
1853: for ( tmat = (Q **)mat2->body; i < row; i++ )
1854: for ( j = 0; j < row; j++ )
1855: STOQ(wmat[i][j+col],tmat[i-col][j]);
1856: MKNODE(node2,mat2,0); MKNODE(node1,mat1,node2); MKLIST(*rp,node1);
1857: }
1858: }
1859:
1860: int gauss_elim_geninv_mod(mat,row,col,md)
1861: unsigned int **mat;
1862: int row,col,md;
1863: {
1864: int i,j,k,inv,a,n,m;
1865: unsigned int *t,*pivot;
1866:
1867: n = col; m = row+col;
1868: for ( j = 0; j < n; j++ ) {
1869: for ( i = j; i < row && !mat[i][j]; i++ );
1870: if ( i == row )
1871: return 1;
1872: if ( i != j ) {
1873: t = mat[i]; mat[i] = mat[j]; mat[j] = t;
1874: }
1875: pivot = mat[j];
1876: inv = invm(pivot[j],md);
1877: for ( k = j; k < m; k++ )
1878: pivot[k] = dmar(pivot[k],inv,0,md);
1879: for ( i = j+1; i < row; i++ ) {
1880: t = mat[i];
1881: if ( a = t[j] )
1882: for ( k = j, a = md - a; k < m; k++ )
1883: t[k] = dmar(pivot[k],a,t[k],md);
1884: }
1885: }
1886: for ( j = n-1; j >= 0; j-- ) {
1887: pivot = mat[j];
1888: for ( i = j-1; i >= 0; i-- ) {
1889: t = mat[i];
1890: if ( a = t[j] )
1891: for ( k = j, a = md - a; k < m; k++ )
1892: t[k] = dmar(pivot[k],a,t[k],md);
1893: }
1894: }
1895: return 0;
1896: }
1897:
1898: void Psolve_by_lu_gfmmat(arg,rp)
1899: NODE arg;
1900: VECT *rp;
1901: {
1902: GFMMAT lu;
1903: Q *perm,*rhs,*v;
1904: int n,i;
1905: unsigned int md;
1906: unsigned int *b,*sol;
1907: VECT r;
1908:
1909: lu = (GFMMAT)ARG0(arg);
1910: perm = (Q *)BDY((VECT)ARG1(arg));
1911: rhs = (Q *)BDY((VECT)ARG2(arg));
1912: md = (unsigned int)QTOS((Q)ARG3(arg));
1913: n = lu->col;
1914: b = (unsigned int *)MALLOC_ATOMIC(n*sizeof(int));
1915: sol = (unsigned int *)MALLOC_ATOMIC(n*sizeof(int));
1916: for ( i = 0; i < n; i++ )
1917: b[i] = QTOS(rhs[QTOS(perm[i])]);
1918: solve_by_lu_gfmmat(lu,md,b,sol);
1919: MKVECT(r,n);
1920: for ( i = 0, v = (Q *)r->body; i < n; i++ )
1921: STOQ(sol[i],v[i]);
1922: *rp = r;
1923: }
1924:
1925: void solve_by_lu_gfmmat(lu,md,b,x)
1926: GFMMAT lu;
1927: unsigned int md;
1928: unsigned int *b;
1929: unsigned int *x;
1930: {
1931: int n;
1932: unsigned int **a;
1933: unsigned int *y;
1934: int i,j;
1935: unsigned int t,m;
1936:
1937: n = lu->col;
1938: a = lu->body;
1939: y = (unsigned int *)MALLOC_ATOMIC(n*sizeof(int));
1940: /* solve Ly=b */
1941: for ( i = 0; i < n; i++ ) {
1942: for ( t = b[i], j = 0; j < i; j++ )
1943: if ( a[i][j] ) {
1944: m = md - a[i][j];
1945: DMAR(m,y[j],t,md,t)
1946: }
1947: y[i] = t;
1948: }
1949: /* solve Ux=y */
1950: for ( i = n-1; i >= 0; i-- ) {
1951: for ( t = y[i], j =i+1; j < n; j++ )
1952: if ( a[i][j] ) {
1953: m = md - a[i][j];
1954: DMAR(m,x[j],t,md,t)
1955: }
1956: /* a[i][i] = 1/U[i][i] */
1957: DMAR(t,a[i][i],0,md,x[i])
1958: }
1959: }
1960:
1961: void Plu_gfmmat(arg,rp)
1962: NODE arg;
1963: LIST *rp;
1964: {
1965: MAT m;
1966: GFMMAT mm;
1967: unsigned int md;
1968: int i,row,col,status;
1969: int *iperm;
1970: Q *v;
1971: VECT perm;
1972: NODE n0;
1973:
1974: asir_assert(ARG0(arg),O_MAT,"mat_to_gfmmat");
1975: asir_assert(ARG1(arg),O_N,"mat_to_gfmmat");
1976: m = (MAT)ARG0(arg); md = (unsigned int)QTOS((Q)ARG1(arg));
1977: mat_to_gfmmat(m,md,&mm);
1978: row = m->row;
1979: col = m->col;
1980: iperm = (int *)MALLOC_ATOMIC(row*sizeof(int));
1981: status = lu_gfmmat(mm,md,iperm);
1982: if ( !status )
1983: n0 = 0;
1984: else {
1985: MKVECT(perm,row);
1986: for ( i = 0, v = (Q *)perm->body; i < row; i++ )
1987: STOQ(iperm[i],v[i]);
1988: n0 = mknode(2,mm,perm);
1989: }
1990: MKLIST(*rp,n0);
1991: }
1992:
1993: void Pmat_to_gfmmat(arg,rp)
1994: NODE arg;
1995: GFMMAT *rp;
1996: {
1997: MAT m;
1998: unsigned int md;
1999:
2000: asir_assert(ARG0(arg),O_MAT,"mat_to_gfmmat");
2001: asir_assert(ARG1(arg),O_N,"mat_to_gfmmat");
2002: m = (MAT)ARG0(arg); md = (unsigned int)QTOS((Q)ARG1(arg));
2003: mat_to_gfmmat(m,md,rp);
2004: }
2005:
2006: void mat_to_gfmmat(m,md,rp)
2007: MAT m;
2008: unsigned int md;
2009: GFMMAT *rp;
2010: {
2011: unsigned int **wmat;
2012: unsigned int t;
2013: Q **mat;
2014: Q q;
2015: int i,j,row,col;
2016:
2017: row = m->row; col = m->col; mat = (Q **)m->body;
2018: wmat = (unsigned int **)almat(row,col);
2019: for ( i = 0; i < row; i++ ) {
2020: bzero((char *)wmat[i],col*sizeof(unsigned int));
2021: for ( j = 0; j < col; j++ )
2022: if ( q = mat[i][j] ) {
2023: t = (unsigned int)rem(NM(q),md);
2024: if ( SGN(q) < 0 )
2025: t = (md - t) % md;
2026: wmat[i][j] = t;
2027: }
2028: }
2029: TOGFMMAT(row,col,wmat,*rp);
2030: }
2031:
2032: void Pgeninvm_swap(arg,rp)
2033: NODE arg;
2034: LIST *rp;
2035: {
2036: MAT m;
2037: pointer **mat;
2038: Q **tmat;
2039: Q *tvect;
2040: Q q;
2041: unsigned int **wmat,**invmat;
2042: int *index;
2043: unsigned int t,md;
2044: int i,j,row,col,status;
2045: MAT mat1;
2046: VECT vect1;
2047: NODE node1,node2;
2048:
2049: asir_assert(ARG0(arg),O_MAT,"geninvm_swap");
2050: asir_assert(ARG1(arg),O_N,"geninvm_swap");
2051: m = (MAT)ARG0(arg); md = QTOS((Q)ARG1(arg));
2052: row = m->row; col = m->col; mat = m->body;
2053: wmat = (unsigned int **)almat(row,col+row);
2054: for ( i = 0; i < row; i++ ) {
2055: bzero((char *)wmat[i],(col+row)*sizeof(int));
2056: for ( j = 0; j < col; j++ )
2057: if ( q = (Q)mat[i][j] ) {
2058: t = (unsigned int)rem(NM(q),md);
2059: if ( SGN(q) < 0 )
2060: t = (md - t) % md;
2061: wmat[i][j] = t;
2062: }
2063: wmat[i][col+i] = 1;
2064: }
2065: status = gauss_elim_geninv_mod_swap(wmat,row,col,md,&invmat,&index);
2066: if ( status > 0 )
2067: *rp = 0;
2068: else {
2069: MKMAT(mat1,col,col);
2070: for ( i = 0, tmat = (Q **)mat1->body; i < col; i++ )
2071: for ( j = 0; j < col; j++ )
2072: UTOQ(invmat[i][j],tmat[i][j]);
2073: MKVECT(vect1,row);
2074: for ( i = 0, tvect = (Q *)vect1->body; i < row; i++ )
2075: STOQ(index[i],tvect[i]);
2076: MKNODE(node2,vect1,0); MKNODE(node1,mat1,node2); MKLIST(*rp,node1);
2077: }
2078: }
2079:
2080: gauss_elim_geninv_mod_swap(mat,row,col,md,invmatp,indexp)
2081: unsigned int **mat;
2082: int row,col;
2083: unsigned int md;
2084: unsigned int ***invmatp;
2085: int **indexp;
2086: {
2087: int i,j,k,inv,a,n,m;
2088: unsigned int *t,*pivot,*s;
2089: int *index;
2090: unsigned int **invmat;
2091:
2092: n = col; m = row+col;
2093: *indexp = index = (int *)MALLOC_ATOMIC(row*sizeof(int));
2094: for ( i = 0; i < row; i++ )
2095: index[i] = i;
2096: for ( j = 0; j < n; j++ ) {
2097: for ( i = j; i < row && !mat[i][j]; i++ );
2098: if ( i == row ) {
2099: *indexp = 0; *invmatp = 0; return 1;
2100: }
2101: if ( i != j ) {
2102: t = mat[i]; mat[i] = mat[j]; mat[j] = t;
2103: k = index[i]; index[i] = index[j]; index[j] = k;
2104: }
2105: pivot = mat[j];
2106: inv = (unsigned int)invm(pivot[j],md);
2107: for ( k = j; k < m; k++ )
2108: if ( pivot[k] )
2109: pivot[k] = (unsigned int)dmar(pivot[k],inv,0,md);
2110: for ( i = j+1; i < row; i++ ) {
2111: t = mat[i];
2112: if ( a = t[j] )
2113: for ( k = j, a = md - a; k < m; k++ )
2114: if ( pivot[k] )
2115: t[k] = dmar(pivot[k],a,t[k],md);
2116: }
2117: }
2118: for ( j = n-1; j >= 0; j-- ) {
2119: pivot = mat[j];
2120: for ( i = j-1; i >= 0; i-- ) {
2121: t = mat[i];
2122: if ( a = t[j] )
2123: for ( k = j, a = md - a; k < m; k++ )
2124: if ( pivot[k] )
2125: t[k] = dmar(pivot[k],a,t[k],md);
2126: }
2127: }
2128: *invmatp = invmat = (unsigned int **)almat(col,col);
2129: for ( i = 0; i < col; i++ )
2130: for ( j = 0, s = invmat[i], t = mat[i]; j < col; j++ )
2131: s[j] = t[col+index[j]];
2132: return 0;
2133: }
2134:
2135: void _addn(N,N,N);
2136: int _subn(N,N,N);
2137: void _muln(N,N,N);
2138:
2139: void inner_product_int(a,b,n,r)
2140: Q *a,*b;
2141: int n;
2142: Q *r;
2143: {
2144: int la,lb,i;
2145: int sgn,sgn1;
2146: N wm,wma,sum,t;
2147:
2148: for ( la = lb = 0, i = 0; i < n; i++ ) {
2149: if ( a[i] )
2150: if ( DN(a[i]) )
2151: error("inner_product_int : invalid argument");
2152: else
2153: la = MAX(PL(NM(a[i])),la);
2154: if ( b[i] )
2155: if ( DN(b[i]) )
2156: error("inner_product_int : invalid argument");
2157: else
2158: lb = MAX(PL(NM(b[i])),lb);
2159: }
2160: sgn = 0;
2161: sum= NALLOC(la+lb+2);
2162: bzero((char *)sum,(la+lb+3)*sizeof(unsigned int));
2163: wm = NALLOC(la+lb+2);
2164: wma = NALLOC(la+lb+2);
2165: for ( i = 0; i < n; i++ ) {
2166: if ( !a[i] || !b[i] )
2167: continue;
2168: _muln(NM(a[i]),NM(b[i]),wm);
2169: sgn1 = SGN(a[i])*SGN(b[i]);
2170: if ( !sgn ) {
2171: sgn = sgn1;
2172: t = wm; wm = sum; sum = t;
2173: } else if ( sgn == sgn1 ) {
2174: _addn(sum,wm,wma);
2175: if ( !PL(wma) )
2176: sgn = 0;
2177: t = wma; wma = sum; sum = t;
2178: } else {
2179: /* sgn*sum+sgn1*wm = sgn*(sum-wm) */
2180: sgn *= _subn(sum,wm,wma);
2181: t = wma; wma = sum; sum = t;
2182: }
2183: }
2184: GC_free(wm);
2185: GC_free(wma);
2186: if ( !sgn ) {
2187: GC_free(sum);
2188: *r = 0;
2189: } else
2190: NTOQ(sum,sgn,*r);
2191: }
2192:
1.3 noro 2193: /* (k,l) element of a*b where a: .x n matrix, b: n x . integer matrix */
2194:
2195: void inner_product_mat_int_mod(a,b,n,k,l,r)
2196: Q **a;
2197: int **b;
2198: int n,k,l;
2199: Q *r;
2200: {
2201: int la,lb,i;
2202: int sgn,sgn1;
2203: N wm,wma,sum,t;
2204: Q aki;
2205: int bil,bilsgn;
2206: struct oN tn;
2207:
2208: for ( la = 0, i = 0; i < n; i++ ) {
2209: if ( aki = a[k][i] )
2210: if ( DN(aki) )
2211: error("inner_product_int : invalid argument");
2212: else
2213: la = MAX(PL(NM(aki)),la);
2214: }
2215: lb = 1;
2216: sgn = 0;
2217: sum= NALLOC(la+lb+2);
2218: bzero((char *)sum,(la+lb+3)*sizeof(unsigned int));
2219: wm = NALLOC(la+lb+2);
2220: wma = NALLOC(la+lb+2);
2221: for ( i = 0; i < n; i++ ) {
2222: if ( !(aki = a[k][i]) || !(bil = b[i][l]) )
2223: continue;
2224: tn.p = 1;
2225: if ( bil > 0 ) {
2226: tn.b[0] = bil; bilsgn = 1;
2227: } else {
2228: tn.b[0] = -bil; bilsgn = -1;
2229: }
2230: _muln(NM(aki),&tn,wm);
2231: sgn1 = SGN(aki)*bilsgn;
2232: if ( !sgn ) {
2233: sgn = sgn1;
2234: t = wm; wm = sum; sum = t;
2235: } else if ( sgn == sgn1 ) {
2236: _addn(sum,wm,wma);
2237: if ( !PL(wma) )
2238: sgn = 0;
2239: t = wma; wma = sum; sum = t;
2240: } else {
2241: /* sgn*sum+sgn1*wm = sgn*(sum-wm) */
2242: sgn *= _subn(sum,wm,wma);
2243: t = wma; wma = sum; sum = t;
2244: }
2245: }
2246: GC_free(wm);
2247: GC_free(wma);
2248: if ( !sgn ) {
2249: GC_free(sum);
2250: *r = 0;
2251: } else
2252: NTOQ(sum,sgn,*r);
2253: }
2254:
1.1 noro 2255: void Pmul_mat_vect_int(arg,rp)
2256: NODE arg;
2257: VECT *rp;
2258: {
2259: MAT mat;
2260: VECT vect,r;
2261: int row,col,i;
2262:
2263: mat = (MAT)ARG0(arg);
2264: vect = (VECT)ARG1(arg);
2265: row = mat->row;
2266: col = mat->col;
2267: MKVECT(r,row);
2268: for ( i = 0; i < row; i++ )
2269: inner_product_int(mat->body[i],vect->body,col,&r->body[i]);
2270: *rp = r;
2271: }
2272:
2273: void Pnbpoly_up2(arg,rp)
2274: NODE arg;
2275: GF2N *rp;
2276: {
2277: int m,type,ret;
2278: UP2 r;
2279:
2280: m = QTOS((Q)ARG0(arg));
2281: type = QTOS((Q)ARG1(arg));
2282: ret = generate_ONB_polynomial(&r,m,type);
2283: if ( ret == 0 )
2284: MKGF2N(r,*rp);
2285: else
2286: *rp = 0;
2287: }
2288:
2289: void Px962_irredpoly_up2(arg,rp)
2290: NODE arg;
2291: GF2N *rp;
2292: {
2293: int m,type,ret,w;
2294: GF2N prev;
2295: UP2 r;
2296:
2297: m = QTOS((Q)ARG0(arg));
2298: prev = (GF2N)ARG1(arg);
2299: if ( !prev ) {
2300: w = (m>>5)+1; NEWUP2(r,w); r->w = 0;
2301: bzero((char *)r->b,w*sizeof(unsigned int));
2302: } else {
2303: r = prev->body;
2304: if ( degup2(r) != m ) {
2305: w = (m>>5)+1; NEWUP2(r,w); r->w = 0;
2306: bzero((char *)r->b,w*sizeof(unsigned int));
2307: }
2308: }
2309: ret = _generate_irreducible_polynomial(r,m,type);
2310: if ( ret == 0 )
2311: MKGF2N(r,*rp);
2312: else
2313: *rp = 0;
2314: }
2315:
2316: void Pirredpoly_up2(arg,rp)
2317: NODE arg;
2318: GF2N *rp;
2319: {
2320: int m,type,ret,w;
2321: GF2N prev;
2322: UP2 r;
2323:
2324: m = QTOS((Q)ARG0(arg));
2325: prev = (GF2N)ARG1(arg);
2326: if ( !prev ) {
2327: w = (m>>5)+1; NEWUP2(r,w); r->w = 0;
2328: bzero((char *)r->b,w*sizeof(unsigned int));
2329: } else {
2330: r = prev->body;
2331: if ( degup2(r) != m ) {
2332: w = (m>>5)+1; NEWUP2(r,w); r->w = 0;
2333: bzero((char *)r->b,w*sizeof(unsigned int));
2334: }
2335: }
2336: ret = _generate_good_irreducible_polynomial(r,m,type);
2337: if ( ret == 0 )
2338: MKGF2N(r,*rp);
2339: else
2340: *rp = 0;
2341: }
2342:
2343: /*
2344: * f = type 'type' normal polynomial of degree m if exists
2345: * IEEE P1363 A.7.2
2346: *
2347: * return value : 0 --- exists
2348: * 1 --- does not exist
2349: * -1 --- failure (memory allocation error)
2350: */
2351:
2352: int generate_ONB_polynomial(UP2 *rp,int m,int type)
2353: {
2354: int i,r;
2355: int w;
2356: UP2 f,f0,f1,f2,t;
2357:
2358: w = (m>>5)+1;
2359: switch ( type ) {
2360: case 1:
2361: if ( !TypeT_NB_check(m,1) ) return 1;
2362: NEWUP2(f,w); *rp = f; f->w = w;
2363: /* set all the bits */
2364: for ( i = 0; i < w; i++ )
2365: f->b[i] = 0xffffffff;
2366: /* mask the top word if necessary */
2367: if ( r = (m+1)&31 )
2368: f->b[w-1] &= (1<<r)-1;
2369: return 0;
2370: break;
2371: case 2:
2372: if ( !TypeT_NB_check(m,2) ) return 1;
2373: NEWUP2(f,w); *rp = f;
2374: W_NEWUP2(f0,w);
2375: W_NEWUP2(f1,w);
2376: W_NEWUP2(f2,w);
2377:
2378: /* recursion for genrating Type II normal polynomial */
2379:
2380: /* f0 = 1, f1 = t+1 */
2381: f0->w = 1; f0->b[0] = 1;
2382: f1->w = 1; f1->b[0] = 3;
2383: for ( i = 2; i <= m; i++ ) {
2384: /* f2 = t*f1+f0 */
2385: _bshiftup2(f1,-1,f2);
2386: _addup2_destructive(f2,f0);
2387: /* cyclic change of the variables */
2388: t = f0; f0 = f1; f1 = f2; f2 = t;
2389: }
2390: _copyup2(f1,f);
2391: return 0;
2392: break;
2393: default:
2394: return -1;
2395: break;
2396: }
2397: }
2398:
2399: /*
2400: * f = an irreducible trinomial or pentanomial of degree d 'after' f
2401: * return value : 0 --- exists
2402: * 1 --- does not exist (exhaustion)
2403: */
2404:
2405: int _generate_irreducible_polynomial(UP2 f,int d)
2406: {
2407: int ret,i,j,k,nz,i0,j0,k0;
2408: int w;
2409: unsigned int *fd;
2410:
2411: /*
2412: * if f = x^d+x^i+1 then i0 <- i, j0 <- 0, k0 <-0.
2413: * if f = x^d+x^k+x^j+x^i+1 (k>j>i) then i0 <- i, j0 <- j, k0 <-k.
2414: * otherwise i0,j0,k0 is set to 0.
2415: */
2416:
2417: fd = f->b;
2418: w = (d>>5)+1;
2419: if ( f->w && (d==degup2(f)) ) {
2420: for ( nz = 0, i = d; i >= 0; i-- )
2421: if ( fd[i>>5]&(1<<(i&31)) ) nz++;
2422: switch ( nz ) {
2423: case 3:
2424: for ( i0 = 1; !(fd[i0>>5]&(1<<(i0&31))) ; i0++ );
2425: /* reset i0-th bit */
2426: fd[i0>>5] &= ~(1<<(i0&31));
2427: j0 = k0 = 0;
2428: break;
2429: case 5:
2430: for ( i0 = 1; !(fd[i0>>5]&(1<<(i0&31))) ; i0++ );
2431: /* reset i0-th bit */
2432: fd[i0>>5] &= ~(1<<(i0&31));
2433: for ( j0 = i0+1; !(fd[j0>>5]&(1<<(j0&31))) ; j0++ );
2434: /* reset j0-th bit */
2435: fd[j0>>5] &= ~(1<<(j0&31));
2436: for ( k0 = j0+1; !(fd[k0>>5]&(1<<(k0&31))) ; k0++ );
2437: /* reset k0-th bit */
2438: fd[k0>>5] &= ~(1<<(k0&31));
2439: break;
2440: default:
2441: f->w = 0; break;
2442: }
2443: } else
2444: f->w = 0;
2445:
2446: if ( !f->w ) {
2447: fd = f->b;
2448: f->w = w; fd[0] |= 1; fd[d>>5] |= (1<<(d&31));
2449: i0 = j0 = k0 = 0;
2450: }
2451: /* if j0 > 0 then f is already a pentanomial */
2452: if ( j0 > 0 ) goto PENTA;
2453:
2454: /* searching for an irreducible trinomial */
2455:
2456: for ( i = 1; 2*i <= d; i++ ) {
2457: /* skip the polynomials 'before' f */
2458: if ( i < i0 ) continue;
2459: if ( i == i0 ) { i0 = 0; continue; }
2460: /* set i-th bit */
2461: fd[i>>5] |= (1<<(i&31));
2462: ret = irredcheck_dddup2(f);
2463: if ( ret == 1 ) return 0;
2464: /* reset i-th bit */
2465: fd[i>>5] &= ~(1<<(i&31));
2466: }
2467:
2468: /* searching for an irreducible pentanomial */
2469: PENTA:
2470: for ( i = 1; i < d; i++ ) {
2471: /* skip the polynomials 'before' f */
2472: if ( i < i0 ) continue;
2473: if ( i == i0 ) i0 = 0;
2474: /* set i-th bit */
2475: fd[i>>5] |= (1<<(i&31));
2476: for ( j = i+1; j < d; j++ ) {
2477: /* skip the polynomials 'before' f */
2478: if ( j < j0 ) continue;
2479: if ( j == j0 ) j0 = 0;
2480: /* set j-th bit */
2481: fd[j>>5] |= (1<<(j&31));
2482: for ( k = j+1; k < d; k++ ) {
2483: /* skip the polynomials 'before' f */
2484: if ( k < k0 ) continue;
2485: else if ( k == k0 ) { k0 = 0; continue; }
2486: /* set k-th bit */
2487: fd[k>>5] |= (1<<(k&31));
2488: ret = irredcheck_dddup2(f);
2489: if ( ret == 1 ) return 0;
2490: /* reset k-th bit */
2491: fd[k>>5] &= ~(1<<(k&31));
2492: }
2493: /* reset j-th bit */
2494: fd[j>>5] &= ~(1<<(j&31));
2495: }
2496: /* reset i-th bit */
2497: fd[i>>5] &= ~(1<<(i&31));
2498: }
2499: /* exhausted */
2500: return 1;
2501: }
2502:
2503: /*
2504: * f = an irreducible trinomial or pentanomial of degree d 'after' f
2505: *
2506: * searching strategy:
2507: * trinomial x^d+x^i+1:
2508: * i is as small as possible.
2509: * trinomial x^d+x^i+x^j+x^k+1:
2510: * i is as small as possible.
2511: * For such i, j is as small as possible.
2512: * For such i and j, 'k' is as small as possible.
2513: *
2514: * return value : 0 --- exists
2515: * 1 --- does not exist (exhaustion)
2516: */
2517:
2518: int _generate_good_irreducible_polynomial(UP2 f,int d)
2519: {
2520: int ret,i,j,k,nz,i0,j0,k0;
2521: int w;
2522: unsigned int *fd;
2523:
2524: /*
2525: * if f = x^d+x^i+1 then i0 <- i, j0 <- 0, k0 <-0.
2526: * if f = x^d+x^k+x^j+x^i+1 (k>j>i) then i0 <- i, j0 <- j, k0 <-k.
2527: * otherwise i0,j0,k0 is set to 0.
2528: */
2529:
2530: fd = f->b;
2531: w = (d>>5)+1;
2532: if ( f->w && (d==degup2(f)) ) {
2533: for ( nz = 0, i = d; i >= 0; i-- )
2534: if ( fd[i>>5]&(1<<(i&31)) ) nz++;
2535: switch ( nz ) {
2536: case 3:
2537: for ( i0 = 1; !(fd[i0>>5]&(1<<(i0&31))) ; i0++ );
2538: /* reset i0-th bit */
2539: fd[i0>>5] &= ~(1<<(i0&31));
2540: j0 = k0 = 0;
2541: break;
2542: case 5:
2543: for ( i0 = 1; !(fd[i0>>5]&(1<<(i0&31))) ; i0++ );
2544: /* reset i0-th bit */
2545: fd[i0>>5] &= ~(1<<(i0&31));
2546: for ( j0 = i0+1; !(fd[j0>>5]&(1<<(j0&31))) ; j0++ );
2547: /* reset j0-th bit */
2548: fd[j0>>5] &= ~(1<<(j0&31));
2549: for ( k0 = j0+1; !(fd[k0>>5]&(1<<(k0&31))) ; k0++ );
2550: /* reset k0-th bit */
2551: fd[k0>>5] &= ~(1<<(k0&31));
2552: break;
2553: default:
2554: f->w = 0; break;
2555: }
2556: } else
2557: f->w = 0;
2558:
2559: if ( !f->w ) {
2560: fd = f->b;
2561: f->w = w; fd[0] |= 1; fd[d>>5] |= (1<<(d&31));
2562: i0 = j0 = k0 = 0;
2563: }
2564: /* if j0 > 0 then f is already a pentanomial */
2565: if ( j0 > 0 ) goto PENTA;
2566:
2567: /* searching for an irreducible trinomial */
2568:
2569: for ( i = 1; 2*i <= d; i++ ) {
2570: /* skip the polynomials 'before' f */
2571: if ( i < i0 ) continue;
2572: if ( i == i0 ) { i0 = 0; continue; }
2573: /* set i-th bit */
2574: fd[i>>5] |= (1<<(i&31));
2575: ret = irredcheck_dddup2(f);
2576: if ( ret == 1 ) return 0;
2577: /* reset i-th bit */
2578: fd[i>>5] &= ~(1<<(i&31));
2579: }
2580:
2581: /* searching for an irreducible pentanomial */
2582: PENTA:
2583: for ( i = 3; i < d; i++ ) {
2584: /* skip the polynomials 'before' f */
2585: if ( i < i0 ) continue;
2586: if ( i == i0 ) i0 = 0;
2587: /* set i-th bit */
2588: fd[i>>5] |= (1<<(i&31));
2589: for ( j = 2; j < i; j++ ) {
2590: /* skip the polynomials 'before' f */
2591: if ( j < j0 ) continue;
2592: if ( j == j0 ) j0 = 0;
2593: /* set j-th bit */
2594: fd[j>>5] |= (1<<(j&31));
2595: for ( k = 1; k < j; k++ ) {
2596: /* skip the polynomials 'before' f */
2597: if ( k < k0 ) continue;
2598: else if ( k == k0 ) { k0 = 0; continue; }
2599: /* set k-th bit */
2600: fd[k>>5] |= (1<<(k&31));
2601: ret = irredcheck_dddup2(f);
2602: if ( ret == 1 ) return 0;
2603: /* reset k-th bit */
2604: fd[k>>5] &= ~(1<<(k&31));
2605: }
2606: /* reset j-th bit */
2607: fd[j>>5] &= ~(1<<(j&31));
2608: }
2609: /* reset i-th bit */
2610: fd[i>>5] &= ~(1<<(i&31));
2611: }
2612: /* exhausted */
2613: return 1;
1.3 noro 2614: }
2615:
2616: printqmat(mat,row,col)
2617: Q **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++ ) {
1.8 noro 2624: printnum((Num)mat[i][j]); printf(" ");
1.3 noro 2625: }
2626: printf("\n");
2627: }
2628: }
2629:
2630: printimat(mat,row,col)
2631: int **mat;
2632: int row,col;
2633: {
2634: int i,j;
2635:
2636: for ( i = 0; i < row; i++ ) {
2637: for ( j = 0; j < col; j++ ) {
2638: printf("%d ",mat[i][j]);
2639: }
2640: printf("\n");
2641: }
1.1 noro 2642: }