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