![[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.7 ! noro 48: * $OpenXM: OpenXM_contrib2/asir2000/builtin/array.c,v 1.6 2000/08/21 08:31:18 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++ ) {
715: /* t[k] = dmar(pivot[k],a,t[k],md); */
716: DMAR(pivot[k],a,t[k],md,t[k])
717: }
718: }
719: }
720: for ( i = n; i < row && !mat[i][n]; i++ );
721: if ( i == row )
722: return 0;
723: else
724: return -1;
725: }
726:
1.4 noro 727: struct oEGT eg_mod,eg_elim,eg_elim1,eg_elim2,eg_chrem,eg_gschk,eg_intrat,eg_symb;
1.1 noro 728:
729: int generic_gauss_elim(mat,nm,dn,rindp,cindp)
730: MAT mat;
731: MAT *nm;
732: Q *dn;
733: int **rindp,**cindp;
734: {
735: int **wmat;
736: Q **bmat;
737: N **tmat;
738: Q *bmi;
739: N *tmi;
740: Q q;
741: int *wmi;
742: int *colstat,*wcolstat,*rind,*cind;
743: int row,col,ind,md,i,j,k,l,t,t1,rank,rank0,inv;
744: N m1,m2,m3,s,u;
745: MAT r,crmat;
746: struct oEGT tmp0,tmp1;
747: struct oEGT eg_mod_split,eg_elim_split,eg_chrem_split;
748: struct oEGT eg_intrat_split,eg_gschk_split;
749: int ret;
750:
751: init_eg(&eg_mod_split); init_eg(&eg_chrem_split);
752: init_eg(&eg_elim_split); init_eg(&eg_intrat_split);
753: init_eg(&eg_gschk_split);
754: bmat = (Q **)mat->body;
755: row = mat->row; col = mat->col;
756: wmat = (int **)almat(row,col);
757: colstat = (int *)MALLOC_ATOMIC(col*sizeof(int));
758: wcolstat = (int *)MALLOC_ATOMIC(col*sizeof(int));
759: for ( ind = 0; ; ind++ ) {
1.2 noro 760: if ( Print ) {
761: fprintf(asir_out,"."); fflush(asir_out);
762: }
1.1 noro 763: md = lprime[ind];
764: get_eg(&tmp0);
765: for ( i = 0; i < row; i++ )
766: for ( j = 0, bmi = bmat[i], wmi = wmat[i]; j < col; j++ )
767: if ( q = (Q)bmi[j] ) {
768: t = rem(NM(q),md);
769: if ( t && SGN(q) < 0 )
770: t = (md - t) % md;
771: wmi[j] = t;
772: } else
773: wmi[j] = 0;
774: get_eg(&tmp1);
775: add_eg(&eg_mod,&tmp0,&tmp1);
776: add_eg(&eg_mod_split,&tmp0,&tmp1);
777: get_eg(&tmp0);
778: rank = generic_gauss_elim_mod(wmat,row,col,md,wcolstat);
779: get_eg(&tmp1);
780: add_eg(&eg_elim,&tmp0,&tmp1);
781: add_eg(&eg_elim_split,&tmp0,&tmp1);
782: if ( !ind ) {
783: RESET:
784: UTON(md,m1);
785: rank0 = rank;
786: bcopy(wcolstat,colstat,col*sizeof(int));
787: MKMAT(crmat,rank,col-rank);
788: MKMAT(r,rank,col-rank); *nm = r;
789: tmat = (N **)crmat->body;
790: for ( i = 0; i < rank; i++ )
791: for ( j = k = 0, tmi = tmat[i], wmi = wmat[i]; j < col; j++ )
792: if ( !colstat[j] ) {
793: UTON(wmi[j],tmi[k]); k++;
794: }
795: } else {
796: if ( rank < rank0 ) {
1.2 noro 797: if ( Print ) {
1.1 noro 798: fprintf(asir_out,"lower rank matrix; continuing...\n");
1.2 noro 799: fflush(asir_out);
800: }
1.1 noro 801: continue;
802: } else if ( rank > rank0 ) {
1.2 noro 803: if ( Print ) {
1.1 noro 804: fprintf(asir_out,"higher rank matrix; resetting...\n");
1.2 noro 805: fflush(asir_out);
806: }
1.1 noro 807: goto RESET;
808: } else {
809: for ( j = 0; (j<col) && (colstat[j]==wcolstat[j]); j++ );
810: if ( j < col ) {
1.2 noro 811: if ( Print ) {
1.1 noro 812: fprintf(asir_out,"inconsitent colstat; resetting...\n");
1.2 noro 813: fflush(asir_out);
814: }
1.1 noro 815: goto RESET;
816: }
817: }
818:
819: get_eg(&tmp0);
820: inv = invm(rem(m1,md),md);
821: UTON(md,m2); muln(m1,m2,&m3);
822: for ( i = 0; i < rank; i++ )
823: for ( j = k = 0, tmi = tmat[i], wmi = wmat[i]; j < col; j++ )
824: if ( !colstat[j] ) {
825: if ( tmi[k] ) {
826: /* f3 = f1+m1*(m1 mod m2)^(-1)*(f2 - f1 mod m2) */
827: t = rem(tmi[k],md);
828: if ( wmi[j] >= t )
829: t = wmi[j]-t;
830: else
831: t = md-(t-wmi[j]);
832: DMAR(t,inv,0,md,t1)
833: UTON(t1,u);
834: muln(m1,u,&s);
835: addn(tmi[k],s,&u); tmi[k] = u;
836: } else if ( wmi[j] ) {
837: /* f3 = m1*(m1 mod m2)^(-1)*f2 */
838: DMAR(wmi[j],inv,0,md,t)
839: UTON(t,u);
840: muln(m1,u,&s); tmi[k] = s;
841: }
842: k++;
843: }
844: m1 = m3;
845: get_eg(&tmp1);
846: add_eg(&eg_chrem,&tmp0,&tmp1);
847: add_eg(&eg_chrem_split,&tmp0,&tmp1);
848:
849: get_eg(&tmp0);
850: ret = intmtoratm(crmat,m1,*nm,dn);
851: get_eg(&tmp1);
852: add_eg(&eg_intrat,&tmp0,&tmp1);
853: add_eg(&eg_intrat_split,&tmp0,&tmp1);
854: if ( ret ) {
855: *rindp = rind = (int *)MALLOC_ATOMIC(rank*sizeof(int));
856: *cindp = cind = (int *)MALLOC_ATOMIC((col-rank)*sizeof(int));
857: for ( j = k = l = 0; j < col; j++ )
858: if ( colstat[j] )
859: rind[k++] = j;
860: else
861: cind[l++] = j;
862: get_eg(&tmp0);
1.3 noro 863: if ( gensolve_check(mat,*nm,*dn,rind,cind) ) {
864: get_eg(&tmp1);
865: add_eg(&eg_gschk,&tmp0,&tmp1);
866: add_eg(&eg_gschk_split,&tmp0,&tmp1);
867: if ( Print ) {
868: print_eg("Mod",&eg_mod_split);
869: print_eg("Elim",&eg_elim_split);
870: print_eg("ChRem",&eg_chrem_split);
871: print_eg("IntRat",&eg_intrat_split);
872: print_eg("Check",&eg_gschk_split);
873: fflush(asir_out);
874: }
875: return rank;
876: }
877: }
878: }
879: }
880: }
881:
882: int generic_gauss_elim_hensel(mat,nmmat,dn,rindp,cindp)
883: MAT mat;
884: MAT *nmmat;
885: Q *dn;
886: int **rindp,**cindp;
887: {
888: MAT bmat,xmat;
889: Q **a0,**a,**b,**x,**nm;
890: Q *ai,*bi,*xi;
891: int row,col;
892: int **w;
893: int *wi;
894: int **wc;
895: Q mdq,q,s,u;
896: N tn;
897: int ind,md,i,j,k,l,li,ri,rank;
898: unsigned int t;
899: int *cinfo,*rinfo;
900: int *rind,*cind;
901: int count;
902: struct oEGT eg_mul,eg_inv,tmp0,tmp1;
903:
904: a0 = (Q **)mat->body;
905: row = mat->row; col = mat->col;
906: w = (int **)almat(row,col);
907: for ( ind = 0; ; ind++ ) {
908: md = lprime[ind];
909: STOQ(md,mdq);
910: for ( i = 0; i < row; i++ )
911: for ( j = 0, ai = a0[i], wi = w[i]; j < col; j++ )
912: if ( q = (Q)ai[j] ) {
913: t = rem(NM(q),md);
914: if ( t && SGN(q) < 0 )
915: t = (md - t) % md;
916: wi[j] = t;
917: } else
918: wi[j] = 0;
919:
920: rank = find_lhs_and_lu_mod(w,row,col,md,&rinfo,&cinfo);
921: a = (Q **)almat_pointer(rank,rank); /* lhs mat */
922: MKMAT(bmat,rank,col-rank); b = (Q **)bmat->body; /* lhs mat */
923: for ( j = li = ri = 0; j < col; j++ )
924: if ( cinfo[j] ) {
925: /* the column is in lhs */
926: for ( i = 0; i < rank; i++ ) {
927: w[i][li] = w[i][j];
928: a[i][li] = a0[rinfo[i]][j];
929: }
930: li++;
931: } else {
932: /* the column is in rhs */
933: for ( i = 0; i < rank; i++ )
934: b[i][ri] = a0[rinfo[i]][j];
935: ri++;
936: }
937:
938: /* solve Ax+B=0; A: rank x rank, B: rank x ri */
939: MKMAT(xmat,rank,ri); x = (Q **)(xmat)->body;
940: MKMAT(*nmmat,rank,ri); nm = (Q **)(*nmmat)->body;
941: /* use the right part of w as work area */
942: /* ri = col - rank */
943: wc = (int **)almat(rank,ri);
944: for ( i = 0; i < rank; i++ )
945: wc[i] = w[i]+rank;
946: *rindp = rind = (int *)MALLOC_ATOMIC(rank*sizeof(int));
947: *cindp = cind = (int *)MALLOC_ATOMIC((ri)*sizeof(int));
948:
949: init_eg(&eg_mul); init_eg(&eg_inv);
950: for ( q = ONE, count = 0; ; count++ ) {
951: fprintf(stderr,".");
952: /* wc = -b mod md */
953: for ( i = 0; i < rank; i++ )
954: for ( j = 0, bi = b[i], wi = wc[i]; j < ri; j++ )
955: if ( u = (Q)bi[j] ) {
956: t = rem(NM(u),md);
957: if ( t && SGN(u) > 0 )
958: t = (md - t) % md;
959: wi[j] = t;
960: } else
961: wi[j] = 0;
962: /* wc = A^(-1)wc; wc is normalized */
963: get_eg(&tmp0);
964: solve_by_lu_mod(w,rank,md,wc,ri);
1.1 noro 965: get_eg(&tmp1);
1.3 noro 966: add_eg(&eg_inv,&tmp0,&tmp1);
967: /* x = x-q*wc */
968: for ( i = 0; i < rank; i++ )
969: for ( j = 0, xi = x[i], wi = wc[i]; j < ri; j++ ) {
970: STOQ(wi[j],u); mulq(q,u,&s);
971: subq(xi[j],s,&u); xi[j] = u;
972: }
973: get_eg(&tmp0);
974: for ( i = 0; i < rank; i++ )
975: for ( j = 0; j < ri; j++ ) {
976: inner_product_mat_int_mod(a,wc,rank,i,j,&u);
977: addq(b[i][j],u,&s);
978: if ( s ) {
979: t = divin(NM(s),md,&tn);
980: if ( t )
981: error("generic_gauss_elim_hensel:incosistent");
982: NTOQ(tn,SGN(s),b[i][j]);
983: } else
984: b[i][j] = 0;
985: }
986: get_eg(&tmp1);
987: add_eg(&eg_mul,&tmp0,&tmp1);
988: /* q = q*md */
989: mulq(q,mdq,&u); q = u;
990: if ( !(count % 2) && intmtoratm_q(xmat,NM(q),*nmmat,dn) ) {
991: for ( j = k = l = 0; j < col; j++ )
992: if ( cinfo[j] )
993: rind[k++] = j;
994: else
995: cind[l++] = j;
996: if ( gensolve_check(mat,*nmmat,*dn,rind,cind) ) {
997: fprintf(stderr,"\n");
998: print_eg("INV",&eg_inv);
999: print_eg("MUL",&eg_mul);
1000: fflush(asir_out);
1001: return rank;
1002: }
1.1 noro 1003: }
1004: }
1005: }
1006: }
1007:
1008: int f4_nocheck;
1009:
1010: int gensolve_check(mat,nm,dn,rind,cind)
1011: MAT mat,nm;
1012: Q dn;
1013: int *rind,*cind;
1014: {
1015: int row,col,rank,clen,i,j,k,l;
1016: Q s,t,u;
1017: Q *w;
1018: Q *mati,*nmk;
1019:
1020: if ( f4_nocheck )
1021: return 1;
1022: row = mat->row; col = mat->col;
1023: rank = nm->row; clen = nm->col;
1024: w = (Q *)MALLOC(clen*sizeof(Q));
1025: for ( i = 0; i < row; i++ ) {
1026: mati = (Q *)mat->body[i];
1027: #if 1
1028: bzero(w,clen*sizeof(Q));
1029: for ( k = 0; k < rank; k++ )
1030: for ( l = 0, nmk = (Q *)nm->body[k]; l < clen; l++ ) {
1031: mulq(mati[rind[k]],nmk[l],&t);
1032: addq(w[l],t,&s); w[l] = s;
1033: }
1034: for ( j = 0; j < clen; j++ ) {
1035: mulq(dn,mati[cind[j]],&t);
1036: if ( cmpq(w[j],t) )
1037: break;
1038: }
1039: #else
1040: for ( j = 0; j < clen; j++ ) {
1041: for ( k = 0, s = 0; k < rank; k++ ) {
1042: mulq(mati[rind[k]],nm->body[k][j],&t);
1043: addq(s,t,&u); s = u;
1044: }
1045: mulq(dn,mati[cind[j]],&t);
1046: if ( cmpq(s,t) )
1047: break;
1048: }
1049: #endif
1050: if ( j != clen )
1051: break;
1052: }
1053: if ( i != row )
1054: return 0;
1055: else
1056: return 1;
1057: }
1058:
1059: /* assuming 0 < c < m */
1060:
1061: int inttorat(c,m,b,sgnp,nmp,dnp)
1062: N c,m,b;
1063: int *sgnp;
1064: N *nmp,*dnp;
1065: {
1066: Q qq,t,u1,v1,r1,nm;
1067: N q,r,u2,v2,r2;
1068:
1069: u1 = 0; v1 = ONE; u2 = m; v2 = c;
1070: while ( cmpn(v2,b) >= 0 ) {
1071: divn(u2,v2,&q,&r2); u2 = v2; v2 = r2;
1072: NTOQ(q,1,qq); mulq(qq,v1,&t); subq(u1,t,&r1); u1 = v1; v1 = r1;
1073: }
1074: if ( cmpn(NM(v1),b) >= 0 )
1075: return 0;
1076: else {
1077: *nmp = v2;
1078: *dnp = NM(v1);
1079: *sgnp = SGN(v1);
1080: return 1;
1081: }
1082: }
1083:
1084: /* mat->body = N ** */
1085:
1086: int intmtoratm(mat,md,nm,dn)
1087: MAT mat;
1088: N md;
1089: MAT nm;
1090: Q *dn;
1091: {
1092: N t,s,b;
1093: Q bound,dn0,dn1,nm1,q,tq;
1094: int i,j,k,l,row,col;
1095: Q **rmat;
1096: N **tmat;
1097: N *tmi;
1098: Q *nmk;
1099: N u,unm,udn;
1100: int sgn,ret;
1101:
1.3 noro 1102: if ( UNIN(md) )
1103: return 0;
1.1 noro 1104: row = mat->row; col = mat->col;
1105: bshiftn(md,1,&t);
1106: isqrt(t,&s);
1107: bshiftn(s,64,&b);
1108: if ( !b )
1109: b = ONEN;
1110: dn0 = ONE;
1111: tmat = (N **)mat->body;
1112: rmat = (Q **)nm->body;
1113: for ( i = 0; i < row; i++ )
1114: for ( j = 0, tmi = tmat[i]; j < col; j++ )
1115: if ( tmi[j] ) {
1116: muln(tmi[j],NM(dn0),&s);
1117: remn(s,md,&u);
1118: ret = inttorat(u,md,b,&sgn,&unm,&udn);
1119: if ( !ret )
1120: return 0;
1121: else {
1122: NTOQ(unm,sgn,nm1);
1123: NTOQ(udn,1,dn1);
1124: if ( !UNIQ(dn1) ) {
1125: for ( k = 0; k < i; k++ )
1126: for ( l = 0, nmk = rmat[k]; l < col; l++ ) {
1127: mulq(nmk[l],dn1,&q); nmk[l] = q;
1128: }
1129: for ( l = 0, nmk = rmat[i]; l < j; l++ ) {
1130: mulq(nmk[l],dn1,&q); nmk[l] = q;
1131: }
1132: }
1133: rmat[i][j] = nm1;
1134: mulq(dn0,dn1,&q); dn0 = q;
1135: }
1136: }
1137: *dn = dn0;
1138: return 1;
1139: }
1140:
1.3 noro 1141: /* mat->body = Q ** */
1142:
1143: int intmtoratm_q(mat,md,nm,dn)
1144: MAT mat;
1145: N md;
1146: MAT nm;
1147: Q *dn;
1148: {
1149: N t,s,b;
1150: Q bound,dn0,dn1,nm1,q,tq;
1151: int i,j,k,l,row,col;
1152: Q **rmat;
1153: Q **tmat;
1154: Q *tmi;
1155: Q *nmk;
1156: N u,unm,udn;
1157: int sgn,ret;
1158:
1159: if ( UNIN(md) )
1160: return 0;
1161: row = mat->row; col = mat->col;
1162: bshiftn(md,1,&t);
1163: isqrt(t,&s);
1164: bshiftn(s,64,&b);
1165: if ( !b )
1166: b = ONEN;
1167: dn0 = ONE;
1168: tmat = (Q **)mat->body;
1169: rmat = (Q **)nm->body;
1170: for ( i = 0; i < row; i++ )
1171: for ( j = 0, tmi = tmat[i]; j < col; j++ )
1172: if ( tmi[j] ) {
1173: muln(NM(tmi[j]),NM(dn0),&s);
1174: remn(s,md,&u);
1175: ret = inttorat(u,md,b,&sgn,&unm,&udn);
1176: if ( !ret )
1177: return 0;
1178: else {
1179: if ( SGN(tmi[j])<0 )
1180: sgn = -sgn;
1181: NTOQ(unm,sgn,nm1);
1182: NTOQ(udn,1,dn1);
1183: if ( !UNIQ(dn1) ) {
1184: for ( k = 0; k < i; k++ )
1185: for ( l = 0, nmk = rmat[k]; l < col; l++ ) {
1186: mulq(nmk[l],dn1,&q); nmk[l] = q;
1187: }
1188: for ( l = 0, nmk = rmat[i]; l < j; l++ ) {
1189: mulq(nmk[l],dn1,&q); nmk[l] = q;
1190: }
1191: }
1192: rmat[i][j] = nm1;
1193: mulq(dn0,dn1,&q); dn0 = q;
1194: }
1195: }
1196: *dn = dn0;
1197: return 1;
1198: }
1199:
1.4 noro 1200: #define ONE_STEP1 if ( zzz = *s ) { DMAR(zzz,hc,*tj,md,*tj) } tj++; s++;
1201:
1202: void reduce_reducers_mod(mat,row,col,md)
1203: int **mat;
1204: int row,col;
1205: int md;
1206: {
1207: int i,j,k,l,hc,zzz;
1208: int *t,*s,*tj,*ind;
1209:
1210: /* reduce the reducers */
1211: ind = (int *)ALLOCA(row*sizeof(int));
1212: for ( i = 0; i < row; i++ ) {
1213: t = mat[i];
1214: for ( j = 0; j < col && !t[j]; j++ );
1215: /* register the position of the head term */
1216: ind[i] = j;
1217: for ( l = i-1; l >= 0; l-- ) {
1218: /* reduce mat[i] by mat[l] */
1219: if ( hc = t[ind[l]] ) {
1220: /* mat[i] = mat[i]-hc*mat[l] */
1221: j = ind[l];
1222: s = mat[l]+j;
1223: tj = t+j;
1224: hc = md-hc;
1225: k = col-j;
1226: for ( ; k >= 64; k -= 64 ) {
1227: ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
1228: ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
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: }
1244: for ( ; k >= 0; k-- ) {
1245: if ( zzz = *s ) { DMAR(zzz,hc,*tj,md,*tj) } tj++; s++;
1246: }
1247: }
1248: }
1249: }
1250: }
1251:
1252: /*
1253: mat[i] : reducers (i=0,...,nred-1)
1254: spolys (i=nred,...,row-1)
1255: mat[0] < mat[1] < ... < mat[nred-1] w.r.t the term order
1256: 1. reduce the reducers
1257: 2. reduce spolys by the reduced reducers
1258: */
1259:
1260: void pre_reduce_mod(mat,row,col,nred,md)
1261: int **mat;
1262: int row,col,nred;
1263: int md;
1264: {
1265: int i,j,k,l,hc,inv;
1266: int *t,*s,*tk,*ind;
1267:
1268: #if 1
1269: /* reduce the reducers */
1270: ind = (int *)ALLOCA(row*sizeof(int));
1271: for ( i = 0; i < nred; i++ ) {
1272: /* make mat[i] monic and mat[i] by mat[0],...,mat[i-1] */
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: inv = invm(t[j],md);
1278: for ( k = j; k < col; k++ )
1279: if ( t[k] )
1280: DMAR(t[k],inv,0,md,t[k])
1281: for ( l = i-1; l >= 0; l-- ) {
1282: /* reduce mat[i] by mat[l] */
1283: if ( hc = t[ind[l]] ) {
1284: /* mat[i] = mat[i]-hc*mat[l] */
1285: for ( k = ind[l], hc = md-hc, s = mat[l]+k, tk = t+k;
1286: k < col; k++, tk++, s++ )
1287: if ( *s )
1288: DMAR(*s,hc,*tk,md,*tk)
1289: }
1290: }
1291: }
1292: /* reduce the spolys */
1293: for ( i = nred; i < row; i++ ) {
1294: t = mat[i];
1295: for ( l = nred-1; l >= 0; l-- ) {
1296: /* reduce mat[i] by mat[l] */
1297: if ( hc = t[ind[l]] ) {
1298: /* mat[i] = mat[i]-hc*mat[l] */
1299: for ( k = ind[l], hc = md-hc, s = mat[l]+k, tk = t+k;
1300: k < col; k++, tk++, s++ )
1301: if ( *s )
1302: DMAR(*s,hc,*tk,md,*tk)
1303: }
1304: }
1305: }
1306: #endif
1307: }
1308: /*
1309: mat[i] : reducers (i=0,...,nred-1)
1310: mat[0] < mat[1] < ... < mat[nred-1] w.r.t the term order
1311: */
1312:
1313: void reduce_sp_by_red_mod(sp,redmat,ind,nred,col,md)
1314: int *sp,**redmat;
1315: int *ind;
1316: int nred,col;
1317: int md;
1318: {
1319: int i,j,k,hc,zzz;
1320: int *t,*s,*tj;
1321:
1322: /* reduce the spolys by redmat */
1323: for ( i = nred-1; i >= 0; i-- ) {
1324: /* reduce sp by redmat[i] */
1325: if ( hc = sp[ind[i]] ) {
1326: /* sp = sp-hc*redmat[i] */
1327: j = ind[i];
1328: hc = md-hc;
1329: s = redmat[i]+j;
1330: tj = sp+j;
1331: for ( k = col-j; k >= 0; k-- ) {
1332: if ( zzz = *s ) { DMAR(zzz,hc,*tj,md,*tj) } tj++; s++;
1333: }
1334: }
1335: }
1336: }
1337:
1338: #define ONE_STEP2 if ( zzz = *pk ) { DMAR(zzz,a,*tk,md,*tk) } pk++; tk++;
1339:
1.1 noro 1340: int generic_gauss_elim_mod(mat,row,col,md,colstat)
1341: int **mat;
1342: int row,col,md;
1343: int *colstat;
1344: {
1.4 noro 1345: int i,j,k,l,inv,a,rank,zzz;
1346: int *t,*pivot,*pk,*tk;
1.1 noro 1347:
1348: for ( rank = 0, j = 0; j < col; j++ ) {
1349: for ( i = rank; i < row && !mat[i][j]; i++ );
1350: if ( i == row ) {
1351: colstat[j] = 0;
1352: continue;
1353: } else
1354: colstat[j] = 1;
1355: if ( i != rank ) {
1356: t = mat[i]; mat[i] = mat[rank]; mat[rank] = t;
1357: }
1358: pivot = mat[rank];
1359: inv = invm(pivot[j],md);
1.4 noro 1360: for ( k = j, pk = pivot+k; k < col; k++, pk++ )
1361: if ( *pk ) {
1362: DMAR(*pk,inv,0,md,*pk)
1.1 noro 1363: }
1364: for ( i = rank+1; i < row; i++ ) {
1365: t = mat[i];
1.4 noro 1366: if ( a = t[j] ) {
1367: a = md - a; pk = pivot+j; tk = t+j;
1368: k = col-j;
1369: for ( ; k >= 64; k -= 64 ) {
1370: ONE_STEP2 ONE_STEP2 ONE_STEP2 ONE_STEP2
1371: ONE_STEP2 ONE_STEP2 ONE_STEP2 ONE_STEP2
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: }
1387: for ( ; k >= 0; k -- ) {
1388: if ( zzz = *pk ) { DMAR(zzz,a,*tk,md,*tk) } pk++; tk++;
1389: }
1390: }
1.1 noro 1391: }
1392: rank++;
1393: }
1394: for ( j = col-1, l = rank-1; j >= 0; j-- )
1395: if ( colstat[j] ) {
1396: pivot = mat[l];
1397: for ( i = 0; i < l; i++ ) {
1398: t = mat[i];
1.4 noro 1399: if ( a = t[j] ) {
1400: a = md-a; pk = pivot+j; tk = t+j;
1401: k = col-j;
1402: for ( ; k >= 64; k -= 64 ) {
1403: ONE_STEP2 ONE_STEP2 ONE_STEP2 ONE_STEP2
1404: ONE_STEP2 ONE_STEP2 ONE_STEP2 ONE_STEP2
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: }
1420: for ( ; k >= 0; k -- ) {
1421: if ( zzz = *pk ) { DMAR(zzz,a,*tk,md,*tk) } pk++; tk++;
1422: }
1423: }
1.1 noro 1424: }
1425: l--;
1426: }
1427: return rank;
1428: }
1429:
1430: /* LU decomposition; a[i][i] = 1/U[i][i] */
1431:
1432: int lu_gfmmat(mat,md,perm)
1433: GFMMAT mat;
1434: unsigned int md;
1435: int *perm;
1436: {
1437: int row,col;
1438: int i,j,k,l;
1439: unsigned int *t,*pivot;
1440: unsigned int **a;
1441: unsigned int inv,m;
1442:
1443: row = mat->row; col = mat->col;
1444: a = mat->body;
1445: bzero(perm,row*sizeof(int));
1446:
1447: for ( i = 0; i < row; i++ )
1448: perm[i] = i;
1449: for ( k = 0; k < col; k++ ) {
1450: for ( i = k; i < row && !a[i][k]; i++ );
1451: if ( i == row )
1452: return 0;
1453: if ( i != k ) {
1454: j = perm[i]; perm[i] = perm[k]; perm[k] = j;
1455: t = a[i]; a[i] = a[k]; a[k] = t;
1456: }
1457: pivot = a[k];
1458: pivot[k] = inv = invm(pivot[k],md);
1459: for ( i = k+1; i < row; i++ ) {
1460: t = a[i];
1461: if ( m = t[k] ) {
1462: DMAR(inv,m,0,md,t[k])
1463: for ( j = k+1, m = md - t[k]; j < col; j++ )
1464: if ( pivot[j] ) {
1465: DMAR(m,pivot[j],t[j],md,t[j])
1466: }
1467: }
1468: }
1469: }
1470: return 1;
1471: }
1472:
1.3 noro 1473: /*
1474: Input
1475: a: a row x col matrix
1476: md : a modulus
1477:
1478: Output:
1479: return : d = the rank of mat
1480: a[0..(d-1)][0..(d-1)] : LU decomposition (a[i][i] = 1/U[i][i])
1481: rinfo: array of length row
1482: cinfo: array of length col
1483: i-th row in new a <-> rinfo[i]-th row in old a
1484: cinfo[j]=1 <=> j-th column is contained in the LU decomp.
1485: */
1486:
1487: int find_lhs_and_lu_mod(a,row,col,md,rinfo,cinfo)
1488: unsigned int **a;
1489: unsigned int md;
1490: int **rinfo,**cinfo;
1491: {
1492: int i,j,k,l,d;
1493: int *rp,*cp;
1494: unsigned int *t,*pivot;
1495: unsigned int inv,m;
1496:
1497: *rinfo = rp = (int *)MALLOC_ATOMIC(row*sizeof(int));
1498: *cinfo = cp = (int *)MALLOC_ATOMIC(col*sizeof(int));
1499: for ( i = 0; i < row; i++ )
1500: rp[i] = i;
1501: for ( k = 0, d = 0; k < col; k++ ) {
1502: for ( i = d; i < row && !a[i][k]; i++ );
1503: if ( i == row ) {
1504: cp[k] = 0;
1505: continue;
1506: } else
1507: cp[k] = 1;
1508: if ( i != d ) {
1509: j = rp[i]; rp[i] = rp[d]; rp[d] = j;
1510: t = a[i]; a[i] = a[d]; a[d] = t;
1511: }
1512: pivot = a[d];
1513: pivot[k] = inv = invm(pivot[k],md);
1514: for ( i = d+1; i < row; i++ ) {
1515: t = a[i];
1516: if ( m = t[k] ) {
1517: DMAR(inv,m,0,md,t[k])
1518: for ( j = k+1, m = md - t[k]; j < col; j++ )
1519: if ( pivot[j] ) {
1520: DMAR(m,pivot[j],t[j],md,t[j])
1521: }
1522: }
1523: }
1524: d++;
1525: }
1526: return d;
1527: }
1528:
1529: /*
1530: Input
1531: a : n x n matrix; a result of LU-decomposition
1532: md : modulus
1533: b : n x l matrix
1534: Output
1535: b = a^(-1)b
1536: */
1537:
1538: void solve_by_lu_mod(a,n,md,b,l)
1539: int **a;
1540: int n;
1541: int md;
1542: int **b;
1543: int l;
1544: {
1545: unsigned int *y,*c;
1546: int i,j,k;
1547: unsigned int t,m,m2;
1548:
1549: y = (int *)MALLOC_ATOMIC(n*sizeof(int));
1550: c = (int *)MALLOC_ATOMIC(n*sizeof(int));
1551: m2 = md>>1;
1552: for ( k = 0; k < l; k++ ) {
1553: /* copy b[.][k] to c */
1554: for ( i = 0; i < n; i++ )
1555: c[i] = (unsigned int)b[i][k];
1556: /* solve Ly=c */
1557: for ( i = 0; i < n; i++ ) {
1558: for ( t = c[i], j = 0; j < i; j++ )
1559: if ( a[i][j] ) {
1560: m = md - a[i][j];
1561: DMAR(m,y[j],t,md,t)
1562: }
1563: y[i] = t;
1564: }
1565: /* solve Uc=y */
1566: for ( i = n-1; i >= 0; i-- ) {
1567: for ( t = y[i], j =i+1; j < n; j++ )
1568: if ( a[i][j] ) {
1569: m = md - a[i][j];
1570: DMAR(m,c[j],t,md,t)
1571: }
1572: /* a[i][i] = 1/U[i][i] */
1573: DMAR(t,a[i][i],0,md,c[i])
1574: }
1575: /* copy c to b[.][k] with normalization */
1576: for ( i = 0; i < n; i++ )
1577: b[i][k] = (int)(c[i]>m2 ? c[i]-md : c[i]);
1578: }
1579: }
1580:
1.1 noro 1581: void Pleqm1(arg,rp)
1582: NODE arg;
1583: VECT *rp;
1584: {
1585: MAT m;
1586: VECT vect;
1587: pointer **mat;
1588: Q *v;
1589: Q q;
1590: int **wmat;
1591: int md,i,j,row,col,t,n,status;
1592:
1593: asir_assert(ARG0(arg),O_MAT,"leqm1");
1594: asir_assert(ARG1(arg),O_N,"leqm1");
1595: m = (MAT)ARG0(arg); md = QTOS((Q)ARG1(arg));
1596: row = m->row; col = m->col; mat = m->body;
1597: wmat = (int **)almat(row,col);
1598: for ( i = 0; i < row; i++ )
1599: for ( j = 0; j < col; j++ )
1600: if ( q = (Q)mat[i][j] ) {
1601: t = rem(NM(q),md);
1602: if ( SGN(q) < 0 )
1603: t = (md - t) % md;
1604: wmat[i][j] = t;
1605: } else
1606: wmat[i][j] = 0;
1607: status = gauss_elim_mod1(wmat,row,col,md);
1608: if ( status < 0 )
1609: *rp = 0;
1610: else if ( status > 0 )
1611: *rp = (VECT)ONE;
1612: else {
1613: n = col - 1;
1614: MKVECT(vect,n);
1615: for ( i = 0, v = (Q *)vect->body; i < n; i++ ) {
1616: t = (md-wmat[i][n])%md; STOQ(t,v[i]);
1617: }
1618: *rp = vect;
1619: }
1620: }
1621:
1622: gauss_elim_mod1(mat,row,col,md)
1623: int **mat;
1624: int row,col,md;
1625: {
1626: int i,j,k,inv,a,n;
1627: int *t,*pivot;
1628:
1629: n = col - 1;
1630: for ( j = 0; j < n; j++ ) {
1631: for ( i = j; i < row && !mat[i][j]; i++ );
1632: if ( i == row )
1633: return 1;
1634: if ( i != j ) {
1635: t = mat[i]; mat[i] = mat[j]; mat[j] = t;
1636: }
1637: pivot = mat[j];
1638: inv = invm(pivot[j],md);
1639: for ( k = j; k <= n; k++ )
1640: pivot[k] = dmar(pivot[k],inv,0,md);
1641: for ( i = j+1; i < row; i++ ) {
1642: t = mat[i];
1643: if ( i != j && (a = t[j]) )
1644: for ( k = j, a = md - a; k <= n; k++ )
1645: t[k] = dmar(pivot[k],a,t[k],md);
1646: }
1647: }
1648: for ( i = n; i < row && !mat[i][n]; i++ );
1649: if ( i == row ) {
1650: for ( j = n-1; j >= 0; j-- ) {
1651: for ( i = j-1, a = (md-mat[j][n])%md; i >= 0; i-- ) {
1652: mat[i][n] = dmar(mat[i][j],a,mat[i][n],md);
1653: mat[i][j] = 0;
1654: }
1655: }
1656: return 0;
1657: } else
1658: return -1;
1659: }
1660:
1661: void Pgeninvm(arg,rp)
1662: NODE arg;
1663: LIST *rp;
1664: {
1665: MAT m;
1666: pointer **mat;
1667: Q **tmat;
1668: Q q;
1669: unsigned int **wmat;
1670: int md,i,j,row,col,t,status;
1671: MAT mat1,mat2;
1672: NODE node1,node2;
1673:
1674: asir_assert(ARG0(arg),O_MAT,"leqm1");
1675: asir_assert(ARG1(arg),O_N,"leqm1");
1676: m = (MAT)ARG0(arg); md = QTOS((Q)ARG1(arg));
1677: row = m->row; col = m->col; mat = m->body;
1678: wmat = (unsigned int **)almat(row,col+row);
1679: for ( i = 0; i < row; i++ ) {
1680: bzero((char *)wmat[i],(col+row)*sizeof(int));
1681: for ( j = 0; j < col; j++ )
1682: if ( q = (Q)mat[i][j] ) {
1683: t = rem(NM(q),md);
1684: if ( SGN(q) < 0 )
1685: t = (md - t) % md;
1686: wmat[i][j] = t;
1687: }
1688: wmat[i][col+i] = 1;
1689: }
1690: status = gauss_elim_geninv_mod(wmat,row,col,md);
1691: if ( status > 0 )
1692: *rp = 0;
1693: else {
1694: MKMAT(mat1,col,row); MKMAT(mat2,row-col,row);
1695: for ( i = 0, tmat = (Q **)mat1->body; i < col; i++ )
1696: for ( j = 0; j < row; j++ )
1697: STOQ(wmat[i][j+col],tmat[i][j]);
1698: for ( tmat = (Q **)mat2->body; i < row; i++ )
1699: for ( j = 0; j < row; j++ )
1700: STOQ(wmat[i][j+col],tmat[i-col][j]);
1701: MKNODE(node2,mat2,0); MKNODE(node1,mat1,node2); MKLIST(*rp,node1);
1702: }
1703: }
1704:
1705: int gauss_elim_geninv_mod(mat,row,col,md)
1706: unsigned int **mat;
1707: int row,col,md;
1708: {
1709: int i,j,k,inv,a,n,m;
1710: unsigned int *t,*pivot;
1711:
1712: n = col; m = row+col;
1713: for ( j = 0; j < n; j++ ) {
1714: for ( i = j; i < row && !mat[i][j]; i++ );
1715: if ( i == row )
1716: return 1;
1717: if ( i != j ) {
1718: t = mat[i]; mat[i] = mat[j]; mat[j] = t;
1719: }
1720: pivot = mat[j];
1721: inv = invm(pivot[j],md);
1722: for ( k = j; k < m; k++ )
1723: pivot[k] = dmar(pivot[k],inv,0,md);
1724: for ( i = j+1; i < row; i++ ) {
1725: t = mat[i];
1726: if ( a = t[j] )
1727: for ( k = j, a = md - a; k < m; k++ )
1728: t[k] = dmar(pivot[k],a,t[k],md);
1729: }
1730: }
1731: for ( j = n-1; j >= 0; j-- ) {
1732: pivot = mat[j];
1733: for ( i = j-1; i >= 0; i-- ) {
1734: t = mat[i];
1735: if ( a = t[j] )
1736: for ( k = j, a = md - a; k < m; k++ )
1737: t[k] = dmar(pivot[k],a,t[k],md);
1738: }
1739: }
1740: return 0;
1741: }
1742:
1743: void Psolve_by_lu_gfmmat(arg,rp)
1744: NODE arg;
1745: VECT *rp;
1746: {
1747: GFMMAT lu;
1748: Q *perm,*rhs,*v;
1749: int n,i;
1750: unsigned int md;
1751: unsigned int *b,*sol;
1752: VECT r;
1753:
1754: lu = (GFMMAT)ARG0(arg);
1755: perm = (Q *)BDY((VECT)ARG1(arg));
1756: rhs = (Q *)BDY((VECT)ARG2(arg));
1757: md = (unsigned int)QTOS((Q)ARG3(arg));
1758: n = lu->col;
1759: b = (unsigned int *)MALLOC_ATOMIC(n*sizeof(int));
1760: sol = (unsigned int *)MALLOC_ATOMIC(n*sizeof(int));
1761: for ( i = 0; i < n; i++ )
1762: b[i] = QTOS(rhs[QTOS(perm[i])]);
1763: solve_by_lu_gfmmat(lu,md,b,sol);
1764: MKVECT(r,n);
1765: for ( i = 0, v = (Q *)r->body; i < n; i++ )
1766: STOQ(sol[i],v[i]);
1767: *rp = r;
1768: }
1769:
1770: void solve_by_lu_gfmmat(lu,md,b,x)
1771: GFMMAT lu;
1772: unsigned int md;
1773: unsigned int *b;
1774: unsigned int *x;
1775: {
1776: int n;
1777: unsigned int **a;
1778: unsigned int *y;
1779: int i,j;
1780: unsigned int t,m;
1781:
1782: n = lu->col;
1783: a = lu->body;
1784: y = (unsigned int *)MALLOC_ATOMIC(n*sizeof(int));
1785: /* solve Ly=b */
1786: for ( i = 0; i < n; i++ ) {
1787: for ( t = b[i], j = 0; j < i; j++ )
1788: if ( a[i][j] ) {
1789: m = md - a[i][j];
1790: DMAR(m,y[j],t,md,t)
1791: }
1792: y[i] = t;
1793: }
1794: /* solve Ux=y */
1795: for ( i = n-1; i >= 0; i-- ) {
1796: for ( t = y[i], j =i+1; j < n; j++ )
1797: if ( a[i][j] ) {
1798: m = md - a[i][j];
1799: DMAR(m,x[j],t,md,t)
1800: }
1801: /* a[i][i] = 1/U[i][i] */
1802: DMAR(t,a[i][i],0,md,x[i])
1803: }
1804: }
1805:
1806: void Plu_gfmmat(arg,rp)
1807: NODE arg;
1808: LIST *rp;
1809: {
1810: MAT m;
1811: GFMMAT mm;
1812: unsigned int md;
1813: int i,row,col,status;
1814: int *iperm;
1815: Q *v;
1816: VECT perm;
1817: NODE n0;
1818:
1819: asir_assert(ARG0(arg),O_MAT,"mat_to_gfmmat");
1820: asir_assert(ARG1(arg),O_N,"mat_to_gfmmat");
1821: m = (MAT)ARG0(arg); md = (unsigned int)QTOS((Q)ARG1(arg));
1822: mat_to_gfmmat(m,md,&mm);
1823: row = m->row;
1824: col = m->col;
1825: iperm = (int *)MALLOC_ATOMIC(row*sizeof(int));
1826: status = lu_gfmmat(mm,md,iperm);
1827: if ( !status )
1828: n0 = 0;
1829: else {
1830: MKVECT(perm,row);
1831: for ( i = 0, v = (Q *)perm->body; i < row; i++ )
1832: STOQ(iperm[i],v[i]);
1833: n0 = mknode(2,mm,perm);
1834: }
1835: MKLIST(*rp,n0);
1836: }
1837:
1838: void Pmat_to_gfmmat(arg,rp)
1839: NODE arg;
1840: GFMMAT *rp;
1841: {
1842: MAT m;
1843: unsigned int md;
1844:
1845: asir_assert(ARG0(arg),O_MAT,"mat_to_gfmmat");
1846: asir_assert(ARG1(arg),O_N,"mat_to_gfmmat");
1847: m = (MAT)ARG0(arg); md = (unsigned int)QTOS((Q)ARG1(arg));
1848: mat_to_gfmmat(m,md,rp);
1849: }
1850:
1851: void mat_to_gfmmat(m,md,rp)
1852: MAT m;
1853: unsigned int md;
1854: GFMMAT *rp;
1855: {
1856: unsigned int **wmat;
1857: unsigned int t;
1858: Q **mat;
1859: Q q;
1860: int i,j,row,col;
1861:
1862: row = m->row; col = m->col; mat = (Q **)m->body;
1863: wmat = (unsigned int **)almat(row,col);
1864: for ( i = 0; i < row; i++ ) {
1865: bzero((char *)wmat[i],col*sizeof(unsigned int));
1866: for ( j = 0; j < col; j++ )
1867: if ( q = mat[i][j] ) {
1868: t = (unsigned int)rem(NM(q),md);
1869: if ( SGN(q) < 0 )
1870: t = (md - t) % md;
1871: wmat[i][j] = t;
1872: }
1873: }
1874: TOGFMMAT(row,col,wmat,*rp);
1875: }
1876:
1877: void Pgeninvm_swap(arg,rp)
1878: NODE arg;
1879: LIST *rp;
1880: {
1881: MAT m;
1882: pointer **mat;
1883: Q **tmat;
1884: Q *tvect;
1885: Q q;
1886: unsigned int **wmat,**invmat;
1887: int *index;
1888: unsigned int t,md;
1889: int i,j,row,col,status;
1890: MAT mat1;
1891: VECT vect1;
1892: NODE node1,node2;
1893:
1894: asir_assert(ARG0(arg),O_MAT,"geninvm_swap");
1895: asir_assert(ARG1(arg),O_N,"geninvm_swap");
1896: m = (MAT)ARG0(arg); md = QTOS((Q)ARG1(arg));
1897: row = m->row; col = m->col; mat = m->body;
1898: wmat = (unsigned int **)almat(row,col+row);
1899: for ( i = 0; i < row; i++ ) {
1900: bzero((char *)wmat[i],(col+row)*sizeof(int));
1901: for ( j = 0; j < col; j++ )
1902: if ( q = (Q)mat[i][j] ) {
1903: t = (unsigned int)rem(NM(q),md);
1904: if ( SGN(q) < 0 )
1905: t = (md - t) % md;
1906: wmat[i][j] = t;
1907: }
1908: wmat[i][col+i] = 1;
1909: }
1910: status = gauss_elim_geninv_mod_swap(wmat,row,col,md,&invmat,&index);
1911: if ( status > 0 )
1912: *rp = 0;
1913: else {
1914: MKMAT(mat1,col,col);
1915: for ( i = 0, tmat = (Q **)mat1->body; i < col; i++ )
1916: for ( j = 0; j < col; j++ )
1917: UTOQ(invmat[i][j],tmat[i][j]);
1918: MKVECT(vect1,row);
1919: for ( i = 0, tvect = (Q *)vect1->body; i < row; i++ )
1920: STOQ(index[i],tvect[i]);
1921: MKNODE(node2,vect1,0); MKNODE(node1,mat1,node2); MKLIST(*rp,node1);
1922: }
1923: }
1924:
1925: gauss_elim_geninv_mod_swap(mat,row,col,md,invmatp,indexp)
1926: unsigned int **mat;
1927: int row,col;
1928: unsigned int md;
1929: unsigned int ***invmatp;
1930: int **indexp;
1931: {
1932: int i,j,k,inv,a,n,m;
1933: unsigned int *t,*pivot,*s;
1934: int *index;
1935: unsigned int **invmat;
1936:
1937: n = col; m = row+col;
1938: *indexp = index = (int *)MALLOC_ATOMIC(row*sizeof(int));
1939: for ( i = 0; i < row; i++ )
1940: index[i] = i;
1941: for ( j = 0; j < n; j++ ) {
1942: for ( i = j; i < row && !mat[i][j]; i++ );
1943: if ( i == row ) {
1944: *indexp = 0; *invmatp = 0; return 1;
1945: }
1946: if ( i != j ) {
1947: t = mat[i]; mat[i] = mat[j]; mat[j] = t;
1948: k = index[i]; index[i] = index[j]; index[j] = k;
1949: }
1950: pivot = mat[j];
1951: inv = (unsigned int)invm(pivot[j],md);
1952: for ( k = j; k < m; k++ )
1953: if ( pivot[k] )
1954: pivot[k] = (unsigned int)dmar(pivot[k],inv,0,md);
1955: for ( i = j+1; i < row; i++ ) {
1956: t = mat[i];
1957: if ( a = t[j] )
1958: for ( k = j, a = md - a; k < m; k++ )
1959: if ( pivot[k] )
1960: t[k] = dmar(pivot[k],a,t[k],md);
1961: }
1962: }
1963: for ( j = n-1; j >= 0; j-- ) {
1964: pivot = mat[j];
1965: for ( i = j-1; i >= 0; i-- ) {
1966: t = mat[i];
1967: if ( a = t[j] )
1968: for ( k = j, a = md - a; k < m; k++ )
1969: if ( pivot[k] )
1970: t[k] = dmar(pivot[k],a,t[k],md);
1971: }
1972: }
1973: *invmatp = invmat = (unsigned int **)almat(col,col);
1974: for ( i = 0; i < col; i++ )
1975: for ( j = 0, s = invmat[i], t = mat[i]; j < col; j++ )
1976: s[j] = t[col+index[j]];
1977: return 0;
1978: }
1979:
1980: void _addn(N,N,N);
1981: int _subn(N,N,N);
1982: void _muln(N,N,N);
1983:
1984: void inner_product_int(a,b,n,r)
1985: Q *a,*b;
1986: int n;
1987: Q *r;
1988: {
1989: int la,lb,i;
1990: int sgn,sgn1;
1991: N wm,wma,sum,t;
1992:
1993: for ( la = lb = 0, i = 0; i < n; i++ ) {
1994: if ( a[i] )
1995: if ( DN(a[i]) )
1996: error("inner_product_int : invalid argument");
1997: else
1998: la = MAX(PL(NM(a[i])),la);
1999: if ( b[i] )
2000: if ( DN(b[i]) )
2001: error("inner_product_int : invalid argument");
2002: else
2003: lb = MAX(PL(NM(b[i])),lb);
2004: }
2005: sgn = 0;
2006: sum= NALLOC(la+lb+2);
2007: bzero((char *)sum,(la+lb+3)*sizeof(unsigned int));
2008: wm = NALLOC(la+lb+2);
2009: wma = NALLOC(la+lb+2);
2010: for ( i = 0; i < n; i++ ) {
2011: if ( !a[i] || !b[i] )
2012: continue;
2013: _muln(NM(a[i]),NM(b[i]),wm);
2014: sgn1 = SGN(a[i])*SGN(b[i]);
2015: if ( !sgn ) {
2016: sgn = sgn1;
2017: t = wm; wm = sum; sum = t;
2018: } else if ( sgn == sgn1 ) {
2019: _addn(sum,wm,wma);
2020: if ( !PL(wma) )
2021: sgn = 0;
2022: t = wma; wma = sum; sum = t;
2023: } else {
2024: /* sgn*sum+sgn1*wm = sgn*(sum-wm) */
2025: sgn *= _subn(sum,wm,wma);
2026: t = wma; wma = sum; sum = t;
2027: }
2028: }
2029: GC_free(wm);
2030: GC_free(wma);
2031: if ( !sgn ) {
2032: GC_free(sum);
2033: *r = 0;
2034: } else
2035: NTOQ(sum,sgn,*r);
2036: }
2037:
1.3 noro 2038: /* (k,l) element of a*b where a: .x n matrix, b: n x . integer matrix */
2039:
2040: void inner_product_mat_int_mod(a,b,n,k,l,r)
2041: Q **a;
2042: int **b;
2043: int n,k,l;
2044: Q *r;
2045: {
2046: int la,lb,i;
2047: int sgn,sgn1;
2048: N wm,wma,sum,t;
2049: Q aki;
2050: int bil,bilsgn;
2051: struct oN tn;
2052:
2053: for ( la = 0, i = 0; i < n; i++ ) {
2054: if ( aki = a[k][i] )
2055: if ( DN(aki) )
2056: error("inner_product_int : invalid argument");
2057: else
2058: la = MAX(PL(NM(aki)),la);
2059: }
2060: lb = 1;
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 ( !(aki = a[k][i]) || !(bil = b[i][l]) )
2068: continue;
2069: tn.p = 1;
2070: if ( bil > 0 ) {
2071: tn.b[0] = bil; bilsgn = 1;
2072: } else {
2073: tn.b[0] = -bil; bilsgn = -1;
2074: }
2075: _muln(NM(aki),&tn,wm);
2076: sgn1 = SGN(aki)*bilsgn;
2077: if ( !sgn ) {
2078: sgn = sgn1;
2079: t = wm; wm = sum; sum = t;
2080: } else if ( sgn == sgn1 ) {
2081: _addn(sum,wm,wma);
2082: if ( !PL(wma) )
2083: sgn = 0;
2084: t = wma; wma = sum; sum = t;
2085: } else {
2086: /* sgn*sum+sgn1*wm = sgn*(sum-wm) */
2087: sgn *= _subn(sum,wm,wma);
2088: t = wma; wma = sum; sum = t;
2089: }
2090: }
2091: GC_free(wm);
2092: GC_free(wma);
2093: if ( !sgn ) {
2094: GC_free(sum);
2095: *r = 0;
2096: } else
2097: NTOQ(sum,sgn,*r);
2098: }
2099:
1.1 noro 2100: void Pmul_mat_vect_int(arg,rp)
2101: NODE arg;
2102: VECT *rp;
2103: {
2104: MAT mat;
2105: VECT vect,r;
2106: int row,col,i;
2107:
2108: mat = (MAT)ARG0(arg);
2109: vect = (VECT)ARG1(arg);
2110: row = mat->row;
2111: col = mat->col;
2112: MKVECT(r,row);
2113: for ( i = 0; i < row; i++ )
2114: inner_product_int(mat->body[i],vect->body,col,&r->body[i]);
2115: *rp = r;
2116: }
2117:
2118: void Pnbpoly_up2(arg,rp)
2119: NODE arg;
2120: GF2N *rp;
2121: {
2122: int m,type,ret;
2123: UP2 r;
2124:
2125: m = QTOS((Q)ARG0(arg));
2126: type = QTOS((Q)ARG1(arg));
2127: ret = generate_ONB_polynomial(&r,m,type);
2128: if ( ret == 0 )
2129: MKGF2N(r,*rp);
2130: else
2131: *rp = 0;
2132: }
2133:
2134: void Px962_irredpoly_up2(arg,rp)
2135: NODE arg;
2136: GF2N *rp;
2137: {
2138: int m,type,ret,w;
2139: GF2N prev;
2140: UP2 r;
2141:
2142: m = QTOS((Q)ARG0(arg));
2143: prev = (GF2N)ARG1(arg);
2144: if ( !prev ) {
2145: w = (m>>5)+1; NEWUP2(r,w); r->w = 0;
2146: bzero((char *)r->b,w*sizeof(unsigned int));
2147: } else {
2148: r = prev->body;
2149: if ( degup2(r) != m ) {
2150: w = (m>>5)+1; NEWUP2(r,w); r->w = 0;
2151: bzero((char *)r->b,w*sizeof(unsigned int));
2152: }
2153: }
2154: ret = _generate_irreducible_polynomial(r,m,type);
2155: if ( ret == 0 )
2156: MKGF2N(r,*rp);
2157: else
2158: *rp = 0;
2159: }
2160:
2161: void Pirredpoly_up2(arg,rp)
2162: NODE arg;
2163: GF2N *rp;
2164: {
2165: int m,type,ret,w;
2166: GF2N prev;
2167: UP2 r;
2168:
2169: m = QTOS((Q)ARG0(arg));
2170: prev = (GF2N)ARG1(arg);
2171: if ( !prev ) {
2172: w = (m>>5)+1; NEWUP2(r,w); r->w = 0;
2173: bzero((char *)r->b,w*sizeof(unsigned int));
2174: } else {
2175: r = prev->body;
2176: if ( degup2(r) != m ) {
2177: w = (m>>5)+1; NEWUP2(r,w); r->w = 0;
2178: bzero((char *)r->b,w*sizeof(unsigned int));
2179: }
2180: }
2181: ret = _generate_good_irreducible_polynomial(r,m,type);
2182: if ( ret == 0 )
2183: MKGF2N(r,*rp);
2184: else
2185: *rp = 0;
2186: }
2187:
2188: /*
2189: * f = type 'type' normal polynomial of degree m if exists
2190: * IEEE P1363 A.7.2
2191: *
2192: * return value : 0 --- exists
2193: * 1 --- does not exist
2194: * -1 --- failure (memory allocation error)
2195: */
2196:
2197: int generate_ONB_polynomial(UP2 *rp,int m,int type)
2198: {
2199: int i,r;
2200: int w;
2201: UP2 f,f0,f1,f2,t;
2202:
2203: w = (m>>5)+1;
2204: switch ( type ) {
2205: case 1:
2206: if ( !TypeT_NB_check(m,1) ) return 1;
2207: NEWUP2(f,w); *rp = f; f->w = w;
2208: /* set all the bits */
2209: for ( i = 0; i < w; i++ )
2210: f->b[i] = 0xffffffff;
2211: /* mask the top word if necessary */
2212: if ( r = (m+1)&31 )
2213: f->b[w-1] &= (1<<r)-1;
2214: return 0;
2215: break;
2216: case 2:
2217: if ( !TypeT_NB_check(m,2) ) return 1;
2218: NEWUP2(f,w); *rp = f;
2219: W_NEWUP2(f0,w);
2220: W_NEWUP2(f1,w);
2221: W_NEWUP2(f2,w);
2222:
2223: /* recursion for genrating Type II normal polynomial */
2224:
2225: /* f0 = 1, f1 = t+1 */
2226: f0->w = 1; f0->b[0] = 1;
2227: f1->w = 1; f1->b[0] = 3;
2228: for ( i = 2; i <= m; i++ ) {
2229: /* f2 = t*f1+f0 */
2230: _bshiftup2(f1,-1,f2);
2231: _addup2_destructive(f2,f0);
2232: /* cyclic change of the variables */
2233: t = f0; f0 = f1; f1 = f2; f2 = t;
2234: }
2235: _copyup2(f1,f);
2236: return 0;
2237: break;
2238: default:
2239: return -1;
2240: break;
2241: }
2242: }
2243:
2244: /*
2245: * f = an irreducible trinomial or pentanomial of degree d 'after' f
2246: * return value : 0 --- exists
2247: * 1 --- does not exist (exhaustion)
2248: */
2249:
2250: int _generate_irreducible_polynomial(UP2 f,int d)
2251: {
2252: int ret,i,j,k,nz,i0,j0,k0;
2253: int w;
2254: unsigned int *fd;
2255:
2256: /*
2257: * if f = x^d+x^i+1 then i0 <- i, j0 <- 0, k0 <-0.
2258: * if f = x^d+x^k+x^j+x^i+1 (k>j>i) then i0 <- i, j0 <- j, k0 <-k.
2259: * otherwise i0,j0,k0 is set to 0.
2260: */
2261:
2262: fd = f->b;
2263: w = (d>>5)+1;
2264: if ( f->w && (d==degup2(f)) ) {
2265: for ( nz = 0, i = d; i >= 0; i-- )
2266: if ( fd[i>>5]&(1<<(i&31)) ) nz++;
2267: switch ( nz ) {
2268: case 3:
2269: for ( i0 = 1; !(fd[i0>>5]&(1<<(i0&31))) ; i0++ );
2270: /* reset i0-th bit */
2271: fd[i0>>5] &= ~(1<<(i0&31));
2272: j0 = k0 = 0;
2273: break;
2274: case 5:
2275: for ( i0 = 1; !(fd[i0>>5]&(1<<(i0&31))) ; i0++ );
2276: /* reset i0-th bit */
2277: fd[i0>>5] &= ~(1<<(i0&31));
2278: for ( j0 = i0+1; !(fd[j0>>5]&(1<<(j0&31))) ; j0++ );
2279: /* reset j0-th bit */
2280: fd[j0>>5] &= ~(1<<(j0&31));
2281: for ( k0 = j0+1; !(fd[k0>>5]&(1<<(k0&31))) ; k0++ );
2282: /* reset k0-th bit */
2283: fd[k0>>5] &= ~(1<<(k0&31));
2284: break;
2285: default:
2286: f->w = 0; break;
2287: }
2288: } else
2289: f->w = 0;
2290:
2291: if ( !f->w ) {
2292: fd = f->b;
2293: f->w = w; fd[0] |= 1; fd[d>>5] |= (1<<(d&31));
2294: i0 = j0 = k0 = 0;
2295: }
2296: /* if j0 > 0 then f is already a pentanomial */
2297: if ( j0 > 0 ) goto PENTA;
2298:
2299: /* searching for an irreducible trinomial */
2300:
2301: for ( i = 1; 2*i <= d; i++ ) {
2302: /* skip the polynomials 'before' f */
2303: if ( i < i0 ) continue;
2304: if ( i == i0 ) { i0 = 0; continue; }
2305: /* set i-th bit */
2306: fd[i>>5] |= (1<<(i&31));
2307: ret = irredcheck_dddup2(f);
2308: if ( ret == 1 ) return 0;
2309: /* reset i-th bit */
2310: fd[i>>5] &= ~(1<<(i&31));
2311: }
2312:
2313: /* searching for an irreducible pentanomial */
2314: PENTA:
2315: for ( i = 1; i < d; i++ ) {
2316: /* skip the polynomials 'before' f */
2317: if ( i < i0 ) continue;
2318: if ( i == i0 ) i0 = 0;
2319: /* set i-th bit */
2320: fd[i>>5] |= (1<<(i&31));
2321: for ( j = i+1; j < d; j++ ) {
2322: /* skip the polynomials 'before' f */
2323: if ( j < j0 ) continue;
2324: if ( j == j0 ) j0 = 0;
2325: /* set j-th bit */
2326: fd[j>>5] |= (1<<(j&31));
2327: for ( k = j+1; k < d; k++ ) {
2328: /* skip the polynomials 'before' f */
2329: if ( k < k0 ) continue;
2330: else if ( k == k0 ) { k0 = 0; continue; }
2331: /* set k-th bit */
2332: fd[k>>5] |= (1<<(k&31));
2333: ret = irredcheck_dddup2(f);
2334: if ( ret == 1 ) return 0;
2335: /* reset k-th bit */
2336: fd[k>>5] &= ~(1<<(k&31));
2337: }
2338: /* reset j-th bit */
2339: fd[j>>5] &= ~(1<<(j&31));
2340: }
2341: /* reset i-th bit */
2342: fd[i>>5] &= ~(1<<(i&31));
2343: }
2344: /* exhausted */
2345: return 1;
2346: }
2347:
2348: /*
2349: * f = an irreducible trinomial or pentanomial of degree d 'after' f
2350: *
2351: * searching strategy:
2352: * trinomial x^d+x^i+1:
2353: * i is as small as possible.
2354: * trinomial x^d+x^i+x^j+x^k+1:
2355: * i is as small as possible.
2356: * For such i, j is as small as possible.
2357: * For such i and j, 'k' is as small as possible.
2358: *
2359: * return value : 0 --- exists
2360: * 1 --- does not exist (exhaustion)
2361: */
2362:
2363: int _generate_good_irreducible_polynomial(UP2 f,int d)
2364: {
2365: int ret,i,j,k,nz,i0,j0,k0;
2366: int w;
2367: unsigned int *fd;
2368:
2369: /*
2370: * if f = x^d+x^i+1 then i0 <- i, j0 <- 0, k0 <-0.
2371: * if f = x^d+x^k+x^j+x^i+1 (k>j>i) then i0 <- i, j0 <- j, k0 <-k.
2372: * otherwise i0,j0,k0 is set to 0.
2373: */
2374:
2375: fd = f->b;
2376: w = (d>>5)+1;
2377: if ( f->w && (d==degup2(f)) ) {
2378: for ( nz = 0, i = d; i >= 0; i-- )
2379: if ( fd[i>>5]&(1<<(i&31)) ) nz++;
2380: switch ( nz ) {
2381: case 3:
2382: for ( i0 = 1; !(fd[i0>>5]&(1<<(i0&31))) ; i0++ );
2383: /* reset i0-th bit */
2384: fd[i0>>5] &= ~(1<<(i0&31));
2385: j0 = k0 = 0;
2386: break;
2387: case 5:
2388: for ( i0 = 1; !(fd[i0>>5]&(1<<(i0&31))) ; i0++ );
2389: /* reset i0-th bit */
2390: fd[i0>>5] &= ~(1<<(i0&31));
2391: for ( j0 = i0+1; !(fd[j0>>5]&(1<<(j0&31))) ; j0++ );
2392: /* reset j0-th bit */
2393: fd[j0>>5] &= ~(1<<(j0&31));
2394: for ( k0 = j0+1; !(fd[k0>>5]&(1<<(k0&31))) ; k0++ );
2395: /* reset k0-th bit */
2396: fd[k0>>5] &= ~(1<<(k0&31));
2397: break;
2398: default:
2399: f->w = 0; break;
2400: }
2401: } else
2402: f->w = 0;
2403:
2404: if ( !f->w ) {
2405: fd = f->b;
2406: f->w = w; fd[0] |= 1; fd[d>>5] |= (1<<(d&31));
2407: i0 = j0 = k0 = 0;
2408: }
2409: /* if j0 > 0 then f is already a pentanomial */
2410: if ( j0 > 0 ) goto PENTA;
2411:
2412: /* searching for an irreducible trinomial */
2413:
2414: for ( i = 1; 2*i <= d; i++ ) {
2415: /* skip the polynomials 'before' f */
2416: if ( i < i0 ) continue;
2417: if ( i == i0 ) { i0 = 0; continue; }
2418: /* set i-th bit */
2419: fd[i>>5] |= (1<<(i&31));
2420: ret = irredcheck_dddup2(f);
2421: if ( ret == 1 ) return 0;
2422: /* reset i-th bit */
2423: fd[i>>5] &= ~(1<<(i&31));
2424: }
2425:
2426: /* searching for an irreducible pentanomial */
2427: PENTA:
2428: for ( i = 3; i < d; i++ ) {
2429: /* skip the polynomials 'before' f */
2430: if ( i < i0 ) continue;
2431: if ( i == i0 ) i0 = 0;
2432: /* set i-th bit */
2433: fd[i>>5] |= (1<<(i&31));
2434: for ( j = 2; j < i; j++ ) {
2435: /* skip the polynomials 'before' f */
2436: if ( j < j0 ) continue;
2437: if ( j == j0 ) j0 = 0;
2438: /* set j-th bit */
2439: fd[j>>5] |= (1<<(j&31));
2440: for ( k = 1; k < j; k++ ) {
2441: /* skip the polynomials 'before' f */
2442: if ( k < k0 ) continue;
2443: else if ( k == k0 ) { k0 = 0; continue; }
2444: /* set k-th bit */
2445: fd[k>>5] |= (1<<(k&31));
2446: ret = irredcheck_dddup2(f);
2447: if ( ret == 1 ) return 0;
2448: /* reset k-th bit */
2449: fd[k>>5] &= ~(1<<(k&31));
2450: }
2451: /* reset j-th bit */
2452: fd[j>>5] &= ~(1<<(j&31));
2453: }
2454: /* reset i-th bit */
2455: fd[i>>5] &= ~(1<<(i&31));
2456: }
2457: /* exhausted */
2458: return 1;
1.3 noro 2459: }
2460:
2461: printqmat(mat,row,col)
2462: Q **mat;
2463: int row,col;
2464: {
2465: int i,j;
2466:
2467: for ( i = 0; i < row; i++ ) {
2468: for ( j = 0; j < col; j++ ) {
2469: printnum(mat[i][j]); printf(" ");
2470: }
2471: printf("\n");
2472: }
2473: }
2474:
2475: printimat(mat,row,col)
2476: int **mat;
2477: int row,col;
2478: {
2479: int i,j;
2480:
2481: for ( i = 0; i < row; i++ ) {
2482: for ( j = 0; j < col; j++ ) {
2483: printf("%d ",mat[i][j]);
2484: }
2485: printf("\n");
2486: }
1.1 noro 2487: }