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