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