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