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