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