![[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.23 ! noro 48: * $OpenXM: OpenXM_contrib2/asir2000/builtin/array.c,v 1.22 2001/09/17 08:37:30 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.23 ! noro 77: void Pinvmat();
1.9 noro 78: void Pnewbytearray();
1.1 noro 79:
80: void Pgeneric_gauss_elim_mod();
81:
82: void Pmat_to_gfmmat(),Plu_gfmmat(),Psolve_by_lu_gfmmat();
83: void Pgeninvm_swap(), Premainder(), Psremainder(), Pvtol();
84: void sepvect();
85: void Pmulmat_gf2n();
86: void Pbconvmat_gf2n();
87: void Pmul_vect_mat_gf2n();
88: void PNBmul_gf2n();
89: void Pmul_mat_vect_int();
90: void Psepmat_destructive();
91: void Px962_irredpoly_up2();
92: void Pirredpoly_up2();
93: void Pnbpoly_up2();
94: void Pqsort();
1.14 noro 95: void Pexponent_vector();
1.1 noro 96:
97: struct ftab array_tab[] = {
98: {"solve_by_lu_gfmmat",Psolve_by_lu_gfmmat,4},
99: {"lu_gfmmat",Plu_gfmmat,2},
100: {"mat_to_gfmmat",Pmat_to_gfmmat,2},
101: {"generic_gauss_elim_mod",Pgeneric_gauss_elim_mod,2},
102: {"newvect",Pnewvect,-2},
1.14 noro 103: {"vector",Pnewvect,-2},
104: {"exponent_vector",Pexponent_vector,-99999999},
1.1 noro 105: {"newmat",Pnewmat,-3},
1.14 noro 106: {"matrix",Pnewmat,-3},
1.9 noro 107: {"newbytearray",Pnewbytearray,-2},
1.1 noro 108: {"sepmat_destructive",Psepmat_destructive,2},
109: {"sepvect",Psepvect,2},
110: {"qsort",Pqsort,-2},
111: {"vtol",Pvtol,1},
112: {"size",Psize,1},
113: {"det",Pdet,-2},
1.23 ! noro 114: {"invmat",Pinvmat,-2},
1.1 noro 115: {"leqm",Pleqm,2},
116: {"leqm1",Pleqm1,2},
117: {"geninvm",Pgeninvm,2},
118: {"geninvm_swap",Pgeninvm_swap,2},
119: {"remainder",Premainder,2},
120: {"sremainder",Psremainder,2},
121: {"mulmat_gf2n",Pmulmat_gf2n,1},
122: {"bconvmat_gf2n",Pbconvmat_gf2n,-4},
123: {"mul_vect_mat_gf2n",Pmul_vect_mat_gf2n,2},
124: {"mul_mat_vect_int",Pmul_mat_vect_int,2},
125: {"nbmul_gf2n",PNBmul_gf2n,3},
126: {"x962_irredpoly_up2",Px962_irredpoly_up2,2},
127: {"irredpoly_up2",Pirredpoly_up2,2},
128: {"nbpoly_up2",Pnbpoly_up2,2},
129: {0,0,0},
130: };
131:
132: int comp_obj(a,b)
133: Obj *a,*b;
134: {
135: return arf_comp(CO,*a,*b);
136: }
137:
138: static FUNC generic_comp_obj_func;
139: static NODE generic_comp_obj_arg;
140:
141: int generic_comp_obj(a,b)
142: Obj *a,*b;
143: {
144: Q r;
145:
146: BDY(generic_comp_obj_arg)=(pointer)(*a);
147: BDY(NEXT(generic_comp_obj_arg))=(pointer)(*b);
148: r = (Q)bevalf(generic_comp_obj_func,generic_comp_obj_arg);
149: if ( !r )
150: return 0;
151: else
152: return SGN(r)>0?1:-1;
153: }
154:
155:
156: void Pqsort(arg,rp)
157: NODE arg;
158: VECT *rp;
159: {
160: VECT vect;
161: char buf[BUFSIZ];
162: char *fname;
163: NODE n;
164: P p;
165: V v;
166:
167: asir_assert(ARG0(arg),O_VECT,"qsort");
168: vect = (VECT)ARG0(arg);
169: if ( argc(arg) == 1 )
170: qsort(BDY(vect),vect->len,sizeof(Obj),(int (*)(const void *,const void *))comp_obj);
171: else {
172: p = (P)ARG1(arg);
173: if ( !p || OID(p)!=2 )
174: error("qsort : invalid argument");
175: v = VR(p);
176: if ( (int)v->attr != V_SR )
177: error("qsort : no such function");
178: generic_comp_obj_func = (FUNC)v->priv;
179: MKNODE(n,0,0); MKNODE(generic_comp_obj_arg,0,n);
180: qsort(BDY(vect),vect->len,sizeof(Obj),(int (*)(const void *,const void *))generic_comp_obj);
181: }
182: *rp = vect;
183: }
184:
185: void PNBmul_gf2n(arg,rp)
186: NODE arg;
187: GF2N *rp;
188: {
189: GF2N a,b;
190: GF2MAT mat;
191: int n,w;
192: unsigned int *ab,*bb;
193: UP2 r;
194:
195: a = (GF2N)ARG0(arg);
196: b = (GF2N)ARG1(arg);
197: mat = (GF2MAT)ARG2(arg);
198: if ( !a || !b )
199: *rp = 0;
200: else {
201: n = mat->row;
202: w = (n+BSH-1)/BSH;
203:
204: ab = (unsigned int *)ALLOCA(w*sizeof(unsigned int));
205: bzero((char *)ab,w*sizeof(unsigned int));
206: bcopy(a->body->b,ab,(a->body->w)*sizeof(unsigned int));
207:
208: bb = (unsigned int *)ALLOCA(w*sizeof(unsigned int));
209: bzero((char *)bb,w*sizeof(unsigned int));
210: bcopy(b->body->b,bb,(b->body->w)*sizeof(unsigned int));
211:
212: NEWUP2(r,w);
213: bzero((char *)r->b,w*sizeof(unsigned int));
214: mul_nb(mat,ab,bb,r->b);
215: r->w = w;
216: _adjup2(r);
217: if ( !r->w )
218: *rp = 0;
219: else
220: MKGF2N(r,*rp);
221: }
222: }
223:
224: void Pmul_vect_mat_gf2n(arg,rp)
225: NODE arg;
226: GF2N *rp;
227: {
228: GF2N a;
229: GF2MAT mat;
230: int n,w;
231: unsigned int *b;
232: UP2 r;
233:
234: a = (GF2N)ARG0(arg);
235: mat = (GF2MAT)ARG1(arg);
236: if ( !a )
237: *rp = 0;
238: else {
239: n = mat->row;
240: w = (n+BSH-1)/BSH;
241: b = (unsigned int *)ALLOCA(w*sizeof(unsigned int));
242: bzero((char *)b,w*sizeof(unsigned int));
243: bcopy(a->body->b,b,(a->body->w)*sizeof(unsigned int));
244: NEWUP2(r,w);
245: bzero((char *)r->b,w*sizeof(unsigned int));
246: mulgf2vectmat(mat->row,b,mat->body,r->b);
247: r->w = w;
248: _adjup2(r);
249: if ( !r->w )
250: *rp = 0;
251: else {
252: MKGF2N(r,*rp);
253: }
254: }
255: }
256:
257: void Pbconvmat_gf2n(arg,rp)
258: NODE arg;
259: LIST *rp;
260: {
261: P p0,p1;
262: int to;
263: GF2MAT p01,p10;
264: GF2N root;
265: NODE n0,n1;
266:
267: p0 = (P)ARG0(arg);
268: p1 = (P)ARG1(arg);
269: to = ARG2(arg)?1:0;
270: if ( argc(arg) == 4 ) {
271: root = (GF2N)ARG3(arg);
272: compute_change_of_basis_matrix_with_root(p0,p1,to,root,&p01,&p10);
273: } else
274: compute_change_of_basis_matrix(p0,p1,to,&p01,&p10);
275: MKNODE(n1,p10,0); MKNODE(n0,p01,n1);
276: MKLIST(*rp,n0);
277: }
278:
279: void Pmulmat_gf2n(arg,rp)
280: NODE arg;
281: GF2MAT *rp;
282: {
283: GF2MAT m;
284:
285: if ( !compute_multiplication_matrix((P)ARG0(arg),&m) )
286: error("mulmat_gf2n : input is not a normal polynomial");
287: *rp = m;
288: }
289:
290: void Psepmat_destructive(arg,rp)
291: NODE arg;
292: LIST *rp;
293: {
294: MAT mat,mat1;
295: int i,j,row,col;
296: Q **a,**a1;
297: Q ent;
298: N nm,mod,rem,quo;
299: int sgn;
300: NODE n0,n1;
301:
302: mat = (MAT)ARG0(arg); mod = NM((Q)ARG1(arg));
303: row = mat->row; col = mat->col;
304: MKMAT(mat1,row,col);
305: a = (Q **)mat->body; a1 = (Q **)mat1->body;
306: for ( i = 0; i < row; i++ )
307: for ( j = 0; j < col; j++ ) {
308: ent = a[i][j];
309: if ( !ent )
310: continue;
311: nm = NM(ent);
312: sgn = SGN(ent);
313: divn(nm,mod,&quo,&rem);
314: /* if ( quo != nm && rem != nm ) */
315: /* GC_free(nm); */
316: /* GC_free(ent); */
317: NTOQ(rem,sgn,a[i][j]); NTOQ(quo,sgn,a1[i][j]);
318: }
319: MKNODE(n1,mat1,0); MKNODE(n0,mat,n1);
320: MKLIST(*rp,n0);
321: }
322:
323: void Psepvect(arg,rp)
324: NODE arg;
325: VECT *rp;
326: {
327: sepvect((VECT)ARG0(arg),QTOS((Q)ARG1(arg)),rp);
328: }
329:
330: void sepvect(v,d,rp)
331: VECT v;
332: int d;
333: VECT *rp;
334: {
335: int i,j,k,n,q,q1,r;
336: pointer *pv,*pw,*pu;
337: VECT w,u;
338:
339: n = v->len;
340: if ( d > n )
341: d = n;
342: q = n/d; r = n%d; q1 = q+1;
343: MKVECT(w,d); *rp = w;
344: pv = BDY(v); pw = BDY(w); k = 0;
345: for ( i = 0; i < r; i++ ) {
346: MKVECT(u,q1); pw[i] = (pointer)u;
347: for ( pu = BDY(u), j = 0; j < q1; j++, k++ )
348: pu[j] = pv[k];
349: }
350: for ( ; i < d; i++ ) {
351: MKVECT(u,q); pw[i] = (pointer)u;
352: for ( pu = BDY(u), j = 0; j < q; j++, k++ )
353: pu[j] = pv[k];
354: }
355: }
356:
357: void Pnewvect(arg,rp)
358: NODE arg;
359: VECT *rp;
360: {
361: int len,i,r;
362: VECT vect;
363: pointer *vb;
364: LIST list;
365: NODE tn;
366:
367: asir_assert(ARG0(arg),O_N,"newvect");
368: len = QTOS((Q)ARG0(arg));
1.5 noro 369: if ( len < 0 )
1.1 noro 370: error("newvect : invalid size");
371: MKVECT(vect,len);
372: if ( argc(arg) == 2 ) {
373: list = (LIST)ARG1(arg);
374: asir_assert(list,O_LIST,"newvect");
375: for ( r = 0, tn = BDY(list); tn; r++, tn = NEXT(tn) );
376: if ( r > len ) {
377: *rp = vect;
378: return;
379: }
380: for ( i = 0, tn = BDY(list), vb = BDY(vect); tn; i++, tn = NEXT(tn) )
381: vb[i] = (pointer)BDY(tn);
382: }
383: *rp = vect;
1.14 noro 384: }
385:
386: void Pexponent_vector(arg,rp)
387: NODE arg;
388: DP *rp;
389: {
390: nodetod(arg,rp);
1.9 noro 391: }
392:
393: void Pnewbytearray(arg,rp)
394: NODE arg;
395: BYTEARRAY *rp;
396: {
397: int len,i,r;
398: BYTEARRAY array;
399: unsigned char *vb;
1.10 noro 400: char *str;
1.9 noro 401: LIST list;
402: NODE tn;
403:
404: asir_assert(ARG0(arg),O_N,"newbytearray");
405: len = QTOS((Q)ARG0(arg));
406: if ( len < 0 )
407: error("newbytearray : invalid size");
408: MKBYTEARRAY(array,len);
409: if ( argc(arg) == 2 ) {
1.10 noro 410: if ( !ARG1(arg) )
411: error("newbytearray : invalid initialization");
412: switch ( OID((Obj)ARG1(arg)) ) {
413: case O_LIST:
414: list = (LIST)ARG1(arg);
415: asir_assert(list,O_LIST,"newbytearray");
416: for ( r = 0, tn = BDY(list); tn; r++, tn = NEXT(tn) );
417: if ( r <= len ) {
418: for ( i = 0, tn = BDY(list), vb = BDY(array); tn;
419: i++, tn = NEXT(tn) )
420: vb[i] = (unsigned char)QTOS((Q)BDY(tn));
421: }
422: break;
423: case O_STR:
424: str = BDY((STRING)ARG1(arg));
425: r = strlen(str);
426: if ( r <= len )
427: bcopy(str,BDY(array),r);
428: break;
429: default:
430: if ( !ARG1(arg) )
431: error("newbytearray : invalid initialization");
1.9 noro 432: }
433: }
434: *rp = array;
1.1 noro 435: }
436:
437: void Pnewmat(arg,rp)
438: NODE arg;
439: MAT *rp;
440: {
441: int row,col;
442: int i,j,r,c;
443: NODE tn,sn;
444: MAT m;
445: pointer **mb;
446: LIST list;
447:
448: asir_assert(ARG0(arg),O_N,"newmat");
449: asir_assert(ARG1(arg),O_N,"newmat");
450: row = QTOS((Q)ARG0(arg)); col = QTOS((Q)ARG1(arg));
1.5 noro 451: if ( row < 0 || col < 0 )
1.1 noro 452: error("newmat : invalid size");
453: MKMAT(m,row,col);
454: if ( argc(arg) == 3 ) {
455: list = (LIST)ARG2(arg);
456: asir_assert(list,O_LIST,"newmat");
457: for ( r = 0, c = 0, tn = BDY(list); tn; r++, tn = NEXT(tn) ) {
458: for ( j = 0, sn = BDY((LIST)BDY(tn)); sn; j++, sn = NEXT(sn) );
459: c = MAX(c,j);
460: }
461: if ( (r > row) || (c > col) ) {
462: *rp = m;
463: return;
464: }
465: for ( i = 0, tn = BDY(list), mb = BDY(m); tn; i++, tn = NEXT(tn) ) {
466: asir_assert(BDY(tn),O_LIST,"newmat");
467: for ( j = 0, sn = BDY((LIST)BDY(tn)); sn; j++, sn = NEXT(sn) )
468: mb[i][j] = (pointer)BDY(sn);
469: }
470: }
471: *rp = m;
472: }
473:
474: void Pvtol(arg,rp)
475: NODE arg;
476: LIST *rp;
477: {
478: NODE n,n1;
479: VECT v;
480: pointer *a;
481: int len,i;
482:
483: asir_assert(ARG0(arg),O_VECT,"vtol");
484: v = (VECT)ARG0(arg); len = v->len; a = BDY(v);
485: for ( i = len - 1, n = 0; i >= 0; i-- ) {
486: MKNODE(n1,a[i],n); n = n1;
487: }
488: MKLIST(*rp,n);
489: }
490:
491: void Premainder(arg,rp)
492: NODE arg;
493: Obj *rp;
494: {
495: Obj a;
496: VECT v,w;
497: MAT m,l;
498: pointer *vb,*wb;
499: pointer **mb,**lb;
500: int id,i,j,n,row,col,t,smd,sgn;
501: Q md,q;
502:
503: a = (Obj)ARG0(arg); md = (Q)ARG1(arg);
504: if ( !a )
505: *rp = 0;
506: else {
507: id = OID(a);
508: switch ( id ) {
509: case O_N:
510: case O_P:
511: cmp(md,(P)a,(P *)rp); break;
512: case O_VECT:
513: smd = QTOS(md);
514: v = (VECT)a; n = v->len; vb = v->body;
515: MKVECT(w,n); wb = w->body;
516: for ( i = 0; i < n; i++ ) {
517: if ( q = (Q)vb[i] ) {
518: sgn = SGN(q); t = rem(NM(q),smd);
519: STOQ(t,q);
520: if ( q )
521: SGN(q) = sgn;
522: }
523: wb[i] = (pointer)q;
524: }
525: *rp = (Obj)w;
526: break;
527: case O_MAT:
528: m = (MAT)a; row = m->row; col = m->col; mb = m->body;
529: MKMAT(l,row,col); lb = l->body;
530: for ( i = 0; i < row; i++ )
531: for ( j = 0, vb = mb[i], wb = lb[i]; j < col; j++ )
532: cmp(md,(P)vb[j],(P *)&wb[j]);
533: *rp = (Obj)l;
534: break;
535: default:
536: error("remainder : invalid argument");
537: }
538: }
539: }
540:
541: void Psremainder(arg,rp)
542: NODE arg;
543: Obj *rp;
544: {
545: Obj a;
546: VECT v,w;
547: MAT m,l;
548: pointer *vb,*wb;
549: pointer **mb,**lb;
550: unsigned int t,smd;
551: int id,i,j,n,row,col;
552: Q md,q;
553:
554: a = (Obj)ARG0(arg); md = (Q)ARG1(arg);
555: if ( !a )
556: *rp = 0;
557: else {
558: id = OID(a);
559: switch ( id ) {
560: case O_N:
561: case O_P:
562: cmp(md,(P)a,(P *)rp); break;
563: case O_VECT:
564: smd = QTOS(md);
565: v = (VECT)a; n = v->len; vb = v->body;
566: MKVECT(w,n); wb = w->body;
567: for ( i = 0; i < n; i++ ) {
568: if ( q = (Q)vb[i] ) {
569: t = (unsigned int)rem(NM(q),smd);
570: if ( SGN(q) < 0 )
571: t = (smd - t) % smd;
572: UTOQ(t,q);
573: }
574: wb[i] = (pointer)q;
575: }
576: *rp = (Obj)w;
577: break;
578: case O_MAT:
579: m = (MAT)a; row = m->row; col = m->col; mb = m->body;
580: MKMAT(l,row,col); lb = l->body;
581: for ( i = 0; i < row; i++ )
582: for ( j = 0, vb = mb[i], wb = lb[i]; j < col; j++ )
583: cmp(md,(P)vb[j],(P *)&wb[j]);
584: *rp = (Obj)l;
585: break;
586: default:
587: error("remainder : invalid argument");
588: }
589: }
590: }
591:
592: void Psize(arg,rp)
593: NODE arg;
594: LIST *rp;
595: {
596:
597: int n,m;
598: Q q;
599: NODE t,s;
600:
601: if ( !ARG0(arg) )
602: t = 0;
603: else {
604: switch (OID(ARG0(arg))) {
605: case O_VECT:
606: n = ((VECT)ARG0(arg))->len;
607: STOQ(n,q); MKNODE(t,q,0);
608: break;
609: case O_MAT:
610: n = ((MAT)ARG0(arg))->row; m = ((MAT)ARG0(arg))->col;
611: STOQ(m,q); MKNODE(s,q,0); STOQ(n,q); MKNODE(t,q,s);
612: break;
613: default:
614: error("size : invalid argument"); break;
615: }
616: }
617: MKLIST(*rp,t);
618: }
619:
620: void Pdet(arg,rp)
621: NODE arg;
622: P *rp;
623: {
624: MAT m;
625: int n,i,j,mod;
626: P d;
627: P **mat,**w;
628:
629: m = (MAT)ARG0(arg);
630: asir_assert(m,O_MAT,"det");
631: if ( m->row != m->col )
632: error("det : non-square matrix");
633: else if ( argc(arg) == 1 )
634: detp(CO,(P **)BDY(m),m->row,rp);
635: else {
636: n = m->row; mod = QTOS((Q)ARG1(arg)); mat = (P **)BDY(m);
637: w = (P **)almat_pointer(n,n);
638: for ( i = 0; i < n; i++ )
639: for ( j = 0; j < n; j++ )
640: ptomp(mod,mat[i][j],&w[i][j]);
641: detmp(CO,mod,w,n,&d);
642: mptop(d,rp);
1.23 ! noro 643: }
! 644: }
! 645:
! 646: void Pinvmat(arg,rp)
! 647: NODE arg;
! 648: LIST *rp;
! 649: {
! 650: MAT m,r;
! 651: int n,i,j,mod;
! 652: P dn;
! 653: P **mat,**imat,**w;
! 654: NODE nd;
! 655:
! 656: m = (MAT)ARG0(arg);
! 657: asir_assert(m,O_MAT,"invmat");
! 658: if ( m->row != m->col )
! 659: error("invmat : non-square matrix");
! 660: else if ( argc(arg) == 1 ) {
! 661: n = m->row;
! 662: invmatp(CO,(P **)BDY(m),n,&imat,&dn);
! 663: NEWMAT(r); r->row = n; r->col = n; r->body = (pointer **)imat;
! 664: nd = mknode(2,r,dn);
! 665: MKLIST(*rp,nd);
! 666: } else {
! 667: n = m->row; mod = QTOS((Q)ARG1(arg)); mat = (P **)BDY(m);
! 668: w = (P **)almat_pointer(n,n);
! 669: for ( i = 0; i < n; i++ )
! 670: for ( j = 0; j < n; j++ )
! 671: ptomp(mod,mat[i][j],&w[i][j]);
! 672: #if 0
! 673: detmp(CO,mod,w,n,&d);
! 674: mptop(d,rp);
! 675: #else
! 676: error("not implemented yet");
! 677: #endif
1.1 noro 678: }
679: }
680:
681: /*
682: input : a row x col matrix A
683: A[I] <-> A[I][0]*x_0+A[I][1]*x_1+...
684:
685: output : [B,R,C]
686: B : a rank(A) x col-rank(A) matrix
687: R : a vector of length rank(A)
688: C : a vector of length col-rank(A)
689: B[I] <-> x_{R[I]}+B[I][0]x_{C[0]}+B[I][1]x_{C[1]}+...
690: */
691:
692: void Pgeneric_gauss_elim_mod(arg,rp)
693: NODE arg;
694: LIST *rp;
695: {
696: NODE n0;
697: MAT m,mat;
698: VECT rind,cind;
699: Q **tmat;
700: int **wmat;
701: Q *rib,*cib;
702: int *colstat;
703: Q q;
704: int md,i,j,k,l,row,col,t,n,rank;
705:
706: asir_assert(ARG0(arg),O_MAT,"generic_gauss_elim_mod");
707: asir_assert(ARG1(arg),O_N,"generic_gauss_elim_mod");
708: m = (MAT)ARG0(arg); md = QTOS((Q)ARG1(arg));
709: row = m->row; col = m->col; tmat = (Q **)m->body;
710: wmat = (int **)almat(row,col);
711: colstat = (int *)MALLOC_ATOMIC(col*sizeof(int));
712: for ( i = 0; i < row; i++ )
713: for ( j = 0; j < col; j++ )
714: if ( q = (Q)tmat[i][j] ) {
715: t = rem(NM(q),md);
716: if ( t && SGN(q) < 0 )
717: t = (md - t) % md;
718: wmat[i][j] = t;
719: } else
720: wmat[i][j] = 0;
721: rank = generic_gauss_elim_mod(wmat,row,col,md,colstat);
722:
723: MKMAT(mat,rank,col-rank);
724: tmat = (Q **)mat->body;
725: for ( i = 0; i < rank; i++ )
726: for ( j = k = 0; j < col; j++ )
727: if ( !colstat[j] ) {
728: UTOQ(wmat[i][j],tmat[i][k]); k++;
729: }
730:
731: MKVECT(rind,rank);
732: MKVECT(cind,col-rank);
733: rib = (Q *)rind->body; cib = (Q *)cind->body;
734: for ( j = k = l = 0; j < col; j++ )
735: if ( colstat[j] ) {
736: STOQ(j,rib[k]); k++;
737: } else {
738: STOQ(j,cib[l]); l++;
739: }
740: n0 = mknode(3,mat,rind,cind);
741: MKLIST(*rp,n0);
742: }
743:
744: void Pleqm(arg,rp)
745: NODE arg;
746: VECT *rp;
747: {
748: MAT m;
749: VECT vect;
750: pointer **mat;
751: Q *v;
752: Q q;
753: int **wmat;
754: int md,i,j,row,col,t,n,status;
755:
756: asir_assert(ARG0(arg),O_MAT,"leqm");
757: asir_assert(ARG1(arg),O_N,"leqm");
758: m = (MAT)ARG0(arg); md = QTOS((Q)ARG1(arg));
759: row = m->row; col = m->col; mat = m->body;
760: wmat = (int **)almat(row,col);
761: for ( i = 0; i < row; i++ )
762: for ( j = 0; j < col; j++ )
763: if ( q = (Q)mat[i][j] ) {
764: t = rem(NM(q),md);
765: if ( SGN(q) < 0 )
766: t = (md - t) % md;
767: wmat[i][j] = t;
768: } else
769: wmat[i][j] = 0;
770: status = gauss_elim_mod(wmat,row,col,md);
771: if ( status < 0 )
772: *rp = 0;
773: else if ( status > 0 )
774: *rp = (VECT)ONE;
775: else {
776: n = col - 1;
777: MKVECT(vect,n);
778: for ( i = 0, v = (Q *)vect->body; i < n; i++ ) {
779: t = (md-wmat[i][n])%md; STOQ(t,v[i]);
780: }
781: *rp = vect;
782: }
783: }
784:
785: int gauss_elim_mod(mat,row,col,md)
786: int **mat;
787: int row,col,md;
788: {
789: int i,j,k,inv,a,n;
790: int *t,*pivot;
791:
792: n = col - 1;
793: for ( j = 0; j < n; j++ ) {
794: for ( i = j; i < row && !mat[i][j]; i++ );
795: if ( i == row )
796: return 1;
797: if ( i != j ) {
798: t = mat[i]; mat[i] = mat[j]; mat[j] = t;
799: }
800: pivot = mat[j];
801: inv = invm(pivot[j],md);
802: for ( k = j; k <= n; k++ ) {
803: /* pivot[k] = dmar(pivot[k],inv,0,md); */
804: DMAR(pivot[k],inv,0,md,pivot[k])
805: }
806: for ( i = 0; i < row; i++ ) {
807: t = mat[i];
808: if ( i != j && (a = t[j]) )
809: for ( k = j, a = md - a; k <= n; k++ ) {
1.8 noro 810: unsigned int tk;
1.1 noro 811: /* t[k] = dmar(pivot[k],a,t[k],md); */
1.8 noro 812: DMAR(pivot[k],a,t[k],md,tk)
813: t[k] = tk;
1.1 noro 814: }
815: }
816: }
817: for ( i = n; i < row && !mat[i][n]; i++ );
818: if ( i == row )
819: return 0;
820: else
821: return -1;
822: }
823:
1.4 noro 824: struct oEGT eg_mod,eg_elim,eg_elim1,eg_elim2,eg_chrem,eg_gschk,eg_intrat,eg_symb;
1.1 noro 825:
826: int generic_gauss_elim(mat,nm,dn,rindp,cindp)
827: MAT mat;
828: MAT *nm;
829: Q *dn;
830: int **rindp,**cindp;
831: {
832: int **wmat;
833: Q **bmat;
834: N **tmat;
835: Q *bmi;
836: N *tmi;
837: Q q;
838: int *wmi;
839: int *colstat,*wcolstat,*rind,*cind;
840: int row,col,ind,md,i,j,k,l,t,t1,rank,rank0,inv;
841: N m1,m2,m3,s,u;
842: MAT r,crmat;
843: struct oEGT tmp0,tmp1;
844: struct oEGT eg_mod_split,eg_elim_split,eg_chrem_split;
845: struct oEGT eg_intrat_split,eg_gschk_split;
846: int ret;
847:
848: init_eg(&eg_mod_split); init_eg(&eg_chrem_split);
849: init_eg(&eg_elim_split); init_eg(&eg_intrat_split);
850: init_eg(&eg_gschk_split);
851: bmat = (Q **)mat->body;
852: row = mat->row; col = mat->col;
853: wmat = (int **)almat(row,col);
854: colstat = (int *)MALLOC_ATOMIC(col*sizeof(int));
855: wcolstat = (int *)MALLOC_ATOMIC(col*sizeof(int));
856: for ( ind = 0; ; ind++ ) {
1.11 noro 857: if ( DP_Print ) {
1.2 noro 858: fprintf(asir_out,"."); fflush(asir_out);
859: }
1.12 noro 860: md = get_lprime(ind);
1.1 noro 861: get_eg(&tmp0);
862: for ( i = 0; i < row; i++ )
863: for ( j = 0, bmi = bmat[i], wmi = wmat[i]; j < col; j++ )
864: if ( q = (Q)bmi[j] ) {
865: t = rem(NM(q),md);
866: if ( t && SGN(q) < 0 )
867: t = (md - t) % md;
868: wmi[j] = t;
869: } else
870: wmi[j] = 0;
871: get_eg(&tmp1);
872: add_eg(&eg_mod,&tmp0,&tmp1);
873: add_eg(&eg_mod_split,&tmp0,&tmp1);
874: get_eg(&tmp0);
875: rank = generic_gauss_elim_mod(wmat,row,col,md,wcolstat);
876: get_eg(&tmp1);
877: add_eg(&eg_elim,&tmp0,&tmp1);
878: add_eg(&eg_elim_split,&tmp0,&tmp1);
879: if ( !ind ) {
880: RESET:
881: UTON(md,m1);
882: rank0 = rank;
883: bcopy(wcolstat,colstat,col*sizeof(int));
884: MKMAT(crmat,rank,col-rank);
885: MKMAT(r,rank,col-rank); *nm = r;
886: tmat = (N **)crmat->body;
887: for ( i = 0; i < rank; i++ )
888: for ( j = k = 0, tmi = tmat[i], wmi = wmat[i]; j < col; j++ )
889: if ( !colstat[j] ) {
890: UTON(wmi[j],tmi[k]); k++;
891: }
892: } else {
893: if ( rank < rank0 ) {
1.11 noro 894: if ( DP_Print ) {
1.1 noro 895: fprintf(asir_out,"lower rank matrix; continuing...\n");
1.2 noro 896: fflush(asir_out);
897: }
1.1 noro 898: continue;
899: } else if ( rank > rank0 ) {
1.11 noro 900: if ( DP_Print ) {
1.1 noro 901: fprintf(asir_out,"higher rank matrix; resetting...\n");
1.2 noro 902: fflush(asir_out);
903: }
1.1 noro 904: goto RESET;
905: } else {
906: for ( j = 0; (j<col) && (colstat[j]==wcolstat[j]); j++ );
907: if ( j < col ) {
1.11 noro 908: if ( DP_Print ) {
1.1 noro 909: fprintf(asir_out,"inconsitent colstat; resetting...\n");
1.2 noro 910: fflush(asir_out);
911: }
1.1 noro 912: goto RESET;
913: }
914: }
915:
916: get_eg(&tmp0);
917: inv = invm(rem(m1,md),md);
918: UTON(md,m2); muln(m1,m2,&m3);
919: for ( i = 0; i < rank; i++ )
920: for ( j = k = 0, tmi = tmat[i], wmi = wmat[i]; j < col; j++ )
921: if ( !colstat[j] ) {
922: if ( tmi[k] ) {
923: /* f3 = f1+m1*(m1 mod m2)^(-1)*(f2 - f1 mod m2) */
924: t = rem(tmi[k],md);
925: if ( wmi[j] >= t )
926: t = wmi[j]-t;
927: else
928: t = md-(t-wmi[j]);
929: DMAR(t,inv,0,md,t1)
930: UTON(t1,u);
931: muln(m1,u,&s);
932: addn(tmi[k],s,&u); tmi[k] = u;
933: } else if ( wmi[j] ) {
934: /* f3 = m1*(m1 mod m2)^(-1)*f2 */
935: DMAR(wmi[j],inv,0,md,t)
936: UTON(t,u);
937: muln(m1,u,&s); tmi[k] = s;
938: }
939: k++;
940: }
941: m1 = m3;
942: get_eg(&tmp1);
943: add_eg(&eg_chrem,&tmp0,&tmp1);
944: add_eg(&eg_chrem_split,&tmp0,&tmp1);
945:
946: get_eg(&tmp0);
1.13 noro 947: if ( ind % 16 )
948: ret = 0;
949: else
950: ret = intmtoratm(crmat,m1,*nm,dn);
1.1 noro 951: get_eg(&tmp1);
952: add_eg(&eg_intrat,&tmp0,&tmp1);
953: add_eg(&eg_intrat_split,&tmp0,&tmp1);
954: if ( ret ) {
955: *rindp = rind = (int *)MALLOC_ATOMIC(rank*sizeof(int));
956: *cindp = cind = (int *)MALLOC_ATOMIC((col-rank)*sizeof(int));
957: for ( j = k = l = 0; j < col; j++ )
958: if ( colstat[j] )
959: rind[k++] = j;
960: else
961: cind[l++] = j;
962: get_eg(&tmp0);
1.3 noro 963: if ( gensolve_check(mat,*nm,*dn,rind,cind) ) {
964: get_eg(&tmp1);
965: add_eg(&eg_gschk,&tmp0,&tmp1);
966: add_eg(&eg_gschk_split,&tmp0,&tmp1);
1.11 noro 967: if ( DP_Print ) {
1.3 noro 968: print_eg("Mod",&eg_mod_split);
969: print_eg("Elim",&eg_elim_split);
970: print_eg("ChRem",&eg_chrem_split);
971: print_eg("IntRat",&eg_intrat_split);
972: print_eg("Check",&eg_gschk_split);
973: fflush(asir_out);
974: }
975: return rank;
976: }
977: }
978: }
979: }
980: }
981:
982: int generic_gauss_elim_hensel(mat,nmmat,dn,rindp,cindp)
983: MAT mat;
984: MAT *nmmat;
985: Q *dn;
986: int **rindp,**cindp;
987: {
988: MAT bmat,xmat;
989: Q **a0,**a,**b,**x,**nm;
990: Q *ai,*bi,*xi;
991: int row,col;
992: int **w;
993: int *wi;
994: int **wc;
995: Q mdq,q,s,u;
996: N tn;
997: int ind,md,i,j,k,l,li,ri,rank;
998: unsigned int t;
999: int *cinfo,*rinfo;
1000: int *rind,*cind;
1001: int count;
1002: struct oEGT eg_mul,eg_inv,tmp0,tmp1;
1003:
1004: a0 = (Q **)mat->body;
1005: row = mat->row; col = mat->col;
1006: w = (int **)almat(row,col);
1007: for ( ind = 0; ; ind++ ) {
1.12 noro 1008: md = get_lprime(ind);
1.3 noro 1009: STOQ(md,mdq);
1010: for ( i = 0; i < row; i++ )
1011: for ( j = 0, ai = a0[i], wi = w[i]; j < col; j++ )
1012: if ( q = (Q)ai[j] ) {
1013: t = rem(NM(q),md);
1014: if ( t && SGN(q) < 0 )
1015: t = (md - t) % md;
1016: wi[j] = t;
1017: } else
1018: wi[j] = 0;
1019:
1020: rank = find_lhs_and_lu_mod(w,row,col,md,&rinfo,&cinfo);
1021: a = (Q **)almat_pointer(rank,rank); /* lhs mat */
1022: MKMAT(bmat,rank,col-rank); b = (Q **)bmat->body; /* lhs mat */
1023: for ( j = li = ri = 0; j < col; j++ )
1024: if ( cinfo[j] ) {
1025: /* the column is in lhs */
1026: for ( i = 0; i < rank; i++ ) {
1027: w[i][li] = w[i][j];
1028: a[i][li] = a0[rinfo[i]][j];
1029: }
1030: li++;
1031: } else {
1032: /* the column is in rhs */
1033: for ( i = 0; i < rank; i++ )
1034: b[i][ri] = a0[rinfo[i]][j];
1035: ri++;
1036: }
1037:
1038: /* solve Ax+B=0; A: rank x rank, B: rank x ri */
1039: MKMAT(xmat,rank,ri); x = (Q **)(xmat)->body;
1040: MKMAT(*nmmat,rank,ri); nm = (Q **)(*nmmat)->body;
1041: /* use the right part of w as work area */
1042: /* ri = col - rank */
1043: wc = (int **)almat(rank,ri);
1044: for ( i = 0; i < rank; i++ )
1045: wc[i] = w[i]+rank;
1046: *rindp = rind = (int *)MALLOC_ATOMIC(rank*sizeof(int));
1047: *cindp = cind = (int *)MALLOC_ATOMIC((ri)*sizeof(int));
1048:
1049: init_eg(&eg_mul); init_eg(&eg_inv);
1050: for ( q = ONE, count = 0; ; count++ ) {
1051: fprintf(stderr,".");
1052: /* wc = -b mod md */
1053: for ( i = 0; i < rank; i++ )
1054: for ( j = 0, bi = b[i], wi = wc[i]; j < ri; j++ )
1055: if ( u = (Q)bi[j] ) {
1056: t = rem(NM(u),md);
1057: if ( t && SGN(u) > 0 )
1058: t = (md - t) % md;
1059: wi[j] = t;
1060: } else
1061: wi[j] = 0;
1062: /* wc = A^(-1)wc; wc is normalized */
1063: get_eg(&tmp0);
1064: solve_by_lu_mod(w,rank,md,wc,ri);
1.1 noro 1065: get_eg(&tmp1);
1.3 noro 1066: add_eg(&eg_inv,&tmp0,&tmp1);
1067: /* x = x-q*wc */
1068: for ( i = 0; i < rank; i++ )
1069: for ( j = 0, xi = x[i], wi = wc[i]; j < ri; j++ ) {
1070: STOQ(wi[j],u); mulq(q,u,&s);
1071: subq(xi[j],s,&u); xi[j] = u;
1072: }
1073: get_eg(&tmp0);
1074: for ( i = 0; i < rank; i++ )
1075: for ( j = 0; j < ri; j++ ) {
1076: inner_product_mat_int_mod(a,wc,rank,i,j,&u);
1077: addq(b[i][j],u,&s);
1078: if ( s ) {
1079: t = divin(NM(s),md,&tn);
1080: if ( t )
1081: error("generic_gauss_elim_hensel:incosistent");
1082: NTOQ(tn,SGN(s),b[i][j]);
1083: } else
1084: b[i][j] = 0;
1085: }
1086: get_eg(&tmp1);
1087: add_eg(&eg_mul,&tmp0,&tmp1);
1088: /* q = q*md */
1089: mulq(q,mdq,&u); q = u;
1.13 noro 1090: if ( !(count % 16) && intmtoratm_q(xmat,NM(q),*nmmat,dn) ) {
1.3 noro 1091: for ( j = k = l = 0; j < col; j++ )
1092: if ( cinfo[j] )
1093: rind[k++] = j;
1094: else
1095: cind[l++] = j;
1096: if ( gensolve_check(mat,*nmmat,*dn,rind,cind) ) {
1097: fprintf(stderr,"\n");
1098: print_eg("INV",&eg_inv);
1099: print_eg("MUL",&eg_mul);
1100: fflush(asir_out);
1101: return rank;
1102: }
1.1 noro 1103: }
1104: }
1105: }
1106: }
1107:
1108: int f4_nocheck;
1109:
1110: int gensolve_check(mat,nm,dn,rind,cind)
1111: MAT mat,nm;
1112: Q dn;
1113: int *rind,*cind;
1114: {
1115: int row,col,rank,clen,i,j,k,l;
1116: Q s,t,u;
1117: Q *w;
1118: Q *mati,*nmk;
1119:
1120: if ( f4_nocheck )
1121: return 1;
1122: row = mat->row; col = mat->col;
1123: rank = nm->row; clen = nm->col;
1124: w = (Q *)MALLOC(clen*sizeof(Q));
1125: for ( i = 0; i < row; i++ ) {
1126: mati = (Q *)mat->body[i];
1127: #if 1
1128: bzero(w,clen*sizeof(Q));
1129: for ( k = 0; k < rank; k++ )
1130: for ( l = 0, nmk = (Q *)nm->body[k]; l < clen; l++ ) {
1131: mulq(mati[rind[k]],nmk[l],&t);
1132: addq(w[l],t,&s); w[l] = s;
1133: }
1134: for ( j = 0; j < clen; j++ ) {
1135: mulq(dn,mati[cind[j]],&t);
1136: if ( cmpq(w[j],t) )
1137: break;
1138: }
1139: #else
1140: for ( j = 0; j < clen; j++ ) {
1141: for ( k = 0, s = 0; k < rank; k++ ) {
1142: mulq(mati[rind[k]],nm->body[k][j],&t);
1143: addq(s,t,&u); s = u;
1144: }
1145: mulq(dn,mati[cind[j]],&t);
1146: if ( cmpq(s,t) )
1147: break;
1148: }
1149: #endif
1150: if ( j != clen )
1151: break;
1152: }
1153: if ( i != row )
1154: return 0;
1155: else
1156: return 1;
1157: }
1158:
1159: /* assuming 0 < c < m */
1160:
1161: int inttorat(c,m,b,sgnp,nmp,dnp)
1162: N c,m,b;
1163: int *sgnp;
1164: N *nmp,*dnp;
1165: {
1166: Q qq,t,u1,v1,r1,nm;
1167: N q,r,u2,v2,r2;
1168:
1169: u1 = 0; v1 = ONE; u2 = m; v2 = c;
1170: while ( cmpn(v2,b) >= 0 ) {
1171: divn(u2,v2,&q,&r2); u2 = v2; v2 = r2;
1172: NTOQ(q,1,qq); mulq(qq,v1,&t); subq(u1,t,&r1); u1 = v1; v1 = r1;
1173: }
1174: if ( cmpn(NM(v1),b) >= 0 )
1175: return 0;
1176: else {
1177: *nmp = v2;
1178: *dnp = NM(v1);
1179: *sgnp = SGN(v1);
1180: return 1;
1181: }
1182: }
1183:
1184: /* mat->body = N ** */
1185:
1186: int intmtoratm(mat,md,nm,dn)
1187: MAT mat;
1188: N md;
1189: MAT nm;
1190: Q *dn;
1191: {
1192: N t,s,b;
1193: Q bound,dn0,dn1,nm1,q,tq;
1194: int i,j,k,l,row,col;
1195: Q **rmat;
1196: N **tmat;
1197: N *tmi;
1198: Q *nmk;
1199: N u,unm,udn;
1200: int sgn,ret;
1201:
1.3 noro 1202: if ( UNIN(md) )
1203: return 0;
1.1 noro 1204: row = mat->row; col = mat->col;
1205: bshiftn(md,1,&t);
1206: isqrt(t,&s);
1207: bshiftn(s,64,&b);
1208: if ( !b )
1209: b = ONEN;
1210: dn0 = ONE;
1211: tmat = (N **)mat->body;
1212: rmat = (Q **)nm->body;
1213: for ( i = 0; i < row; i++ )
1214: for ( j = 0, tmi = tmat[i]; j < col; j++ )
1215: if ( tmi[j] ) {
1216: muln(tmi[j],NM(dn0),&s);
1217: remn(s,md,&u);
1218: ret = inttorat(u,md,b,&sgn,&unm,&udn);
1219: if ( !ret )
1220: return 0;
1221: else {
1222: NTOQ(unm,sgn,nm1);
1223: NTOQ(udn,1,dn1);
1224: if ( !UNIQ(dn1) ) {
1225: for ( k = 0; k < i; k++ )
1226: for ( l = 0, nmk = rmat[k]; l < col; l++ ) {
1227: mulq(nmk[l],dn1,&q); nmk[l] = q;
1228: }
1229: for ( l = 0, nmk = rmat[i]; l < j; l++ ) {
1230: mulq(nmk[l],dn1,&q); nmk[l] = q;
1231: }
1232: }
1233: rmat[i][j] = nm1;
1234: mulq(dn0,dn1,&q); dn0 = q;
1235: }
1236: }
1237: *dn = dn0;
1238: return 1;
1239: }
1240:
1.3 noro 1241: /* mat->body = Q ** */
1242:
1243: int intmtoratm_q(mat,md,nm,dn)
1244: MAT mat;
1245: N md;
1246: MAT nm;
1247: Q *dn;
1248: {
1249: N t,s,b;
1250: Q bound,dn0,dn1,nm1,q,tq;
1251: int i,j,k,l,row,col;
1252: Q **rmat;
1253: Q **tmat;
1254: Q *tmi;
1255: Q *nmk;
1256: N u,unm,udn;
1257: int sgn,ret;
1258:
1259: if ( UNIN(md) )
1260: return 0;
1261: row = mat->row; col = mat->col;
1262: bshiftn(md,1,&t);
1263: isqrt(t,&s);
1264: bshiftn(s,64,&b);
1265: if ( !b )
1266: b = ONEN;
1267: dn0 = ONE;
1268: tmat = (Q **)mat->body;
1269: rmat = (Q **)nm->body;
1270: for ( i = 0; i < row; i++ )
1271: for ( j = 0, tmi = tmat[i]; j < col; j++ )
1272: if ( tmi[j] ) {
1273: muln(NM(tmi[j]),NM(dn0),&s);
1274: remn(s,md,&u);
1275: ret = inttorat(u,md,b,&sgn,&unm,&udn);
1276: if ( !ret )
1277: return 0;
1278: else {
1279: if ( SGN(tmi[j])<0 )
1280: sgn = -sgn;
1281: NTOQ(unm,sgn,nm1);
1282: NTOQ(udn,1,dn1);
1283: if ( !UNIQ(dn1) ) {
1284: for ( k = 0; k < i; k++ )
1285: for ( l = 0, nmk = rmat[k]; l < col; l++ ) {
1286: mulq(nmk[l],dn1,&q); nmk[l] = q;
1287: }
1288: for ( l = 0, nmk = rmat[i]; l < j; l++ ) {
1289: mulq(nmk[l],dn1,&q); nmk[l] = q;
1290: }
1291: }
1292: rmat[i][j] = nm1;
1293: mulq(dn0,dn1,&q); dn0 = q;
1294: }
1295: }
1296: *dn = dn0;
1297: return 1;
1298: }
1299:
1.4 noro 1300: #define ONE_STEP1 if ( zzz = *s ) { DMAR(zzz,hc,*tj,md,*tj) } tj++; s++;
1301:
1302: void reduce_reducers_mod(mat,row,col,md)
1303: int **mat;
1304: int row,col;
1305: int md;
1306: {
1307: int i,j,k,l,hc,zzz;
1308: int *t,*s,*tj,*ind;
1309:
1310: /* reduce the reducers */
1311: ind = (int *)ALLOCA(row*sizeof(int));
1312: for ( i = 0; i < row; i++ ) {
1313: t = mat[i];
1314: for ( j = 0; j < col && !t[j]; j++ );
1315: /* register the position of the head term */
1316: ind[i] = j;
1317: for ( l = i-1; l >= 0; l-- ) {
1318: /* reduce mat[i] by mat[l] */
1319: if ( hc = t[ind[l]] ) {
1320: /* mat[i] = mat[i]-hc*mat[l] */
1321: j = ind[l];
1322: s = mat[l]+j;
1323: tj = t+j;
1324: hc = md-hc;
1325: k = col-j;
1326: for ( ; k >= 64; k -= 64 ) {
1327: ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
1328: ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
1329: ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
1330: ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
1331: ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
1332: ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
1333: ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
1334: ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
1335: ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
1336: ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
1337: ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
1338: ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
1339: ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
1340: ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
1341: ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
1342: ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
1343: }
1.16 noro 1344: for ( ; k > 0; k-- ) {
1.4 noro 1345: if ( zzz = *s ) { DMAR(zzz,hc,*tj,md,*tj) } tj++; s++;
1346: }
1347: }
1348: }
1349: }
1350: }
1351:
1352: /*
1353: mat[i] : reducers (i=0,...,nred-1)
1354: spolys (i=nred,...,row-1)
1355: mat[0] < mat[1] < ... < mat[nred-1] w.r.t the term order
1356: 1. reduce the reducers
1357: 2. reduce spolys by the reduced reducers
1358: */
1359:
1360: void pre_reduce_mod(mat,row,col,nred,md)
1361: int **mat;
1362: int row,col,nred;
1363: int md;
1364: {
1365: int i,j,k,l,hc,inv;
1366: int *t,*s,*tk,*ind;
1367:
1368: #if 1
1369: /* reduce the reducers */
1370: ind = (int *)ALLOCA(row*sizeof(int));
1371: for ( i = 0; i < nred; i++ ) {
1372: /* make mat[i] monic and mat[i] by mat[0],...,mat[i-1] */
1373: t = mat[i];
1374: for ( j = 0; j < col && !t[j]; j++ );
1375: /* register the position of the head term */
1376: ind[i] = j;
1377: inv = invm(t[j],md);
1378: for ( k = j; k < col; k++ )
1379: if ( t[k] )
1380: DMAR(t[k],inv,0,md,t[k])
1381: for ( l = i-1; l >= 0; l-- ) {
1382: /* reduce mat[i] by mat[l] */
1383: if ( hc = t[ind[l]] ) {
1384: /* mat[i] = mat[i]-hc*mat[l] */
1385: for ( k = ind[l], hc = md-hc, s = mat[l]+k, tk = t+k;
1386: k < col; k++, tk++, s++ )
1387: if ( *s )
1388: DMAR(*s,hc,*tk,md,*tk)
1389: }
1390: }
1391: }
1392: /* reduce the spolys */
1393: for ( i = nred; i < row; i++ ) {
1394: t = mat[i];
1395: for ( l = nred-1; l >= 0; l-- ) {
1396: /* reduce mat[i] by mat[l] */
1397: if ( hc = t[ind[l]] ) {
1398: /* mat[i] = mat[i]-hc*mat[l] */
1399: for ( k = ind[l], hc = md-hc, s = mat[l]+k, tk = t+k;
1400: k < col; k++, tk++, s++ )
1401: if ( *s )
1402: DMAR(*s,hc,*tk,md,*tk)
1403: }
1404: }
1405: }
1406: #endif
1407: }
1408: /*
1409: mat[i] : reducers (i=0,...,nred-1)
1410: mat[0] < mat[1] < ... < mat[nred-1] w.r.t the term order
1411: */
1412:
1413: void reduce_sp_by_red_mod(sp,redmat,ind,nred,col,md)
1414: int *sp,**redmat;
1415: int *ind;
1416: int nred,col;
1417: int md;
1418: {
1419: int i,j,k,hc,zzz;
1420: int *t,*s,*tj;
1421:
1422: /* reduce the spolys by redmat */
1423: for ( i = nred-1; i >= 0; i-- ) {
1424: /* reduce sp by redmat[i] */
1425: if ( hc = sp[ind[i]] ) {
1426: /* sp = sp-hc*redmat[i] */
1427: j = ind[i];
1428: hc = md-hc;
1429: s = redmat[i]+j;
1430: tj = sp+j;
1.16 noro 1431: for ( k = col-j; k > 0; k-- ) {
1.4 noro 1432: if ( zzz = *s ) { DMAR(zzz,hc,*tj,md,*tj) } tj++; s++;
1.15 noro 1433: }
1434: }
1.17 noro 1435: }
1436: }
1437:
1438: /*
1.15 noro 1439: mat[i] : compressed reducers (i=0,...,nred-1)
1440: mat[0] < mat[1] < ... < mat[nred-1] w.r.t the term order
1441: */
1442:
1.21 noro 1443: int red_by_compress(m,p,r,ri,hc,len)
1.18 noro 1444: int m;
1445: unsigned int *p;
1.21 noro 1446: register unsigned int *r;
1.22 noro 1447: register unsigned int *ri;
1.18 noro 1448: unsigned int hc;
1.19 noro 1449: register int len;
1.18 noro 1450: {
1.19 noro 1451: unsigned int up,lo;
1.18 noro 1452: unsigned int dmy;
1453: unsigned int *pj;
1454:
1.21 noro 1455: p[*ri] = 0; r++; ri++;
1456: for ( len--; len; len--, r++, ri++ ) {
1457: pj = p+ *ri;
1458: DMA(*r,hc,*pj,up,lo);
1.18 noro 1459: if ( up ) {
1460: DSAB(m,up,lo,dmy,*pj);
1461: } else
1462: *pj = lo;
1463: }
1464: }
1465:
1466: /* p -= hc*r */
1467:
1468: int red_by_vect(m,p,r,hc,len)
1469: int m;
1470: unsigned int *p,*r;
1471: unsigned int hc;
1472: int len;
1473: {
1474: register unsigned int up,lo;
1475: unsigned int dmy;
1476:
1477: *p++ = 0; r++; len--;
1478: for ( ; len; len--, r++, p++ )
1479: if ( *r ) {
1.20 noro 1480: DMA(*r,hc,*p,up,lo);
1.18 noro 1481: if ( up ) {
1482: DSAB(m,up,lo,dmy,*p);
1483: } else
1484: *p = lo;
1485: }
1486: }
1487:
1.21 noro 1488: extern unsigned int **psca;
1489:
1.15 noro 1490: void reduce_sp_by_red_mod_compress (sp,redmat,ind,nred,col,md)
1491: int *sp;
1492: CDP *redmat;
1493: int *ind;
1494: int nred,col;
1495: int md;
1496: {
1.18 noro 1497: int i,j,k,len;
1498: unsigned int *tj;
1.15 noro 1499: CDP ri;
1.18 noro 1500: unsigned int hc,up,lo,up1,lo1,c;
1501: unsigned int *usp;
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.21 noro 1513: red_by_compress(md,usp,psca[ri->psindex],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: }