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