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