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