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