Annotation of OpenXM_contrib2/asir2000/builtin/array.c, Revision 1.48
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.48 ! noro 48: * $OpenXM: OpenXM_contrib2/asir2000/builtin/array.c,v 1.47 2005/11/27 00:07:05 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.46 saito 161: void Pqsort(NODE arg,LIST *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: }
1.46 saito 208: MKLIST(*rp,n);
1.35 ohara 209: }else {
1.46 saito 210: *rp = (LIST)vect;
1.35 ohara 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;
1.43 saito 761: STOQ(m,q); MKNODE(s,q,0); STOQ(n,q); MKNODE(t,q,s);
762: break;
763: case O_IMAT:
764: n = ((IMAT)ARG0(arg))->row; m = ((IMAT)ARG0(arg))->col;
1.1 noro 765: STOQ(m,q); MKNODE(s,q,0); STOQ(n,q); MKNODE(t,q,s);
766: break;
767: default:
768: error("size : invalid argument"); break;
769: }
770: }
771: MKLIST(*rp,t);
772: }
773:
1.24 noro 774: void Pdet(NODE arg,P *rp)
1.1 noro 775: {
776: MAT m;
777: int n,i,j,mod;
778: P d;
779: P **mat,**w;
780:
781: m = (MAT)ARG0(arg);
782: asir_assert(m,O_MAT,"det");
783: if ( m->row != m->col )
784: error("det : non-square matrix");
785: else if ( argc(arg) == 1 )
786: detp(CO,(P **)BDY(m),m->row,rp);
787: else {
788: n = m->row; mod = QTOS((Q)ARG1(arg)); mat = (P **)BDY(m);
789: w = (P **)almat_pointer(n,n);
790: for ( i = 0; i < n; i++ )
791: for ( j = 0; j < n; j++ )
792: ptomp(mod,mat[i][j],&w[i][j]);
793: detmp(CO,mod,w,n,&d);
794: mptop(d,rp);
1.23 noro 795: }
796: }
797:
1.24 noro 798: void Pinvmat(NODE arg,LIST *rp)
1.23 noro 799: {
800: MAT m,r;
801: int n,i,j,mod;
802: P dn;
803: P **mat,**imat,**w;
804: NODE nd;
805:
806: m = (MAT)ARG0(arg);
807: asir_assert(m,O_MAT,"invmat");
808: if ( m->row != m->col )
809: error("invmat : non-square matrix");
810: else if ( argc(arg) == 1 ) {
811: n = m->row;
812: invmatp(CO,(P **)BDY(m),n,&imat,&dn);
813: NEWMAT(r); r->row = n; r->col = n; r->body = (pointer **)imat;
814: nd = mknode(2,r,dn);
815: MKLIST(*rp,nd);
816: } else {
817: n = m->row; mod = QTOS((Q)ARG1(arg)); mat = (P **)BDY(m);
818: w = (P **)almat_pointer(n,n);
819: for ( i = 0; i < n; i++ )
820: for ( j = 0; j < n; j++ )
821: ptomp(mod,mat[i][j],&w[i][j]);
822: #if 0
823: detmp(CO,mod,w,n,&d);
824: mptop(d,rp);
825: #else
826: error("not implemented yet");
827: #endif
1.1 noro 828: }
1.25 noro 829: }
830:
831: /*
832: input : a row x col matrix A
833: A[I] <-> A[I][0]*x_0+A[I][1]*x_1+...
834:
1.45 noro 835: output : [B,D,R,C]
1.25 noro 836: B : a rank(A) x col-rank(A) matrix
1.45 noro 837: D : the denominator
1.25 noro 838: R : a vector of length rank(A)
839: C : a vector of length col-rank(A)
1.45 noro 840: B[I] <-> D*x_{R[I]}+B[I][0]x_{C[0]}+B[I][1]x_{C[1]}+...
1.25 noro 841: */
842:
843: void Pgeneric_gauss_elim(NODE arg,LIST *rp)
844: {
1.48 ! noro 845: NODE n0,opt,p;
1.25 noro 846: MAT m,nm;
847: int *ri,*ci;
848: VECT rind,cind;
849: Q dn,q;
850: int i,j,k,l,row,col,t,rank;
1.48 ! noro 851: int is_hensel = 0;
! 852: char *key;
! 853: Obj value;
! 854:
! 855: if ( current_option ) {
! 856: for ( opt = current_option; opt; opt = NEXT(opt) ) {
! 857: p = BDY((LIST)BDY(opt));
! 858: key = BDY((STRING)BDY(p));
! 859: value = (Obj)BDY(NEXT(p));
! 860: if ( !strcmp(key,"hensel") && value ) {
! 861: is_hensel = value ? 1 : 0;
! 862: break;
! 863: }
! 864: }
! 865: }
1.25 noro 866: asir_assert(ARG0(arg),O_MAT,"generic_gauss_elim");
867: m = (MAT)ARG0(arg);
868: row = m->row; col = m->col;
1.48 ! noro 869: if ( is_hensel )
! 870: rank = generic_gauss_elim_hensel(m,&nm,&dn,&ri,&ci);
! 871: else
! 872: rank = generic_gauss_elim(m,&nm,&dn,&ri,&ci);
1.25 noro 873: t = col-rank;
874: MKVECT(rind,rank);
875: MKVECT(cind,t);
876: for ( i = 0; i < rank; i++ ) {
877: STOQ(ri[i],q);
878: BDY(rind)[i] = (pointer)q;
879: }
880: for ( i = 0; i < t; i++ ) {
881: STOQ(ci[i],q);
882: BDY(cind)[i] = (pointer)q;
883: }
884: n0 = mknode(4,nm,dn,rind,cind);
885: MKLIST(*rp,n0);
1.1 noro 886: }
887:
888: /*
889: input : a row x col matrix A
890: A[I] <-> A[I][0]*x_0+A[I][1]*x_1+...
891:
892: output : [B,R,C]
893: B : a rank(A) x col-rank(A) matrix
894: R : a vector of length rank(A)
895: C : a vector of length col-rank(A)
1.47 noro 896: RN : a vector of length rank(A) indicating useful rows
897:
1.1 noro 898: B[I] <-> x_{R[I]}+B[I][0]x_{C[0]}+B[I][1]x_{C[1]}+...
899: */
900:
1.24 noro 901: void Pgeneric_gauss_elim_mod(NODE arg,LIST *rp)
1.1 noro 902: {
903: NODE n0;
904: MAT m,mat;
1.47 noro 905: VECT rind,cind,rnum;
1.1 noro 906: Q **tmat;
1.47 noro 907: int **wmat,**row0;
908: Q *rib,*cib,*rnb;
909: int *colstat,*p;
1.1 noro 910: Q q;
1.24 noro 911: int md,i,j,k,l,row,col,t,rank;
1.1 noro 912:
913: asir_assert(ARG0(arg),O_MAT,"generic_gauss_elim_mod");
914: asir_assert(ARG1(arg),O_N,"generic_gauss_elim_mod");
915: m = (MAT)ARG0(arg); md = QTOS((Q)ARG1(arg));
916: row = m->row; col = m->col; tmat = (Q **)m->body;
917: wmat = (int **)almat(row,col);
1.47 noro 918:
919: row0 = (int **)ALLOCA(row*sizeof(int *));
920: for ( i = 0; i < row; i++ ) row0[i] = wmat[i];
921:
1.1 noro 922: colstat = (int *)MALLOC_ATOMIC(col*sizeof(int));
923: for ( i = 0; i < row; i++ )
924: for ( j = 0; j < col; j++ )
925: if ( q = (Q)tmat[i][j] ) {
926: t = rem(NM(q),md);
927: if ( t && SGN(q) < 0 )
928: t = (md - t) % md;
929: wmat[i][j] = t;
930: } else
931: wmat[i][j] = 0;
932: rank = generic_gauss_elim_mod(wmat,row,col,md,colstat);
933:
1.47 noro 934: MKVECT(rnum,rank);
935: rnb = (Q *)rnum->body;
936: for ( i = 0; i < rank; i++ )
937: for ( j = 0, p = wmat[i]; j < row; j++ )
938: if ( p == row0[j] )
939: STOQ(j,rnb[i]);
940:
1.1 noro 941: MKMAT(mat,rank,col-rank);
942: tmat = (Q **)mat->body;
943: for ( i = 0; i < rank; i++ )
944: for ( j = k = 0; j < col; j++ )
945: if ( !colstat[j] ) {
946: UTOQ(wmat[i][j],tmat[i][k]); k++;
947: }
948:
949: MKVECT(rind,rank);
950: MKVECT(cind,col-rank);
951: rib = (Q *)rind->body; cib = (Q *)cind->body;
952: for ( j = k = l = 0; j < col; j++ )
953: if ( colstat[j] ) {
954: STOQ(j,rib[k]); k++;
955: } else {
956: STOQ(j,cib[l]); l++;
957: }
1.47 noro 958: n0 = mknode(4,mat,rind,cind,rnum);
1.1 noro 959: MKLIST(*rp,n0);
960: }
961:
1.24 noro 962: void Pleqm(NODE arg,VECT *rp)
1.1 noro 963: {
964: MAT m;
965: VECT vect;
966: pointer **mat;
967: Q *v;
968: Q q;
969: int **wmat;
970: int md,i,j,row,col,t,n,status;
971:
972: asir_assert(ARG0(arg),O_MAT,"leqm");
973: asir_assert(ARG1(arg),O_N,"leqm");
974: m = (MAT)ARG0(arg); md = QTOS((Q)ARG1(arg));
975: row = m->row; col = m->col; mat = m->body;
976: wmat = (int **)almat(row,col);
977: for ( i = 0; i < row; i++ )
978: for ( j = 0; j < col; j++ )
979: if ( q = (Q)mat[i][j] ) {
980: t = rem(NM(q),md);
981: if ( SGN(q) < 0 )
982: t = (md - t) % md;
983: wmat[i][j] = t;
984: } else
985: wmat[i][j] = 0;
986: status = gauss_elim_mod(wmat,row,col,md);
987: if ( status < 0 )
988: *rp = 0;
989: else if ( status > 0 )
990: *rp = (VECT)ONE;
991: else {
992: n = col - 1;
993: MKVECT(vect,n);
994: for ( i = 0, v = (Q *)vect->body; i < n; i++ ) {
995: t = (md-wmat[i][n])%md; STOQ(t,v[i]);
996: }
997: *rp = vect;
998: }
999: }
1000:
1.24 noro 1001: int gauss_elim_mod(int **mat,int row,int col,int md)
1.1 noro 1002: {
1003: int i,j,k,inv,a,n;
1004: int *t,*pivot;
1005:
1006: n = col - 1;
1007: for ( j = 0; j < n; j++ ) {
1008: for ( i = j; i < row && !mat[i][j]; i++ );
1009: if ( i == row )
1010: return 1;
1011: if ( i != j ) {
1012: t = mat[i]; mat[i] = mat[j]; mat[j] = t;
1013: }
1014: pivot = mat[j];
1015: inv = invm(pivot[j],md);
1016: for ( k = j; k <= n; k++ ) {
1017: /* pivot[k] = dmar(pivot[k],inv,0,md); */
1018: DMAR(pivot[k],inv,0,md,pivot[k])
1019: }
1020: for ( i = 0; i < row; i++ ) {
1021: t = mat[i];
1022: if ( i != j && (a = t[j]) )
1023: for ( k = j, a = md - a; k <= n; k++ ) {
1.8 noro 1024: unsigned int tk;
1.1 noro 1025: /* t[k] = dmar(pivot[k],a,t[k],md); */
1.8 noro 1026: DMAR(pivot[k],a,t[k],md,tk)
1027: t[k] = tk;
1.1 noro 1028: }
1029: }
1030: }
1031: for ( i = n; i < row && !mat[i][n]; i++ );
1032: if ( i == row )
1033: return 0;
1034: else
1035: return -1;
1036: }
1037:
1.4 noro 1038: struct oEGT eg_mod,eg_elim,eg_elim1,eg_elim2,eg_chrem,eg_gschk,eg_intrat,eg_symb;
1.31 noro 1039: struct oEGT eg_conv;
1.1 noro 1040:
1.24 noro 1041: int generic_gauss_elim(MAT mat,MAT *nm,Q *dn,int **rindp,int **cindp)
1.1 noro 1042: {
1043: int **wmat;
1044: Q **bmat;
1045: N **tmat;
1046: Q *bmi;
1047: N *tmi;
1048: Q q;
1049: int *wmi;
1050: int *colstat,*wcolstat,*rind,*cind;
1051: int row,col,ind,md,i,j,k,l,t,t1,rank,rank0,inv;
1052: N m1,m2,m3,s,u;
1053: MAT r,crmat;
1054: struct oEGT tmp0,tmp1;
1055: struct oEGT eg_mod_split,eg_elim_split,eg_chrem_split;
1056: struct oEGT eg_intrat_split,eg_gschk_split;
1057: int ret;
1058:
1059: init_eg(&eg_mod_split); init_eg(&eg_chrem_split);
1060: init_eg(&eg_elim_split); init_eg(&eg_intrat_split);
1061: init_eg(&eg_gschk_split);
1062: bmat = (Q **)mat->body;
1063: row = mat->row; col = mat->col;
1064: wmat = (int **)almat(row,col);
1065: colstat = (int *)MALLOC_ATOMIC(col*sizeof(int));
1066: wcolstat = (int *)MALLOC_ATOMIC(col*sizeof(int));
1067: for ( ind = 0; ; ind++ ) {
1.11 noro 1068: if ( DP_Print ) {
1.2 noro 1069: fprintf(asir_out,"."); fflush(asir_out);
1070: }
1.12 noro 1071: md = get_lprime(ind);
1.1 noro 1072: get_eg(&tmp0);
1073: for ( i = 0; i < row; i++ )
1074: for ( j = 0, bmi = bmat[i], wmi = wmat[i]; j < col; j++ )
1075: if ( q = (Q)bmi[j] ) {
1076: t = rem(NM(q),md);
1077: if ( t && SGN(q) < 0 )
1078: t = (md - t) % md;
1079: wmi[j] = t;
1080: } else
1081: wmi[j] = 0;
1082: get_eg(&tmp1);
1083: add_eg(&eg_mod,&tmp0,&tmp1);
1084: add_eg(&eg_mod_split,&tmp0,&tmp1);
1085: get_eg(&tmp0);
1086: rank = generic_gauss_elim_mod(wmat,row,col,md,wcolstat);
1087: get_eg(&tmp1);
1088: add_eg(&eg_elim,&tmp0,&tmp1);
1089: add_eg(&eg_elim_split,&tmp0,&tmp1);
1090: if ( !ind ) {
1091: RESET:
1092: UTON(md,m1);
1093: rank0 = rank;
1094: bcopy(wcolstat,colstat,col*sizeof(int));
1095: MKMAT(crmat,rank,col-rank);
1096: MKMAT(r,rank,col-rank); *nm = r;
1097: tmat = (N **)crmat->body;
1098: for ( i = 0; i < rank; i++ )
1099: for ( j = k = 0, tmi = tmat[i], wmi = wmat[i]; j < col; j++ )
1100: if ( !colstat[j] ) {
1101: UTON(wmi[j],tmi[k]); k++;
1102: }
1103: } else {
1104: if ( rank < rank0 ) {
1.11 noro 1105: if ( DP_Print ) {
1.1 noro 1106: fprintf(asir_out,"lower rank matrix; continuing...\n");
1.2 noro 1107: fflush(asir_out);
1108: }
1.1 noro 1109: continue;
1110: } else if ( rank > rank0 ) {
1.11 noro 1111: if ( DP_Print ) {
1.1 noro 1112: fprintf(asir_out,"higher rank matrix; resetting...\n");
1.2 noro 1113: fflush(asir_out);
1114: }
1.1 noro 1115: goto RESET;
1116: } else {
1117: for ( j = 0; (j<col) && (colstat[j]==wcolstat[j]); j++ );
1118: if ( j < col ) {
1.11 noro 1119: if ( DP_Print ) {
1.1 noro 1120: fprintf(asir_out,"inconsitent colstat; resetting...\n");
1.2 noro 1121: fflush(asir_out);
1122: }
1.1 noro 1123: goto RESET;
1124: }
1125: }
1126:
1127: get_eg(&tmp0);
1128: inv = invm(rem(m1,md),md);
1129: UTON(md,m2); muln(m1,m2,&m3);
1130: for ( i = 0; i < rank; i++ )
1131: for ( j = k = 0, tmi = tmat[i], wmi = wmat[i]; j < col; j++ )
1132: if ( !colstat[j] ) {
1133: if ( tmi[k] ) {
1134: /* f3 = f1+m1*(m1 mod m2)^(-1)*(f2 - f1 mod m2) */
1135: t = rem(tmi[k],md);
1136: if ( wmi[j] >= t )
1137: t = wmi[j]-t;
1138: else
1139: t = md-(t-wmi[j]);
1140: DMAR(t,inv,0,md,t1)
1141: UTON(t1,u);
1142: muln(m1,u,&s);
1143: addn(tmi[k],s,&u); tmi[k] = u;
1144: } else if ( wmi[j] ) {
1145: /* f3 = m1*(m1 mod m2)^(-1)*f2 */
1146: DMAR(wmi[j],inv,0,md,t)
1147: UTON(t,u);
1148: muln(m1,u,&s); tmi[k] = s;
1149: }
1150: k++;
1151: }
1152: m1 = m3;
1153: get_eg(&tmp1);
1154: add_eg(&eg_chrem,&tmp0,&tmp1);
1155: add_eg(&eg_chrem_split,&tmp0,&tmp1);
1156:
1157: get_eg(&tmp0);
1.38 noro 1158: if ( ind % F4_INTRAT_PERIOD )
1.13 noro 1159: ret = 0;
1160: else
1161: ret = intmtoratm(crmat,m1,*nm,dn);
1.1 noro 1162: get_eg(&tmp1);
1163: add_eg(&eg_intrat,&tmp0,&tmp1);
1164: add_eg(&eg_intrat_split,&tmp0,&tmp1);
1165: if ( ret ) {
1166: *rindp = rind = (int *)MALLOC_ATOMIC(rank*sizeof(int));
1167: *cindp = cind = (int *)MALLOC_ATOMIC((col-rank)*sizeof(int));
1168: for ( j = k = l = 0; j < col; j++ )
1169: if ( colstat[j] )
1170: rind[k++] = j;
1171: else
1172: cind[l++] = j;
1173: get_eg(&tmp0);
1.3 noro 1174: if ( gensolve_check(mat,*nm,*dn,rind,cind) ) {
1175: get_eg(&tmp1);
1176: add_eg(&eg_gschk,&tmp0,&tmp1);
1177: add_eg(&eg_gschk_split,&tmp0,&tmp1);
1.11 noro 1178: if ( DP_Print ) {
1.3 noro 1179: print_eg("Mod",&eg_mod_split);
1180: print_eg("Elim",&eg_elim_split);
1181: print_eg("ChRem",&eg_chrem_split);
1182: print_eg("IntRat",&eg_intrat_split);
1183: print_eg("Check",&eg_gschk_split);
1184: fflush(asir_out);
1185: }
1186: return rank;
1187: }
1188: }
1189: }
1190: }
1191: }
1192:
1.24 noro 1193: int generic_gauss_elim_hensel(MAT mat,MAT *nmmat,Q *dn,int **rindp,int **cindp)
1.3 noro 1194: {
1195: MAT bmat,xmat;
1196: Q **a0,**a,**b,**x,**nm;
1197: Q *ai,*bi,*xi;
1198: int row,col;
1199: int **w;
1200: int *wi;
1201: int **wc;
1202: Q mdq,q,s,u;
1203: N tn;
1204: int ind,md,i,j,k,l,li,ri,rank;
1205: unsigned int t;
1206: int *cinfo,*rinfo;
1207: int *rind,*cind;
1208: int count;
1.41 noro 1209: int ret;
1210: struct oEGT eg_mul,eg_inv,eg_intrat,eg_check,tmp0,tmp1;
1.39 noro 1211: int period;
1.44 noro 1212: int *wx,*ptr;
1213: int wxsize,nsize;
1214: N wn;
1215: Q wq;
1.3 noro 1216:
1217: a0 = (Q **)mat->body;
1218: row = mat->row; col = mat->col;
1219: w = (int **)almat(row,col);
1220: for ( ind = 0; ; ind++ ) {
1.12 noro 1221: md = get_lprime(ind);
1.3 noro 1222: STOQ(md,mdq);
1223: for ( i = 0; i < row; i++ )
1224: for ( j = 0, ai = a0[i], wi = w[i]; j < col; j++ )
1225: if ( q = (Q)ai[j] ) {
1226: t = rem(NM(q),md);
1227: if ( t && SGN(q) < 0 )
1228: t = (md - t) % md;
1229: wi[j] = t;
1230: } else
1231: wi[j] = 0;
1232:
1.48 ! noro 1233: if ( DP_Print ) {
! 1234: fprintf(asir_out,"LU decomposition.."); fflush(asir_out);
! 1235: }
1.27 noro 1236: rank = find_lhs_and_lu_mod((unsigned int **)w,row,col,md,&rinfo,&cinfo);
1.48 ! noro 1237: if ( DP_Print ) {
! 1238: fprintf(asir_out,"done.\n"); fflush(asir_out);
! 1239: }
1.3 noro 1240: a = (Q **)almat_pointer(rank,rank); /* lhs mat */
1241: MKMAT(bmat,rank,col-rank); b = (Q **)bmat->body; /* lhs mat */
1242: for ( j = li = ri = 0; j < col; j++ )
1243: if ( cinfo[j] ) {
1244: /* the column is in lhs */
1245: for ( i = 0; i < rank; i++ ) {
1246: w[i][li] = w[i][j];
1247: a[i][li] = a0[rinfo[i]][j];
1248: }
1249: li++;
1250: } else {
1251: /* the column is in rhs */
1252: for ( i = 0; i < rank; i++ )
1253: b[i][ri] = a0[rinfo[i]][j];
1254: ri++;
1255: }
1256:
1257: /* solve Ax+B=0; A: rank x rank, B: rank x ri */
1258: MKMAT(xmat,rank,ri); x = (Q **)(xmat)->body;
1259: MKMAT(*nmmat,rank,ri); nm = (Q **)(*nmmat)->body;
1260: /* use the right part of w as work area */
1261: /* ri = col - rank */
1262: wc = (int **)almat(rank,ri);
1263: for ( i = 0; i < rank; i++ )
1264: wc[i] = w[i]+rank;
1265: *rindp = rind = (int *)MALLOC_ATOMIC(rank*sizeof(int));
1266: *cindp = cind = (int *)MALLOC_ATOMIC((ri)*sizeof(int));
1267:
1268: init_eg(&eg_mul); init_eg(&eg_inv);
1.41 noro 1269: init_eg(&eg_check); init_eg(&eg_intrat);
1.39 noro 1270: period = F4_INTRAT_PERIOD;
1.44 noro 1271: nsize = period;
1272: wxsize = rank*ri*nsize;
1273: wx = (int *)MALLOC_ATOMIC(wxsize*sizeof(int));
1274: for ( i = 0; i < wxsize; i++ ) wx[i] = 0;
1275: for ( q = ONE, count = 0; ; ) {
1.48 ! noro 1276: if ( DP_Print )
1.41 noro 1277: fprintf(stderr,"o");
1.3 noro 1278: /* wc = -b mod md */
1.44 noro 1279: get_eg(&tmp0);
1.3 noro 1280: for ( i = 0; i < rank; i++ )
1281: for ( j = 0, bi = b[i], wi = wc[i]; j < ri; j++ )
1282: if ( u = (Q)bi[j] ) {
1283: t = rem(NM(u),md);
1284: if ( t && SGN(u) > 0 )
1285: t = (md - t) % md;
1286: wi[j] = t;
1287: } else
1288: wi[j] = 0;
1.44 noro 1289: /* wc = A^(-1)wc; wc is not normalized */
1290: solve_by_lu_mod(w,rank,md,wc,ri,0);
1291: /* wx += q*wc */
1292: ptr = wx;
1293: for ( i = 0; i < rank; i++ )
1294: for ( j = 0, wi = wc[i]; j < ri; j++ ) {
1295: if ( wi[j] )
1296: muln_1(BD(NM(q)),PL(NM(q)),wi[j],ptr);
1297: ptr += nsize;
1298: }
1299: count++;
1.1 noro 1300: get_eg(&tmp1);
1.3 noro 1301: add_eg(&eg_inv,&tmp0,&tmp1);
1302: get_eg(&tmp0);
1303: for ( i = 0; i < rank; i++ )
1304: for ( j = 0; j < ri; j++ ) {
1305: inner_product_mat_int_mod(a,wc,rank,i,j,&u);
1306: addq(b[i][j],u,&s);
1307: if ( s ) {
1308: t = divin(NM(s),md,&tn);
1309: if ( t )
1310: error("generic_gauss_elim_hensel:incosistent");
1311: NTOQ(tn,SGN(s),b[i][j]);
1312: } else
1313: b[i][j] = 0;
1314: }
1315: get_eg(&tmp1);
1316: add_eg(&eg_mul,&tmp0,&tmp1);
1317: /* q = q*md */
1318: mulq(q,mdq,&u); q = u;
1.44 noro 1319: if ( count == period ) {
1.41 noro 1320: get_eg(&tmp0);
1.44 noro 1321: ptr = wx;
1322: for ( i = 0; i < rank; i++ )
1323: for ( j = 0, xi = x[i]; j < ri;
1324: j++, ptr += nsize ) {
1325: for ( k = nsize-1; k >= 0 && !ptr[k]; k-- );
1326: if ( k >= 0 ) {
1327: wn = NALLOC(k+1);
1328: PL(wn) = k+1;
1329: for ( l = 0; l <= k; l++ ) BD(wn)[l] = (unsigned int)ptr[l];
1330: NTOQ(wn,1,wq);
1331: subq(xi[j],wq,&u); xi[j] = u;
1332: }
1333: }
1.41 noro 1334: ret = intmtoratm_q(xmat,NM(q),*nmmat,dn);
1335: get_eg(&tmp1); add_eg(&eg_intrat,&tmp0,&tmp1);
1336: if ( ret ) {
1.39 noro 1337: for ( j = k = l = 0; j < col; j++ )
1338: if ( cinfo[j] )
1339: rind[k++] = j;
1340: else
1341: cind[l++] = j;
1.41 noro 1342: get_eg(&tmp0);
1343: ret = gensolve_check(mat,*nmmat,*dn,rind,cind);
1344: get_eg(&tmp1); add_eg(&eg_check,&tmp0,&tmp1);
1345: if ( ret ) {
1.42 noro 1346: if ( DP_Print > 3 ) {
1.40 noro 1347: fprintf(stderr,"\n");
1348: print_eg("INV",&eg_inv);
1349: print_eg("MUL",&eg_mul);
1.41 noro 1350: print_eg("INTRAT",&eg_intrat);
1351: print_eg("CHECK",&eg_check);
1.40 noro 1352: fflush(asir_out);
1353: }
1.39 noro 1354: return rank;
1355: }
1.44 noro 1356: } else {
1357: period = period*3/2;
1358: count = 0;
1359: nsize += period;
1360: wxsize += rank*ri*nsize;
1361: wx = (int *)REALLOC(wx,wxsize*sizeof(int));
1362: for ( i = 0; i < wxsize; i++ ) wx[i] = 0;
1363: }
1.41 noro 1364: }
1.1 noro 1365: }
1366: }
1367: }
1368:
1369: int f4_nocheck;
1370:
1.24 noro 1371: int gensolve_check(MAT mat,MAT nm,Q dn,int *rind,int *cind)
1.1 noro 1372: {
1373: int row,col,rank,clen,i,j,k,l;
1.24 noro 1374: Q s,t;
1.1 noro 1375: Q *w;
1376: Q *mati,*nmk;
1377:
1378: if ( f4_nocheck )
1379: return 1;
1380: row = mat->row; col = mat->col;
1381: rank = nm->row; clen = nm->col;
1382: w = (Q *)MALLOC(clen*sizeof(Q));
1383: for ( i = 0; i < row; i++ ) {
1384: mati = (Q *)mat->body[i];
1385: #if 1
1386: bzero(w,clen*sizeof(Q));
1387: for ( k = 0; k < rank; k++ )
1388: for ( l = 0, nmk = (Q *)nm->body[k]; l < clen; l++ ) {
1389: mulq(mati[rind[k]],nmk[l],&t);
1390: addq(w[l],t,&s); w[l] = s;
1391: }
1392: for ( j = 0; j < clen; j++ ) {
1393: mulq(dn,mati[cind[j]],&t);
1394: if ( cmpq(w[j],t) )
1395: break;
1396: }
1397: #else
1398: for ( j = 0; j < clen; j++ ) {
1399: for ( k = 0, s = 0; k < rank; k++ ) {
1400: mulq(mati[rind[k]],nm->body[k][j],&t);
1401: addq(s,t,&u); s = u;
1402: }
1403: mulq(dn,mati[cind[j]],&t);
1404: if ( cmpq(s,t) )
1405: break;
1406: }
1407: #endif
1408: if ( j != clen )
1409: break;
1410: }
1411: if ( i != row )
1412: return 0;
1413: else
1414: return 1;
1415: }
1416:
1417: /* assuming 0 < c < m */
1418:
1.24 noro 1419: int inttorat(N c,N m,N b,int *sgnp,N *nmp,N *dnp)
1.1 noro 1420: {
1.24 noro 1421: Q qq,t,u1,v1,r1;
1422: N q,u2,v2,r2;
1.1 noro 1423:
1424: u1 = 0; v1 = ONE; u2 = m; v2 = c;
1425: while ( cmpn(v2,b) >= 0 ) {
1426: divn(u2,v2,&q,&r2); u2 = v2; v2 = r2;
1427: NTOQ(q,1,qq); mulq(qq,v1,&t); subq(u1,t,&r1); u1 = v1; v1 = r1;
1428: }
1429: if ( cmpn(NM(v1),b) >= 0 )
1430: return 0;
1431: else {
1432: *nmp = v2;
1433: *dnp = NM(v1);
1434: *sgnp = SGN(v1);
1435: return 1;
1436: }
1437: }
1438:
1439: /* mat->body = N ** */
1440:
1.24 noro 1441: int intmtoratm(MAT mat,N md,MAT nm,Q *dn)
1.1 noro 1442: {
1443: N t,s,b;
1.24 noro 1444: Q dn0,dn1,nm1,q;
1.1 noro 1445: int i,j,k,l,row,col;
1446: Q **rmat;
1447: N **tmat;
1448: N *tmi;
1449: Q *nmk;
1450: N u,unm,udn;
1451: int sgn,ret;
1452:
1.3 noro 1453: if ( UNIN(md) )
1454: return 0;
1.1 noro 1455: row = mat->row; col = mat->col;
1456: bshiftn(md,1,&t);
1457: isqrt(t,&s);
1458: bshiftn(s,64,&b);
1459: if ( !b )
1460: b = ONEN;
1461: dn0 = ONE;
1462: tmat = (N **)mat->body;
1463: rmat = (Q **)nm->body;
1464: for ( i = 0; i < row; i++ )
1465: for ( j = 0, tmi = tmat[i]; j < col; j++ )
1466: if ( tmi[j] ) {
1467: muln(tmi[j],NM(dn0),&s);
1468: remn(s,md,&u);
1469: ret = inttorat(u,md,b,&sgn,&unm,&udn);
1470: if ( !ret )
1471: return 0;
1472: else {
1473: NTOQ(unm,sgn,nm1);
1474: NTOQ(udn,1,dn1);
1475: if ( !UNIQ(dn1) ) {
1476: for ( k = 0; k < i; k++ )
1477: for ( l = 0, nmk = rmat[k]; l < col; l++ ) {
1478: mulq(nmk[l],dn1,&q); nmk[l] = q;
1479: }
1480: for ( l = 0, nmk = rmat[i]; l < j; l++ ) {
1481: mulq(nmk[l],dn1,&q); nmk[l] = q;
1482: }
1483: }
1484: rmat[i][j] = nm1;
1485: mulq(dn0,dn1,&q); dn0 = q;
1486: }
1487: }
1488: *dn = dn0;
1489: return 1;
1490: }
1491:
1.3 noro 1492: /* mat->body = Q ** */
1493:
1.24 noro 1494: int intmtoratm_q(MAT mat,N md,MAT nm,Q *dn)
1.3 noro 1495: {
1496: N t,s,b;
1.24 noro 1497: Q dn0,dn1,nm1,q;
1.3 noro 1498: int i,j,k,l,row,col;
1499: Q **rmat;
1500: Q **tmat;
1501: Q *tmi;
1502: Q *nmk;
1503: N u,unm,udn;
1504: int sgn,ret;
1505:
1506: if ( UNIN(md) )
1507: return 0;
1508: row = mat->row; col = mat->col;
1509: bshiftn(md,1,&t);
1510: isqrt(t,&s);
1511: bshiftn(s,64,&b);
1512: if ( !b )
1513: b = ONEN;
1514: dn0 = ONE;
1515: tmat = (Q **)mat->body;
1516: rmat = (Q **)nm->body;
1517: for ( i = 0; i < row; i++ )
1518: for ( j = 0, tmi = tmat[i]; j < col; j++ )
1519: if ( tmi[j] ) {
1520: muln(NM(tmi[j]),NM(dn0),&s);
1521: remn(s,md,&u);
1522: ret = inttorat(u,md,b,&sgn,&unm,&udn);
1523: if ( !ret )
1524: return 0;
1525: else {
1526: if ( SGN(tmi[j])<0 )
1527: sgn = -sgn;
1528: NTOQ(unm,sgn,nm1);
1529: NTOQ(udn,1,dn1);
1530: if ( !UNIQ(dn1) ) {
1531: for ( k = 0; k < i; k++ )
1532: for ( l = 0, nmk = rmat[k]; l < col; l++ ) {
1533: mulq(nmk[l],dn1,&q); nmk[l] = q;
1534: }
1535: for ( l = 0, nmk = rmat[i]; l < j; l++ ) {
1536: mulq(nmk[l],dn1,&q); nmk[l] = q;
1537: }
1538: }
1539: rmat[i][j] = nm1;
1540: mulq(dn0,dn1,&q); dn0 = q;
1541: }
1542: }
1543: *dn = dn0;
1544: return 1;
1545: }
1546:
1.4 noro 1547: #define ONE_STEP1 if ( zzz = *s ) { DMAR(zzz,hc,*tj,md,*tj) } tj++; s++;
1548:
1.24 noro 1549: void reduce_reducers_mod(int **mat,int row,int col,int md)
1.4 noro 1550: {
1551: int i,j,k,l,hc,zzz;
1552: int *t,*s,*tj,*ind;
1553:
1554: /* reduce the reducers */
1555: ind = (int *)ALLOCA(row*sizeof(int));
1556: for ( i = 0; i < row; i++ ) {
1557: t = mat[i];
1558: for ( j = 0; j < col && !t[j]; j++ );
1559: /* register the position of the head term */
1560: ind[i] = j;
1561: for ( l = i-1; l >= 0; l-- ) {
1562: /* reduce mat[i] by mat[l] */
1563: if ( hc = t[ind[l]] ) {
1564: /* mat[i] = mat[i]-hc*mat[l] */
1565: j = ind[l];
1566: s = mat[l]+j;
1567: tj = t+j;
1568: hc = md-hc;
1569: k = col-j;
1570: for ( ; k >= 64; k -= 64 ) {
1571: ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
1572: ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
1573: ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
1574: ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
1575: ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
1576: ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
1577: ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
1578: ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
1579: ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
1580: ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
1581: ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
1582: ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
1583: ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
1584: ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
1585: ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
1586: ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
1587: }
1.16 noro 1588: for ( ; k > 0; k-- ) {
1.4 noro 1589: if ( zzz = *s ) { DMAR(zzz,hc,*tj,md,*tj) } tj++; s++;
1590: }
1591: }
1592: }
1593: }
1594: }
1595:
1596: /*
1597: mat[i] : reducers (i=0,...,nred-1)
1598: spolys (i=nred,...,row-1)
1599: mat[0] < mat[1] < ... < mat[nred-1] w.r.t the term order
1600: 1. reduce the reducers
1601: 2. reduce spolys by the reduced reducers
1602: */
1603:
1.24 noro 1604: void pre_reduce_mod(int **mat,int row,int col,int nred,int md)
1.4 noro 1605: {
1606: int i,j,k,l,hc,inv;
1607: int *t,*s,*tk,*ind;
1608:
1609: #if 1
1610: /* reduce the reducers */
1611: ind = (int *)ALLOCA(row*sizeof(int));
1612: for ( i = 0; i < nred; i++ ) {
1613: /* make mat[i] monic and mat[i] by mat[0],...,mat[i-1] */
1614: t = mat[i];
1615: for ( j = 0; j < col && !t[j]; j++ );
1616: /* register the position of the head term */
1617: ind[i] = j;
1618: inv = invm(t[j],md);
1619: for ( k = j; k < col; k++ )
1620: if ( t[k] )
1621: DMAR(t[k],inv,0,md,t[k])
1622: for ( l = i-1; l >= 0; l-- ) {
1623: /* reduce mat[i] by mat[l] */
1624: if ( hc = t[ind[l]] ) {
1625: /* mat[i] = mat[i]-hc*mat[l] */
1626: for ( k = ind[l], hc = md-hc, s = mat[l]+k, tk = t+k;
1627: k < col; k++, tk++, s++ )
1628: if ( *s )
1629: DMAR(*s,hc,*tk,md,*tk)
1630: }
1631: }
1632: }
1633: /* reduce the spolys */
1634: for ( i = nred; i < row; i++ ) {
1635: t = mat[i];
1636: for ( l = nred-1; l >= 0; l-- ) {
1637: /* reduce mat[i] by mat[l] */
1638: if ( hc = t[ind[l]] ) {
1639: /* mat[i] = mat[i]-hc*mat[l] */
1640: for ( k = ind[l], hc = md-hc, s = mat[l]+k, tk = t+k;
1641: k < col; k++, tk++, s++ )
1642: if ( *s )
1643: DMAR(*s,hc,*tk,md,*tk)
1644: }
1645: }
1646: }
1647: #endif
1648: }
1649: /*
1650: mat[i] : reducers (i=0,...,nred-1)
1651: mat[0] < mat[1] < ... < mat[nred-1] w.r.t the term order
1652: */
1653:
1.24 noro 1654: void reduce_sp_by_red_mod(int *sp,int **redmat,int *ind,int nred,int col,int md)
1.4 noro 1655: {
1656: int i,j,k,hc,zzz;
1.24 noro 1657: int *s,*tj;
1.4 noro 1658:
1659: /* reduce the spolys by redmat */
1660: for ( i = nred-1; i >= 0; i-- ) {
1661: /* reduce sp by redmat[i] */
1662: if ( hc = sp[ind[i]] ) {
1663: /* sp = sp-hc*redmat[i] */
1664: j = ind[i];
1665: hc = md-hc;
1666: s = redmat[i]+j;
1667: tj = sp+j;
1.16 noro 1668: for ( k = col-j; k > 0; k-- ) {
1.4 noro 1669: if ( zzz = *s ) { DMAR(zzz,hc,*tj,md,*tj) } tj++; s++;
1.15 noro 1670: }
1671: }
1.17 noro 1672: }
1673: }
1674:
1675: /*
1.15 noro 1676: mat[i] : compressed reducers (i=0,...,nred-1)
1677: mat[0] < mat[1] < ... < mat[nred-1] w.r.t the term order
1678: */
1679:
1.24 noro 1680: void red_by_compress(int m,unsigned int *p,unsigned int *r,
1681: unsigned int *ri,unsigned int hc,int len)
1.18 noro 1682: {
1.19 noro 1683: unsigned int up,lo;
1.18 noro 1684: unsigned int dmy;
1685: unsigned int *pj;
1686:
1.21 noro 1687: p[*ri] = 0; r++; ri++;
1688: for ( len--; len; len--, r++, ri++ ) {
1689: pj = p+ *ri;
1690: DMA(*r,hc,*pj,up,lo);
1.18 noro 1691: if ( up ) {
1692: DSAB(m,up,lo,dmy,*pj);
1693: } else
1694: *pj = lo;
1695: }
1696: }
1697:
1698: /* p -= hc*r */
1699:
1.24 noro 1700: void red_by_vect(int m,unsigned int *p,unsigned int *r,unsigned int hc,int len)
1.18 noro 1701: {
1702: register unsigned int up,lo;
1703: unsigned int dmy;
1704:
1705: *p++ = 0; r++; len--;
1706: for ( ; len; len--, r++, p++ )
1707: if ( *r ) {
1.20 noro 1708: DMA(*r,hc,*p,up,lo);
1.18 noro 1709: if ( up ) {
1710: DSAB(m,up,lo,dmy,*p);
1711: } else
1712: *p = lo;
1713: }
1714: }
1715:
1.32 noro 1716: void red_by_vect_sf(int m,unsigned int *p,unsigned int *r,unsigned int hc,int len)
1717: {
1718: *p++ = 0; r++; len--;
1719: for ( ; len; len--, r++, p++ )
1720: if ( *r )
1721: *p = _addsf(_mulsf(*r,hc),*p);
1722: }
1723:
1.21 noro 1724: extern unsigned int **psca;
1725:
1.24 noro 1726: void reduce_sp_by_red_mod_compress (int *sp,CDP *redmat,int *ind,
1727: int nred,int col,int md)
1.15 noro 1728: {
1.24 noro 1729: int i,len;
1.15 noro 1730: CDP ri;
1.24 noro 1731: unsigned int hc;
1.18 noro 1732: unsigned int *usp;
1.15 noro 1733:
1.18 noro 1734: usp = (unsigned int *)sp;
1.15 noro 1735: /* reduce the spolys by redmat */
1736: for ( i = nred-1; i >= 0; i-- ) {
1737: /* reduce sp by redmat[i] */
1.18 noro 1738: usp[ind[i]] %= md;
1739: if ( hc = usp[ind[i]] ) {
1.15 noro 1740: /* sp = sp-hc*redmat[i] */
1741: hc = md-hc;
1742: ri = redmat[i];
1743: len = ri->len;
1.21 noro 1744: red_by_compress(md,usp,psca[ri->psindex],ri->body,hc,len);
1.4 noro 1745: }
1746: }
1.18 noro 1747: for ( i = 0; i < col; i++ )
1.24 noro 1748: if ( usp[i] >= (unsigned int)md )
1.18 noro 1749: usp[i] %= md;
1.4 noro 1750: }
1751:
1752: #define ONE_STEP2 if ( zzz = *pk ) { DMAR(zzz,a,*tk,md,*tk) } pk++; tk++;
1753:
1.24 noro 1754: int generic_gauss_elim_mod(int **mat0,int row,int col,int md,int *colstat)
1.1 noro 1755: {
1.24 noro 1756: int i,j,k,l,inv,a,rank;
1757: unsigned int *t,*pivot,*pk;
1.18 noro 1758: unsigned int **mat;
1.1 noro 1759:
1.18 noro 1760: mat = (unsigned int **)mat0;
1.1 noro 1761: for ( rank = 0, j = 0; j < col; j++ ) {
1.18 noro 1762: for ( i = rank; i < row; i++ )
1763: mat[i][j] %= md;
1764: for ( i = rank; i < row; i++ )
1765: if ( mat[i][j] )
1766: break;
1.1 noro 1767: if ( i == row ) {
1768: colstat[j] = 0;
1769: continue;
1770: } else
1771: colstat[j] = 1;
1772: if ( i != rank ) {
1773: t = mat[i]; mat[i] = mat[rank]; mat[rank] = t;
1774: }
1775: pivot = mat[rank];
1776: inv = invm(pivot[j],md);
1.4 noro 1777: for ( k = j, pk = pivot+k; k < col; k++, pk++ )
1778: if ( *pk ) {
1.24 noro 1779: if ( *pk >= (unsigned int)md )
1.18 noro 1780: *pk %= md;
1.4 noro 1781: DMAR(*pk,inv,0,md,*pk)
1.1 noro 1782: }
1783: for ( i = rank+1; i < row; i++ ) {
1784: t = mat[i];
1.18 noro 1785: if ( a = t[j] )
1786: red_by_vect(md,t+j,pivot+j,md-a,col-j);
1.1 noro 1787: }
1788: rank++;
1789: }
1790: for ( j = col-1, l = rank-1; j >= 0; j-- )
1791: if ( colstat[j] ) {
1792: pivot = mat[l];
1793: for ( i = 0; i < l; i++ ) {
1794: t = mat[i];
1.18 noro 1795: t[j] %= md;
1796: if ( a = t[j] )
1797: red_by_vect(md,t+j,pivot+j,md-a,col-j);
1.1 noro 1798: }
1799: l--;
1.18 noro 1800: }
1801: for ( j = 0, l = 0; l < rank; j++ )
1802: if ( colstat[j] ) {
1803: t = mat[l];
1804: for ( k = j; k < col; k++ )
1.24 noro 1805: if ( t[k] >= (unsigned int)md )
1.18 noro 1806: t[k] %= md;
1807: l++;
1.32 noro 1808: }
1809: return rank;
1810: }
1811:
1812: int generic_gauss_elim_sf(int **mat0,int row,int col,int md,int *colstat)
1813: {
1814: int i,j,k,l,inv,a,rank;
1815: unsigned int *t,*pivot,*pk;
1816: unsigned int **mat;
1817:
1818: mat = (unsigned int **)mat0;
1819: for ( rank = 0, j = 0; j < col; j++ ) {
1820: for ( i = rank; i < row; i++ )
1821: if ( mat[i][j] )
1822: break;
1823: if ( i == row ) {
1824: colstat[j] = 0;
1825: continue;
1826: } else
1827: colstat[j] = 1;
1828: if ( i != rank ) {
1829: t = mat[i]; mat[i] = mat[rank]; mat[rank] = t;
1830: }
1831: pivot = mat[rank];
1832: inv = _invsf(pivot[j]);
1833: for ( k = j, pk = pivot+k; k < col; k++, pk++ )
1834: if ( *pk )
1835: *pk = _mulsf(*pk,inv);
1836: for ( i = rank+1; i < row; i++ ) {
1837: t = mat[i];
1838: if ( a = t[j] )
1839: red_by_vect_sf(md,t+j,pivot+j,_chsgnsf(a),col-j);
1840: }
1841: rank++;
1842: }
1843: for ( j = col-1, l = rank-1; j >= 0; j-- )
1844: if ( colstat[j] ) {
1845: pivot = mat[l];
1846: for ( i = 0; i < l; i++ ) {
1847: t = mat[i];
1848: if ( a = t[j] )
1849: red_by_vect_sf(md,t+j,pivot+j,_chsgnsf(a),col-j);
1850: }
1851: l--;
1.1 noro 1852: }
1853: return rank;
1854: }
1855:
1856: /* LU decomposition; a[i][i] = 1/U[i][i] */
1857:
1.24 noro 1858: int lu_gfmmat(GFMMAT mat,unsigned int md,int *perm)
1.1 noro 1859: {
1860: int row,col;
1.24 noro 1861: int i,j,k;
1.1 noro 1862: unsigned int *t,*pivot;
1863: unsigned int **a;
1864: unsigned int inv,m;
1865:
1866: row = mat->row; col = mat->col;
1867: a = mat->body;
1868: bzero(perm,row*sizeof(int));
1869:
1870: for ( i = 0; i < row; i++ )
1871: perm[i] = i;
1872: for ( k = 0; k < col; k++ ) {
1873: for ( i = k; i < row && !a[i][k]; i++ );
1874: if ( i == row )
1875: return 0;
1876: if ( i != k ) {
1877: j = perm[i]; perm[i] = perm[k]; perm[k] = j;
1878: t = a[i]; a[i] = a[k]; a[k] = t;
1879: }
1880: pivot = a[k];
1881: pivot[k] = inv = invm(pivot[k],md);
1882: for ( i = k+1; i < row; i++ ) {
1883: t = a[i];
1884: if ( m = t[k] ) {
1885: DMAR(inv,m,0,md,t[k])
1886: for ( j = k+1, m = md - t[k]; j < col; j++ )
1887: if ( pivot[j] ) {
1.8 noro 1888: unsigned int tj;
1889:
1890: DMAR(m,pivot[j],t[j],md,tj)
1891: t[j] = tj;
1.1 noro 1892: }
1893: }
1894: }
1895: }
1896: return 1;
1897: }
1898:
1.3 noro 1899: /*
1900: Input
1901: a: a row x col matrix
1902: md : a modulus
1903:
1904: Output:
1905: return : d = the rank of mat
1906: a[0..(d-1)][0..(d-1)] : LU decomposition (a[i][i] = 1/U[i][i])
1907: rinfo: array of length row
1908: cinfo: array of length col
1909: i-th row in new a <-> rinfo[i]-th row in old a
1910: cinfo[j]=1 <=> j-th column is contained in the LU decomp.
1911: */
1912:
1.24 noro 1913: int find_lhs_and_lu_mod(unsigned int **a,int row,int col,
1914: unsigned int md,int **rinfo,int **cinfo)
1.3 noro 1915: {
1.24 noro 1916: int i,j,k,d;
1.3 noro 1917: int *rp,*cp;
1918: unsigned int *t,*pivot;
1919: unsigned int inv,m;
1920:
1921: *rinfo = rp = (int *)MALLOC_ATOMIC(row*sizeof(int));
1922: *cinfo = cp = (int *)MALLOC_ATOMIC(col*sizeof(int));
1923: for ( i = 0; i < row; i++ )
1924: rp[i] = i;
1925: for ( k = 0, d = 0; k < col; k++ ) {
1926: for ( i = d; i < row && !a[i][k]; i++ );
1927: if ( i == row ) {
1928: cp[k] = 0;
1929: continue;
1930: } else
1931: cp[k] = 1;
1932: if ( i != d ) {
1933: j = rp[i]; rp[i] = rp[d]; rp[d] = j;
1934: t = a[i]; a[i] = a[d]; a[d] = t;
1935: }
1936: pivot = a[d];
1937: pivot[k] = inv = invm(pivot[k],md);
1938: for ( i = d+1; i < row; i++ ) {
1939: t = a[i];
1940: if ( m = t[k] ) {
1941: DMAR(inv,m,0,md,t[k])
1942: for ( j = k+1, m = md - t[k]; j < col; j++ )
1943: if ( pivot[j] ) {
1.8 noro 1944: unsigned int tj;
1945: DMAR(m,pivot[j],t[j],md,tj)
1946: t[j] = tj;
1.3 noro 1947: }
1948: }
1949: }
1950: d++;
1951: }
1952: return d;
1953: }
1954:
1955: /*
1956: Input
1957: a : n x n matrix; a result of LU-decomposition
1958: md : modulus
1959: b : n x l matrix
1960: Output
1961: b = a^(-1)b
1962: */
1963:
1.44 noro 1964: void solve_by_lu_mod(int **a,int n,int md,int **b,int l,int normalize)
1.3 noro 1965: {
1966: unsigned int *y,*c;
1967: int i,j,k;
1968: unsigned int t,m,m2;
1969:
1970: y = (int *)MALLOC_ATOMIC(n*sizeof(int));
1971: c = (int *)MALLOC_ATOMIC(n*sizeof(int));
1972: m2 = md>>1;
1973: for ( k = 0; k < l; k++ ) {
1974: /* copy b[.][k] to c */
1975: for ( i = 0; i < n; i++ )
1976: c[i] = (unsigned int)b[i][k];
1977: /* solve Ly=c */
1978: for ( i = 0; i < n; i++ ) {
1979: for ( t = c[i], j = 0; j < i; j++ )
1980: if ( a[i][j] ) {
1981: m = md - a[i][j];
1982: DMAR(m,y[j],t,md,t)
1983: }
1984: y[i] = t;
1985: }
1986: /* solve Uc=y */
1987: for ( i = n-1; i >= 0; i-- ) {
1988: for ( t = y[i], j =i+1; j < n; j++ )
1989: if ( a[i][j] ) {
1990: m = md - a[i][j];
1991: DMAR(m,c[j],t,md,t)
1992: }
1993: /* a[i][i] = 1/U[i][i] */
1994: DMAR(t,a[i][i],0,md,c[i])
1995: }
1996: /* copy c to b[.][k] with normalization */
1.44 noro 1997: if ( normalize )
1998: for ( i = 0; i < n; i++ )
1999: b[i][k] = (int)(c[i]>m2 ? c[i]-md : c[i]);
2000: else
2001: for ( i = 0; i < n; i++ )
2002: b[i][k] = c[i];
1.3 noro 2003: }
2004: }
2005:
1.24 noro 2006: void Pleqm1(NODE arg,VECT *rp)
1.1 noro 2007: {
2008: MAT m;
2009: VECT vect;
2010: pointer **mat;
2011: Q *v;
2012: Q q;
2013: int **wmat;
2014: int md,i,j,row,col,t,n,status;
2015:
2016: asir_assert(ARG0(arg),O_MAT,"leqm1");
2017: asir_assert(ARG1(arg),O_N,"leqm1");
2018: m = (MAT)ARG0(arg); md = QTOS((Q)ARG1(arg));
2019: row = m->row; col = m->col; mat = m->body;
2020: wmat = (int **)almat(row,col);
2021: for ( i = 0; i < row; i++ )
2022: for ( j = 0; j < col; j++ )
2023: if ( q = (Q)mat[i][j] ) {
2024: t = rem(NM(q),md);
2025: if ( SGN(q) < 0 )
2026: t = (md - t) % md;
2027: wmat[i][j] = t;
2028: } else
2029: wmat[i][j] = 0;
2030: status = gauss_elim_mod1(wmat,row,col,md);
2031: if ( status < 0 )
2032: *rp = 0;
2033: else if ( status > 0 )
2034: *rp = (VECT)ONE;
2035: else {
2036: n = col - 1;
2037: MKVECT(vect,n);
2038: for ( i = 0, v = (Q *)vect->body; i < n; i++ ) {
2039: t = (md-wmat[i][n])%md; STOQ(t,v[i]);
2040: }
2041: *rp = vect;
2042: }
2043: }
2044:
1.24 noro 2045: int gauss_elim_mod1(int **mat,int row,int col,int md)
1.1 noro 2046: {
2047: int i,j,k,inv,a,n;
2048: int *t,*pivot;
2049:
2050: n = col - 1;
2051: for ( j = 0; j < n; j++ ) {
2052: for ( i = j; i < row && !mat[i][j]; i++ );
2053: if ( i == row )
2054: return 1;
2055: if ( i != j ) {
2056: t = mat[i]; mat[i] = mat[j]; mat[j] = t;
2057: }
2058: pivot = mat[j];
2059: inv = invm(pivot[j],md);
2060: for ( k = j; k <= n; k++ )
2061: pivot[k] = dmar(pivot[k],inv,0,md);
2062: for ( i = j+1; i < row; i++ ) {
2063: t = mat[i];
2064: if ( i != j && (a = t[j]) )
2065: for ( k = j, a = md - a; k <= n; k++ )
2066: t[k] = dmar(pivot[k],a,t[k],md);
2067: }
2068: }
2069: for ( i = n; i < row && !mat[i][n]; i++ );
2070: if ( i == row ) {
2071: for ( j = n-1; j >= 0; j-- ) {
2072: for ( i = j-1, a = (md-mat[j][n])%md; i >= 0; i-- ) {
2073: mat[i][n] = dmar(mat[i][j],a,mat[i][n],md);
2074: mat[i][j] = 0;
2075: }
2076: }
2077: return 0;
2078: } else
2079: return -1;
2080: }
2081:
1.24 noro 2082: void Pgeninvm(NODE arg,LIST *rp)
1.1 noro 2083: {
2084: MAT m;
2085: pointer **mat;
2086: Q **tmat;
2087: Q q;
2088: unsigned int **wmat;
2089: int md,i,j,row,col,t,status;
2090: MAT mat1,mat2;
2091: NODE node1,node2;
2092:
2093: asir_assert(ARG0(arg),O_MAT,"leqm1");
2094: asir_assert(ARG1(arg),O_N,"leqm1");
2095: m = (MAT)ARG0(arg); md = QTOS((Q)ARG1(arg));
2096: row = m->row; col = m->col; mat = m->body;
2097: wmat = (unsigned int **)almat(row,col+row);
2098: for ( i = 0; i < row; i++ ) {
2099: bzero((char *)wmat[i],(col+row)*sizeof(int));
2100: for ( j = 0; j < col; j++ )
2101: if ( q = (Q)mat[i][j] ) {
2102: t = rem(NM(q),md);
2103: if ( SGN(q) < 0 )
2104: t = (md - t) % md;
2105: wmat[i][j] = t;
2106: }
2107: wmat[i][col+i] = 1;
2108: }
2109: status = gauss_elim_geninv_mod(wmat,row,col,md);
2110: if ( status > 0 )
2111: *rp = 0;
2112: else {
2113: MKMAT(mat1,col,row); MKMAT(mat2,row-col,row);
2114: for ( i = 0, tmat = (Q **)mat1->body; i < col; i++ )
2115: for ( j = 0; j < row; j++ )
1.24 noro 2116: UTOQ(wmat[i][j+col],tmat[i][j]);
1.1 noro 2117: for ( tmat = (Q **)mat2->body; i < row; i++ )
2118: for ( j = 0; j < row; j++ )
1.24 noro 2119: UTOQ(wmat[i][j+col],tmat[i-col][j]);
1.1 noro 2120: MKNODE(node2,mat2,0); MKNODE(node1,mat1,node2); MKLIST(*rp,node1);
2121: }
2122: }
2123:
1.24 noro 2124: int gauss_elim_geninv_mod(unsigned int **mat,int row,int col,int md)
1.1 noro 2125: {
2126: int i,j,k,inv,a,n,m;
2127: unsigned int *t,*pivot;
2128:
2129: n = col; m = row+col;
2130: for ( j = 0; j < n; j++ ) {
2131: for ( i = j; i < row && !mat[i][j]; i++ );
2132: if ( i == row )
2133: return 1;
2134: if ( i != j ) {
2135: t = mat[i]; mat[i] = mat[j]; mat[j] = t;
2136: }
2137: pivot = mat[j];
2138: inv = invm(pivot[j],md);
2139: for ( k = j; k < m; k++ )
2140: pivot[k] = dmar(pivot[k],inv,0,md);
2141: for ( i = j+1; i < row; i++ ) {
2142: t = mat[i];
2143: if ( a = t[j] )
2144: for ( k = j, a = md - a; k < m; k++ )
2145: t[k] = dmar(pivot[k],a,t[k],md);
2146: }
2147: }
2148: for ( j = n-1; j >= 0; j-- ) {
2149: pivot = mat[j];
2150: for ( i = j-1; i >= 0; i-- ) {
2151: t = mat[i];
2152: if ( a = t[j] )
2153: for ( k = j, a = md - a; k < m; k++ )
2154: t[k] = dmar(pivot[k],a,t[k],md);
2155: }
2156: }
2157: return 0;
2158: }
2159:
1.24 noro 2160: void Psolve_by_lu_gfmmat(NODE arg,VECT *rp)
1.1 noro 2161: {
2162: GFMMAT lu;
2163: Q *perm,*rhs,*v;
2164: int n,i;
2165: unsigned int md;
2166: unsigned int *b,*sol;
2167: VECT r;
2168:
2169: lu = (GFMMAT)ARG0(arg);
2170: perm = (Q *)BDY((VECT)ARG1(arg));
2171: rhs = (Q *)BDY((VECT)ARG2(arg));
2172: md = (unsigned int)QTOS((Q)ARG3(arg));
2173: n = lu->col;
2174: b = (unsigned int *)MALLOC_ATOMIC(n*sizeof(int));
2175: sol = (unsigned int *)MALLOC_ATOMIC(n*sizeof(int));
2176: for ( i = 0; i < n; i++ )
2177: b[i] = QTOS(rhs[QTOS(perm[i])]);
2178: solve_by_lu_gfmmat(lu,md,b,sol);
2179: MKVECT(r,n);
2180: for ( i = 0, v = (Q *)r->body; i < n; i++ )
1.24 noro 2181: UTOQ(sol[i],v[i]);
1.1 noro 2182: *rp = r;
2183: }
2184:
1.24 noro 2185: void solve_by_lu_gfmmat(GFMMAT lu,unsigned int md,
2186: unsigned int *b,unsigned int *x)
1.1 noro 2187: {
2188: int n;
2189: unsigned int **a;
2190: unsigned int *y;
2191: int i,j;
2192: unsigned int t,m;
2193:
2194: n = lu->col;
2195: a = lu->body;
2196: y = (unsigned int *)MALLOC_ATOMIC(n*sizeof(int));
2197: /* solve Ly=b */
2198: for ( i = 0; i < n; i++ ) {
2199: for ( t = b[i], j = 0; j < i; j++ )
2200: if ( a[i][j] ) {
2201: m = md - a[i][j];
2202: DMAR(m,y[j],t,md,t)
2203: }
2204: y[i] = t;
2205: }
2206: /* solve Ux=y */
2207: for ( i = n-1; i >= 0; i-- ) {
2208: for ( t = y[i], j =i+1; j < n; j++ )
2209: if ( a[i][j] ) {
2210: m = md - a[i][j];
2211: DMAR(m,x[j],t,md,t)
2212: }
2213: /* a[i][i] = 1/U[i][i] */
2214: DMAR(t,a[i][i],0,md,x[i])
2215: }
2216: }
2217:
1.24 noro 2218: void Plu_gfmmat(NODE arg,LIST *rp)
1.1 noro 2219: {
2220: MAT m;
2221: GFMMAT mm;
2222: unsigned int md;
2223: int i,row,col,status;
2224: int *iperm;
2225: Q *v;
2226: VECT perm;
2227: NODE n0;
2228:
2229: asir_assert(ARG0(arg),O_MAT,"mat_to_gfmmat");
2230: asir_assert(ARG1(arg),O_N,"mat_to_gfmmat");
2231: m = (MAT)ARG0(arg); md = (unsigned int)QTOS((Q)ARG1(arg));
2232: mat_to_gfmmat(m,md,&mm);
2233: row = m->row;
2234: col = m->col;
2235: iperm = (int *)MALLOC_ATOMIC(row*sizeof(int));
2236: status = lu_gfmmat(mm,md,iperm);
2237: if ( !status )
2238: n0 = 0;
2239: else {
2240: MKVECT(perm,row);
2241: for ( i = 0, v = (Q *)perm->body; i < row; i++ )
2242: STOQ(iperm[i],v[i]);
2243: n0 = mknode(2,mm,perm);
2244: }
2245: MKLIST(*rp,n0);
2246: }
2247:
1.24 noro 2248: void Pmat_to_gfmmat(NODE arg,GFMMAT *rp)
1.1 noro 2249: {
2250: MAT m;
2251: unsigned int md;
2252:
2253: asir_assert(ARG0(arg),O_MAT,"mat_to_gfmmat");
2254: asir_assert(ARG1(arg),O_N,"mat_to_gfmmat");
2255: m = (MAT)ARG0(arg); md = (unsigned int)QTOS((Q)ARG1(arg));
2256: mat_to_gfmmat(m,md,rp);
2257: }
2258:
1.24 noro 2259: void mat_to_gfmmat(MAT m,unsigned int md,GFMMAT *rp)
1.1 noro 2260: {
2261: unsigned int **wmat;
2262: unsigned int t;
2263: Q **mat;
2264: Q q;
2265: int i,j,row,col;
2266:
2267: row = m->row; col = m->col; mat = (Q **)m->body;
2268: wmat = (unsigned int **)almat(row,col);
2269: for ( i = 0; i < row; i++ ) {
2270: bzero((char *)wmat[i],col*sizeof(unsigned int));
2271: for ( j = 0; j < col; j++ )
2272: if ( q = mat[i][j] ) {
2273: t = (unsigned int)rem(NM(q),md);
2274: if ( SGN(q) < 0 )
2275: t = (md - t) % md;
2276: wmat[i][j] = t;
2277: }
2278: }
2279: TOGFMMAT(row,col,wmat,*rp);
2280: }
2281:
1.27 noro 2282: void Pgeninvm_swap(arg,rp)
2283: NODE arg;
2284: LIST *rp;
1.1 noro 2285: {
2286: MAT m;
2287: pointer **mat;
2288: Q **tmat;
2289: Q *tvect;
2290: Q q;
2291: unsigned int **wmat,**invmat;
2292: int *index;
2293: unsigned int t,md;
2294: int i,j,row,col,status;
2295: MAT mat1;
2296: VECT vect1;
2297: NODE node1,node2;
2298:
2299: asir_assert(ARG0(arg),O_MAT,"geninvm_swap");
2300: asir_assert(ARG1(arg),O_N,"geninvm_swap");
2301: m = (MAT)ARG0(arg); md = QTOS((Q)ARG1(arg));
2302: row = m->row; col = m->col; mat = m->body;
2303: wmat = (unsigned int **)almat(row,col+row);
2304: for ( i = 0; i < row; i++ ) {
2305: bzero((char *)wmat[i],(col+row)*sizeof(int));
2306: for ( j = 0; j < col; j++ )
2307: if ( q = (Q)mat[i][j] ) {
2308: t = (unsigned int)rem(NM(q),md);
2309: if ( SGN(q) < 0 )
2310: t = (md - t) % md;
2311: wmat[i][j] = t;
2312: }
2313: wmat[i][col+i] = 1;
2314: }
2315: status = gauss_elim_geninv_mod_swap(wmat,row,col,md,&invmat,&index);
2316: if ( status > 0 )
2317: *rp = 0;
2318: else {
2319: MKMAT(mat1,col,col);
2320: for ( i = 0, tmat = (Q **)mat1->body; i < col; i++ )
2321: for ( j = 0; j < col; j++ )
2322: UTOQ(invmat[i][j],tmat[i][j]);
2323: MKVECT(vect1,row);
2324: for ( i = 0, tvect = (Q *)vect1->body; i < row; i++ )
2325: STOQ(index[i],tvect[i]);
2326: MKNODE(node2,vect1,0); MKNODE(node1,mat1,node2); MKLIST(*rp,node1);
2327: }
2328: }
2329:
1.27 noro 2330: gauss_elim_geninv_mod_swap(mat,row,col,md,invmatp,indexp)
2331: unsigned int **mat;
2332: int row,col;
2333: unsigned int md;
2334: unsigned int ***invmatp;
2335: int **indexp;
1.1 noro 2336: {
2337: int i,j,k,inv,a,n,m;
2338: unsigned int *t,*pivot,*s;
2339: int *index;
2340: unsigned int **invmat;
2341:
2342: n = col; m = row+col;
2343: *indexp = index = (int *)MALLOC_ATOMIC(row*sizeof(int));
2344: for ( i = 0; i < row; i++ )
2345: index[i] = i;
2346: for ( j = 0; j < n; j++ ) {
2347: for ( i = j; i < row && !mat[i][j]; i++ );
2348: if ( i == row ) {
2349: *indexp = 0; *invmatp = 0; return 1;
2350: }
2351: if ( i != j ) {
2352: t = mat[i]; mat[i] = mat[j]; mat[j] = t;
2353: k = index[i]; index[i] = index[j]; index[j] = k;
2354: }
2355: pivot = mat[j];
2356: inv = (unsigned int)invm(pivot[j],md);
2357: for ( k = j; k < m; k++ )
2358: if ( pivot[k] )
2359: pivot[k] = (unsigned int)dmar(pivot[k],inv,0,md);
2360: for ( i = j+1; i < row; i++ ) {
2361: t = mat[i];
2362: if ( a = t[j] )
2363: for ( k = j, a = md - a; k < m; k++ )
2364: if ( pivot[k] )
2365: t[k] = dmar(pivot[k],a,t[k],md);
2366: }
2367: }
2368: for ( j = n-1; j >= 0; j-- ) {
2369: pivot = mat[j];
2370: for ( i = j-1; i >= 0; i-- ) {
2371: t = mat[i];
2372: if ( a = t[j] )
2373: for ( k = j, a = md - a; k < m; k++ )
2374: if ( pivot[k] )
2375: t[k] = dmar(pivot[k],a,t[k],md);
2376: }
2377: }
2378: *invmatp = invmat = (unsigned int **)almat(col,col);
1.27 noro 2379: for ( i = 0; i < col; i++ )
2380: for ( j = 0, s = invmat[i], t = mat[i]; j < col; j++ )
2381: s[j] = t[col+index[j]];
2382: return 0;
2383: }
2384:
2385: void Pgeninv_sf_swap(NODE arg,LIST *rp)
2386: {
2387: MAT m;
2388: GFS **mat,**tmat;
2389: Q *tvect;
2390: GFS q;
2391: int **wmat,**invmat;
2392: int *index;
2393: unsigned int t;
2394: int i,j,row,col,status;
2395: MAT mat1;
2396: VECT vect1;
2397: NODE node1,node2;
2398:
2399: asir_assert(ARG0(arg),O_MAT,"geninv_sf_swap");
2400: m = (MAT)ARG0(arg);
2401: row = m->row; col = m->col; mat = (GFS **)m->body;
2402: wmat = (int **)almat(row,col+row);
2403: for ( i = 0; i < row; i++ ) {
2404: bzero((char *)wmat[i],(col+row)*sizeof(int));
2405: for ( j = 0; j < col; j++ )
2406: if ( q = (GFS)mat[i][j] )
2407: wmat[i][j] = FTOIF(CONT(q));
2408: wmat[i][col+i] = _onesf();
2409: }
2410: status = gauss_elim_geninv_sf_swap(wmat,row,col,&invmat,&index);
2411: if ( status > 0 )
2412: *rp = 0;
2413: else {
2414: MKMAT(mat1,col,col);
2415: for ( i = 0, tmat = (GFS **)mat1->body; i < col; i++ )
2416: for ( j = 0; j < col; j++ )
2417: if ( t = invmat[i][j] ) {
2418: MKGFS(IFTOF(t),tmat[i][j]);
2419: }
2420: MKVECT(vect1,row);
2421: for ( i = 0, tvect = (Q *)vect1->body; i < row; i++ )
2422: STOQ(index[i],tvect[i]);
2423: MKNODE(node2,vect1,0); MKNODE(node1,mat1,node2); MKLIST(*rp,node1);
2424: }
2425: }
2426:
2427: int gauss_elim_geninv_sf_swap(int **mat,int row,int col,
2428: int ***invmatp,int **indexp)
2429: {
2430: int i,j,k,inv,a,n,m,u;
2431: int *t,*pivot,*s;
2432: int *index;
2433: int **invmat;
2434:
2435: n = col; m = row+col;
2436: *indexp = index = (int *)MALLOC_ATOMIC(row*sizeof(int));
2437: for ( i = 0; i < row; i++ )
2438: index[i] = i;
2439: for ( j = 0; j < n; j++ ) {
2440: for ( i = j; i < row && !mat[i][j]; i++ );
2441: if ( i == row ) {
2442: *indexp = 0; *invmatp = 0; return 1;
2443: }
2444: if ( i != j ) {
2445: t = mat[i]; mat[i] = mat[j]; mat[j] = t;
2446: k = index[i]; index[i] = index[j]; index[j] = k;
2447: }
2448: pivot = mat[j];
2449: inv = _invsf(pivot[j]);
2450: for ( k = j; k < m; k++ )
2451: if ( pivot[k] )
2452: pivot[k] = _mulsf(pivot[k],inv);
2453: for ( i = j+1; i < row; i++ ) {
2454: t = mat[i];
2455: if ( a = t[j] )
2456: for ( k = j, a = _chsgnsf(a); k < m; k++ )
2457: if ( pivot[k] ) {
2458: u = _mulsf(pivot[k],a);
2459: t[k] = _addsf(u,t[k]);
2460: }
2461: }
2462: }
2463: for ( j = n-1; j >= 0; j-- ) {
2464: pivot = mat[j];
2465: for ( i = j-1; i >= 0; i-- ) {
2466: t = mat[i];
2467: if ( a = t[j] )
2468: for ( k = j, a = _chsgnsf(a); k < m; k++ )
2469: if ( pivot[k] ) {
2470: u = _mulsf(pivot[k],a);
2471: t[k] = _addsf(u,t[k]);
2472: }
2473: }
2474: }
2475: *invmatp = invmat = (int **)almat(col,col);
1.1 noro 2476: for ( i = 0; i < col; i++ )
2477: for ( j = 0, s = invmat[i], t = mat[i]; j < col; j++ )
2478: s[j] = t[col+index[j]];
2479: return 0;
2480: }
2481:
2482: void _addn(N,N,N);
2483: int _subn(N,N,N);
2484: void _muln(N,N,N);
2485:
1.24 noro 2486: void inner_product_int(Q *a,Q *b,int n,Q *r)
1.1 noro 2487: {
2488: int la,lb,i;
2489: int sgn,sgn1;
2490: N wm,wma,sum,t;
2491:
2492: for ( la = lb = 0, i = 0; i < n; i++ ) {
2493: if ( a[i] )
2494: if ( DN(a[i]) )
2495: error("inner_product_int : invalid argument");
2496: else
2497: la = MAX(PL(NM(a[i])),la);
2498: if ( b[i] )
2499: if ( DN(b[i]) )
2500: error("inner_product_int : invalid argument");
2501: else
2502: lb = MAX(PL(NM(b[i])),lb);
2503: }
2504: sgn = 0;
2505: sum= NALLOC(la+lb+2);
2506: bzero((char *)sum,(la+lb+3)*sizeof(unsigned int));
2507: wm = NALLOC(la+lb+2);
2508: wma = NALLOC(la+lb+2);
2509: for ( i = 0; i < n; i++ ) {
2510: if ( !a[i] || !b[i] )
2511: continue;
2512: _muln(NM(a[i]),NM(b[i]),wm);
2513: sgn1 = SGN(a[i])*SGN(b[i]);
2514: if ( !sgn ) {
2515: sgn = sgn1;
2516: t = wm; wm = sum; sum = t;
2517: } else if ( sgn == sgn1 ) {
2518: _addn(sum,wm,wma);
2519: if ( !PL(wma) )
2520: sgn = 0;
2521: t = wma; wma = sum; sum = t;
2522: } else {
2523: /* sgn*sum+sgn1*wm = sgn*(sum-wm) */
2524: sgn *= _subn(sum,wm,wma);
2525: t = wma; wma = sum; sum = t;
2526: }
2527: }
2528: GC_free(wm);
2529: GC_free(wma);
2530: if ( !sgn ) {
2531: GC_free(sum);
2532: *r = 0;
2533: } else
2534: NTOQ(sum,sgn,*r);
2535: }
2536:
1.3 noro 2537: /* (k,l) element of a*b where a: .x n matrix, b: n x . integer matrix */
2538:
1.24 noro 2539: void inner_product_mat_int_mod(Q **a,int **b,int n,int k,int l,Q *r)
1.3 noro 2540: {
2541: int la,lb,i;
2542: int sgn,sgn1;
2543: N wm,wma,sum,t;
2544: Q aki;
2545: int bil,bilsgn;
2546: struct oN tn;
2547:
2548: for ( la = 0, i = 0; i < n; i++ ) {
2549: if ( aki = a[k][i] )
2550: if ( DN(aki) )
2551: error("inner_product_int : invalid argument");
2552: else
2553: la = MAX(PL(NM(aki)),la);
2554: }
2555: lb = 1;
2556: sgn = 0;
2557: sum= NALLOC(la+lb+2);
2558: bzero((char *)sum,(la+lb+3)*sizeof(unsigned int));
2559: wm = NALLOC(la+lb+2);
2560: wma = NALLOC(la+lb+2);
2561: for ( i = 0; i < n; i++ ) {
2562: if ( !(aki = a[k][i]) || !(bil = b[i][l]) )
2563: continue;
2564: tn.p = 1;
2565: if ( bil > 0 ) {
2566: tn.b[0] = bil; bilsgn = 1;
2567: } else {
2568: tn.b[0] = -bil; bilsgn = -1;
2569: }
2570: _muln(NM(aki),&tn,wm);
2571: sgn1 = SGN(aki)*bilsgn;
2572: if ( !sgn ) {
2573: sgn = sgn1;
2574: t = wm; wm = sum; sum = t;
2575: } else if ( sgn == sgn1 ) {
2576: _addn(sum,wm,wma);
2577: if ( !PL(wma) )
2578: sgn = 0;
2579: t = wma; wma = sum; sum = t;
2580: } else {
2581: /* sgn*sum+sgn1*wm = sgn*(sum-wm) */
2582: sgn *= _subn(sum,wm,wma);
2583: t = wma; wma = sum; sum = t;
2584: }
2585: }
2586: GC_free(wm);
2587: GC_free(wma);
2588: if ( !sgn ) {
2589: GC_free(sum);
2590: *r = 0;
2591: } else
2592: NTOQ(sum,sgn,*r);
2593: }
2594:
1.24 noro 2595: void Pmul_mat_vect_int(NODE arg,VECT *rp)
1.1 noro 2596: {
2597: MAT mat;
2598: VECT vect,r;
2599: int row,col,i;
2600:
2601: mat = (MAT)ARG0(arg);
2602: vect = (VECT)ARG1(arg);
2603: row = mat->row;
2604: col = mat->col;
2605: MKVECT(r,row);
1.24 noro 2606: for ( i = 0; i < row; i++ ) {
2607: inner_product_int((Q *)mat->body[i],(Q *)vect->body,col,(Q *)&r->body[i]);
2608: }
1.1 noro 2609: *rp = r;
2610: }
2611:
1.24 noro 2612: void Pnbpoly_up2(NODE arg,GF2N *rp)
1.1 noro 2613: {
2614: int m,type,ret;
2615: UP2 r;
2616:
2617: m = QTOS((Q)ARG0(arg));
2618: type = QTOS((Q)ARG1(arg));
2619: ret = generate_ONB_polynomial(&r,m,type);
2620: if ( ret == 0 )
2621: MKGF2N(r,*rp);
2622: else
2623: *rp = 0;
2624: }
2625:
1.24 noro 2626: void Px962_irredpoly_up2(NODE arg,GF2N *rp)
1.1 noro 2627: {
1.24 noro 2628: int m,ret,w;
1.1 noro 2629: GF2N prev;
2630: UP2 r;
2631:
2632: m = QTOS((Q)ARG0(arg));
2633: prev = (GF2N)ARG1(arg);
2634: if ( !prev ) {
2635: w = (m>>5)+1; NEWUP2(r,w); r->w = 0;
2636: bzero((char *)r->b,w*sizeof(unsigned int));
2637: } else {
2638: r = prev->body;
2639: if ( degup2(r) != m ) {
2640: w = (m>>5)+1; NEWUP2(r,w); r->w = 0;
2641: bzero((char *)r->b,w*sizeof(unsigned int));
2642: }
2643: }
1.24 noro 2644: ret = _generate_irreducible_polynomial(r,m);
1.1 noro 2645: if ( ret == 0 )
2646: MKGF2N(r,*rp);
2647: else
2648: *rp = 0;
2649: }
2650:
1.24 noro 2651: void Pirredpoly_up2(NODE arg,GF2N *rp)
1.1 noro 2652: {
1.24 noro 2653: int m,ret,w;
1.1 noro 2654: GF2N prev;
2655: UP2 r;
2656:
2657: m = QTOS((Q)ARG0(arg));
2658: prev = (GF2N)ARG1(arg);
2659: if ( !prev ) {
2660: w = (m>>5)+1; NEWUP2(r,w); r->w = 0;
2661: bzero((char *)r->b,w*sizeof(unsigned int));
2662: } else {
2663: r = prev->body;
2664: if ( degup2(r) != m ) {
2665: w = (m>>5)+1; NEWUP2(r,w); r->w = 0;
2666: bzero((char *)r->b,w*sizeof(unsigned int));
2667: }
2668: }
1.24 noro 2669: ret = _generate_good_irreducible_polynomial(r,m);
1.1 noro 2670: if ( ret == 0 )
2671: MKGF2N(r,*rp);
2672: else
2673: *rp = 0;
2674: }
2675:
1.26 noro 2676: void Pmat_swap_row_destructive(NODE arg, MAT *m)
2677: {
2678: int i1,i2;
2679: pointer *t;
2680: MAT mat;
2681:
2682: asir_assert(ARG0(arg),O_MAT,"mat_swap_row_destructive");
2683: asir_assert(ARG1(arg),O_N,"mat_swap_row_destructive");
2684: asir_assert(ARG2(arg),O_N,"mat_swap_row_destructive");
2685: mat = (MAT)ARG0(arg);
2686: i1 = QTOS((Q)ARG1(arg));
2687: i2 = QTOS((Q)ARG2(arg));
2688: if ( i1 < 0 || i2 < 0 || i1 >= mat->row || i2 >= mat->row )
2689: error("mat_swap_row_destructive : Out of range");
2690: t = mat->body[i1];
2691: mat->body[i1] = mat->body[i2];
2692: mat->body[i2] = t;
2693: *m = mat;
2694: }
2695:
2696: void Pmat_swap_col_destructive(NODE arg, MAT *m)
2697: {
2698: int j1,j2,i,n;
2699: pointer *mi;
2700: pointer t;
2701: MAT mat;
2702:
2703: asir_assert(ARG0(arg),O_MAT,"mat_swap_col_destructive");
2704: asir_assert(ARG1(arg),O_N,"mat_swap_col_destructive");
2705: asir_assert(ARG2(arg),O_N,"mat_swap_col_destructive");
2706: mat = (MAT)ARG0(arg);
2707: j1 = QTOS((Q)ARG1(arg));
2708: j2 = QTOS((Q)ARG2(arg));
2709: if ( j1 < 0 || j2 < 0 || j1 >= mat->col || j2 >= mat->col )
2710: error("mat_swap_col_destructive : Out of range");
2711: n = mat->row;
2712: for ( i = 0; i < n; i++ ) {
2713: mi = mat->body[i];
2714: t = mi[j1]; mi[j1] = mi[j2]; mi[j2] = t;
2715: }
2716: *m = mat;
2717: }
1.1 noro 2718: /*
2719: * f = type 'type' normal polynomial of degree m if exists
2720: * IEEE P1363 A.7.2
2721: *
2722: * return value : 0 --- exists
2723: * 1 --- does not exist
2724: * -1 --- failure (memory allocation error)
2725: */
2726:
2727: int generate_ONB_polynomial(UP2 *rp,int m,int type)
2728: {
2729: int i,r;
2730: int w;
2731: UP2 f,f0,f1,f2,t;
2732:
2733: w = (m>>5)+1;
2734: switch ( type ) {
2735: case 1:
2736: if ( !TypeT_NB_check(m,1) ) return 1;
2737: NEWUP2(f,w); *rp = f; f->w = w;
2738: /* set all the bits */
2739: for ( i = 0; i < w; i++ )
2740: f->b[i] = 0xffffffff;
2741: /* mask the top word if necessary */
2742: if ( r = (m+1)&31 )
2743: f->b[w-1] &= (1<<r)-1;
2744: return 0;
2745: break;
2746: case 2:
2747: if ( !TypeT_NB_check(m,2) ) return 1;
2748: NEWUP2(f,w); *rp = f;
2749: W_NEWUP2(f0,w);
2750: W_NEWUP2(f1,w);
2751: W_NEWUP2(f2,w);
2752:
2753: /* recursion for genrating Type II normal polynomial */
2754:
2755: /* f0 = 1, f1 = t+1 */
2756: f0->w = 1; f0->b[0] = 1;
2757: f1->w = 1; f1->b[0] = 3;
2758: for ( i = 2; i <= m; i++ ) {
2759: /* f2 = t*f1+f0 */
2760: _bshiftup2(f1,-1,f2);
2761: _addup2_destructive(f2,f0);
2762: /* cyclic change of the variables */
2763: t = f0; f0 = f1; f1 = f2; f2 = t;
2764: }
2765: _copyup2(f1,f);
2766: return 0;
2767: break;
2768: default:
2769: return -1;
2770: break;
2771: }
2772: }
2773:
2774: /*
2775: * f = an irreducible trinomial or pentanomial of degree d 'after' f
2776: * return value : 0 --- exists
2777: * 1 --- does not exist (exhaustion)
2778: */
2779:
2780: int _generate_irreducible_polynomial(UP2 f,int d)
2781: {
2782: int ret,i,j,k,nz,i0,j0,k0;
2783: int w;
2784: unsigned int *fd;
2785:
2786: /*
2787: * if f = x^d+x^i+1 then i0 <- i, j0 <- 0, k0 <-0.
2788: * if f = x^d+x^k+x^j+x^i+1 (k>j>i) then i0 <- i, j0 <- j, k0 <-k.
2789: * otherwise i0,j0,k0 is set to 0.
2790: */
2791:
2792: fd = f->b;
2793: w = (d>>5)+1;
2794: if ( f->w && (d==degup2(f)) ) {
2795: for ( nz = 0, i = d; i >= 0; i-- )
2796: if ( fd[i>>5]&(1<<(i&31)) ) nz++;
2797: switch ( nz ) {
2798: case 3:
2799: for ( i0 = 1; !(fd[i0>>5]&(1<<(i0&31))) ; i0++ );
2800: /* reset i0-th bit */
2801: fd[i0>>5] &= ~(1<<(i0&31));
2802: j0 = k0 = 0;
2803: break;
2804: case 5:
2805: for ( i0 = 1; !(fd[i0>>5]&(1<<(i0&31))) ; i0++ );
2806: /* reset i0-th bit */
2807: fd[i0>>5] &= ~(1<<(i0&31));
2808: for ( j0 = i0+1; !(fd[j0>>5]&(1<<(j0&31))) ; j0++ );
2809: /* reset j0-th bit */
2810: fd[j0>>5] &= ~(1<<(j0&31));
2811: for ( k0 = j0+1; !(fd[k0>>5]&(1<<(k0&31))) ; k0++ );
2812: /* reset k0-th bit */
2813: fd[k0>>5] &= ~(1<<(k0&31));
2814: break;
2815: default:
2816: f->w = 0; break;
2817: }
2818: } else
2819: f->w = 0;
2820:
2821: if ( !f->w ) {
2822: fd = f->b;
2823: f->w = w; fd[0] |= 1; fd[d>>5] |= (1<<(d&31));
2824: i0 = j0 = k0 = 0;
2825: }
2826: /* if j0 > 0 then f is already a pentanomial */
2827: if ( j0 > 0 ) goto PENTA;
2828:
2829: /* searching for an irreducible trinomial */
2830:
2831: for ( i = 1; 2*i <= d; i++ ) {
2832: /* skip the polynomials 'before' f */
2833: if ( i < i0 ) continue;
2834: if ( i == i0 ) { i0 = 0; continue; }
2835: /* set i-th bit */
2836: fd[i>>5] |= (1<<(i&31));
2837: ret = irredcheck_dddup2(f);
2838: if ( ret == 1 ) return 0;
2839: /* reset i-th bit */
2840: fd[i>>5] &= ~(1<<(i&31));
2841: }
2842:
2843: /* searching for an irreducible pentanomial */
2844: PENTA:
2845: for ( i = 1; i < d; i++ ) {
2846: /* skip the polynomials 'before' f */
2847: if ( i < i0 ) continue;
2848: if ( i == i0 ) i0 = 0;
2849: /* set i-th bit */
2850: fd[i>>5] |= (1<<(i&31));
2851: for ( j = i+1; j < d; j++ ) {
2852: /* skip the polynomials 'before' f */
2853: if ( j < j0 ) continue;
2854: if ( j == j0 ) j0 = 0;
2855: /* set j-th bit */
2856: fd[j>>5] |= (1<<(j&31));
2857: for ( k = j+1; k < d; k++ ) {
2858: /* skip the polynomials 'before' f */
2859: if ( k < k0 ) continue;
2860: else if ( k == k0 ) { k0 = 0; continue; }
2861: /* set k-th bit */
2862: fd[k>>5] |= (1<<(k&31));
2863: ret = irredcheck_dddup2(f);
2864: if ( ret == 1 ) return 0;
2865: /* reset k-th bit */
2866: fd[k>>5] &= ~(1<<(k&31));
2867: }
2868: /* reset j-th bit */
2869: fd[j>>5] &= ~(1<<(j&31));
2870: }
2871: /* reset i-th bit */
2872: fd[i>>5] &= ~(1<<(i&31));
2873: }
2874: /* exhausted */
2875: return 1;
2876: }
2877:
2878: /*
2879: * f = an irreducible trinomial or pentanomial of degree d 'after' f
2880: *
2881: * searching strategy:
2882: * trinomial x^d+x^i+1:
2883: * i is as small as possible.
2884: * trinomial x^d+x^i+x^j+x^k+1:
2885: * i is as small as possible.
2886: * For such i, j is as small as possible.
2887: * For such i and j, 'k' is as small as possible.
2888: *
2889: * return value : 0 --- exists
2890: * 1 --- does not exist (exhaustion)
2891: */
2892:
2893: int _generate_good_irreducible_polynomial(UP2 f,int d)
2894: {
2895: int ret,i,j,k,nz,i0,j0,k0;
2896: int w;
2897: unsigned int *fd;
2898:
2899: /*
2900: * if f = x^d+x^i+1 then i0 <- i, j0 <- 0, k0 <-0.
2901: * if f = x^d+x^k+x^j+x^i+1 (k>j>i) then i0 <- i, j0 <- j, k0 <-k.
2902: * otherwise i0,j0,k0 is set to 0.
2903: */
2904:
2905: fd = f->b;
2906: w = (d>>5)+1;
2907: if ( f->w && (d==degup2(f)) ) {
2908: for ( nz = 0, i = d; i >= 0; i-- )
2909: if ( fd[i>>5]&(1<<(i&31)) ) nz++;
2910: switch ( nz ) {
2911: case 3:
2912: for ( i0 = 1; !(fd[i0>>5]&(1<<(i0&31))) ; i0++ );
2913: /* reset i0-th bit */
2914: fd[i0>>5] &= ~(1<<(i0&31));
2915: j0 = k0 = 0;
2916: break;
2917: case 5:
2918: for ( i0 = 1; !(fd[i0>>5]&(1<<(i0&31))) ; i0++ );
2919: /* reset i0-th bit */
2920: fd[i0>>5] &= ~(1<<(i0&31));
2921: for ( j0 = i0+1; !(fd[j0>>5]&(1<<(j0&31))) ; j0++ );
2922: /* reset j0-th bit */
2923: fd[j0>>5] &= ~(1<<(j0&31));
2924: for ( k0 = j0+1; !(fd[k0>>5]&(1<<(k0&31))) ; k0++ );
2925: /* reset k0-th bit */
2926: fd[k0>>5] &= ~(1<<(k0&31));
2927: break;
2928: default:
2929: f->w = 0; break;
2930: }
2931: } else
2932: f->w = 0;
2933:
2934: if ( !f->w ) {
2935: fd = f->b;
2936: f->w = w; fd[0] |= 1; fd[d>>5] |= (1<<(d&31));
2937: i0 = j0 = k0 = 0;
2938: }
2939: /* if j0 > 0 then f is already a pentanomial */
2940: if ( j0 > 0 ) goto PENTA;
2941:
2942: /* searching for an irreducible trinomial */
2943:
2944: for ( i = 1; 2*i <= d; i++ ) {
2945: /* skip the polynomials 'before' f */
2946: if ( i < i0 ) continue;
2947: if ( i == i0 ) { i0 = 0; continue; }
2948: /* set i-th bit */
2949: fd[i>>5] |= (1<<(i&31));
2950: ret = irredcheck_dddup2(f);
2951: if ( ret == 1 ) return 0;
2952: /* reset i-th bit */
2953: fd[i>>5] &= ~(1<<(i&31));
2954: }
2955:
2956: /* searching for an irreducible pentanomial */
2957: PENTA:
2958: for ( i = 3; i < d; i++ ) {
2959: /* skip the polynomials 'before' f */
2960: if ( i < i0 ) continue;
2961: if ( i == i0 ) i0 = 0;
2962: /* set i-th bit */
2963: fd[i>>5] |= (1<<(i&31));
2964: for ( j = 2; j < i; j++ ) {
2965: /* skip the polynomials 'before' f */
2966: if ( j < j0 ) continue;
2967: if ( j == j0 ) j0 = 0;
2968: /* set j-th bit */
2969: fd[j>>5] |= (1<<(j&31));
2970: for ( k = 1; k < j; k++ ) {
2971: /* skip the polynomials 'before' f */
2972: if ( k < k0 ) continue;
2973: else if ( k == k0 ) { k0 = 0; continue; }
2974: /* set k-th bit */
2975: fd[k>>5] |= (1<<(k&31));
2976: ret = irredcheck_dddup2(f);
2977: if ( ret == 1 ) return 0;
2978: /* reset k-th bit */
2979: fd[k>>5] &= ~(1<<(k&31));
2980: }
2981: /* reset j-th bit */
2982: fd[j>>5] &= ~(1<<(j&31));
2983: }
2984: /* reset i-th bit */
2985: fd[i>>5] &= ~(1<<(i&31));
2986: }
2987: /* exhausted */
2988: return 1;
1.3 noro 2989: }
2990:
1.24 noro 2991: void printqmat(Q **mat,int row,int col)
1.3 noro 2992: {
2993: int i,j;
2994:
2995: for ( i = 0; i < row; i++ ) {
2996: for ( j = 0; j < col; j++ ) {
1.8 noro 2997: printnum((Num)mat[i][j]); printf(" ");
1.3 noro 2998: }
2999: printf("\n");
3000: }
3001: }
3002:
1.24 noro 3003: void printimat(int **mat,int row,int col)
1.3 noro 3004: {
3005: int i,j;
3006:
3007: for ( i = 0; i < row; i++ ) {
3008: for ( j = 0; j < col; j++ ) {
3009: printf("%d ",mat[i][j]);
3010: }
3011: printf("\n");
3012: }
1.36 noro 3013: }
3014:
3015: void Pnd_det(NODE arg,P *rp)
3016: {
1.37 noro 3017: if ( argc(arg) == 1 )
3018: nd_det(0,ARG0(arg),rp);
3019: else
3020: nd_det(QTOS((Q)ARG1(arg)),ARG0(arg),rp);
1.1 noro 3021: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>