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