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