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