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