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