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