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