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