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