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