![[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.71 ! noro 48: * $OpenXM: OpenXM_contrib2/asir2000/builtin/array.c,v 1.70 2017/01/08 03:05:39 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:
1.51 noro 55: #include <sys/types.h>
56: #include <sys/stat.h>
1.58 ohara 57: #if !defined(_MSC_VER)
1.51 noro 58: #include <unistd.h>
1.58 ohara 59: #endif
1.51 noro 60:
1.38 noro 61: #define F4_INTRAT_PERIOD 8
62:
1.4 noro 63: #if 0
1.1 noro 64: #undef DMAR
65: #define DMAR(a1,a2,a3,d,r) (r)=dmar(a1,a2,a3,d);
1.4 noro 66: #endif
1.1 noro 67:
1.11 noro 68: extern int DP_Print; /* XXX */
1.1 noro 69:
1.24 noro 70:
1.71 ! noro 71: void Pnewvect(), Pnewmat(), Psepvect(), Psize(), Pdet(), Pleqm(), Pleqm1(), Pgeninvm(), Ptriangleq();
1.23 noro 72: void Pinvmat();
1.49 noro 73: void Pnewbytearray(),Pmemoryplot_to_coord();
1.1 noro 74:
1.25 noro 75: void Pgeneric_gauss_elim();
1.1 noro 76: void Pgeneric_gauss_elim_mod();
77:
1.69 noro 78: void Pindep_rows_mod();
79:
1.1 noro 80: void Pmat_to_gfmmat(),Plu_gfmmat(),Psolve_by_lu_gfmmat();
1.33 noro 81: void Pgeninvm_swap(), Premainder(), Psremainder(), Pvtol(), Pltov();
1.27 noro 82: void Pgeninv_sf_swap();
1.1 noro 83: void sepvect();
84: void Pmulmat_gf2n();
85: void Pbconvmat_gf2n();
86: void Pmul_vect_mat_gf2n();
87: void PNBmul_gf2n();
88: void Pmul_mat_vect_int();
89: void Psepmat_destructive();
90: void Px962_irredpoly_up2();
91: void Pirredpoly_up2();
92: void Pnbpoly_up2();
93: void Pqsort();
1.14 noro 94: void Pexponent_vector();
1.26 noro 95: void Pmat_swap_row_destructive();
96: void Pmat_swap_col_destructive();
1.28 saito 97: void Pvect();
98: void Pmat();
1.29 saito 99: void Pmatc();
1.36 noro 100: void Pnd_det();
1.53 noro 101: void Plu_mat();
1.59 ohara 102: void Pmat_col();
1.63 noro 103: void Plusolve_prep();
104: void Plusolve_main();
1.1 noro 105:
106: struct ftab array_tab[] = {
1.53 noro 107: {"lu_mat",Plu_mat,1},
1.1 noro 108: {"solve_by_lu_gfmmat",Psolve_by_lu_gfmmat,4},
109: {"lu_gfmmat",Plu_gfmmat,2},
110: {"mat_to_gfmmat",Pmat_to_gfmmat,2},
1.25 noro 111: {"generic_gauss_elim",Pgeneric_gauss_elim,1},
1.1 noro 112: {"generic_gauss_elim_mod",Pgeneric_gauss_elim_mod,2},
1.69 noro 113: {"indep_rows_mod",Pindep_rows_mod,2},
1.1 noro 114: {"newvect",Pnewvect,-2},
1.28 saito 115: {"vect",Pvect,-99999999},
1.14 noro 116: {"vector",Pnewvect,-2},
117: {"exponent_vector",Pexponent_vector,-99999999},
1.1 noro 118: {"newmat",Pnewmat,-3},
1.14 noro 119: {"matrix",Pnewmat,-3},
1.28 saito 120: {"mat",Pmat,-99999999},
1.29 saito 121: {"matr",Pmat,-99999999},
122: {"matc",Pmatc,-99999999},
1.9 noro 123: {"newbytearray",Pnewbytearray,-2},
1.49 noro 124: {"memoryplot_to_coord",Pmemoryplot_to_coord,1},
1.1 noro 125: {"sepmat_destructive",Psepmat_destructive,2},
126: {"sepvect",Psepvect,2},
127: {"qsort",Pqsort,-2},
128: {"vtol",Pvtol,1},
1.33 noro 129: {"ltov",Pltov,1},
1.1 noro 130: {"size",Psize,1},
131: {"det",Pdet,-2},
1.37 noro 132: {"nd_det",Pnd_det,-2},
1.23 noro 133: {"invmat",Pinvmat,-2},
1.1 noro 134: {"leqm",Pleqm,2},
135: {"leqm1",Pleqm1,2},
136: {"geninvm",Pgeninvm,2},
137: {"geninvm_swap",Pgeninvm_swap,2},
1.27 noro 138: {"geninv_sf_swap",Pgeninv_sf_swap,1},
1.1 noro 139: {"remainder",Premainder,2},
140: {"sremainder",Psremainder,2},
141: {"mulmat_gf2n",Pmulmat_gf2n,1},
142: {"bconvmat_gf2n",Pbconvmat_gf2n,-4},
143: {"mul_vect_mat_gf2n",Pmul_vect_mat_gf2n,2},
144: {"mul_mat_vect_int",Pmul_mat_vect_int,2},
145: {"nbmul_gf2n",PNBmul_gf2n,3},
146: {"x962_irredpoly_up2",Px962_irredpoly_up2,2},
147: {"irredpoly_up2",Pirredpoly_up2,2},
148: {"nbpoly_up2",Pnbpoly_up2,2},
1.26 noro 149: {"mat_swap_row_destructive",Pmat_swap_row_destructive,3},
150: {"mat_swap_col_destructive",Pmat_swap_col_destructive,3},
1.59 ohara 151: {"mat_col",Pmat_col,2},
1.63 noro 152: {"lusolve_prep",Plusolve_prep,1},
153: {"lusolve_main",Plusolve_main,1},
1.71 ! noro 154: {"triangleq",Ptriangleq,1},
1.1 noro 155: {0,0,0},
156: };
157:
1.63 noro 158: typedef struct _ent { int j; unsigned int e; } ent;
159:
160: ent *get_row(FILE *,int *l);
161: void put_row(FILE *out,int l,ent *a);
162: int lu_elim(int *l,ent **a,int k,int i,int mul,int mod);
163:
164: static int *ul,*ll;
165: static ent **u,**l;
166: static int modulus;
167:
168: void Plusolve_prep(NODE arg,Q *rp)
169: {
170: char *fname;
171: FILE *in;
172: int len,i,rank;
173: int *rhs;
174:
175: fname = BDY((STRING)ARG0(arg));
176: in = fopen(fname,"r");
177: modulus = getw(in);
178: len = getw(in);
179: ul = (int *)MALLOC_ATOMIC(len*sizeof(int));
180: u = (ent **)MALLOC(len*sizeof(ent *));
181: ll = (int *)MALLOC_ATOMIC(len*sizeof(int));
182: l = (ent **)MALLOC(len*sizeof(ent *));
183: for ( i = 0; i < len; i++ ) {
184: u[i] = get_row(in,&ul[i]);
185: }
186: for ( i = 0; i < len; i++ ) {
187: l[i] = get_row(in,&ll[i]);
188: }
189: fclose(in);
190: *rp = ONE;
191: }
192:
193: void Plusolve_main(NODE arg,VECT *rp)
194: {
195: Q *d,*p;
196: VECT v,r;
197: int len,i;
198: int *rhs;
199:
200: v = (VECT)ARG0(arg); len = v->len;
201: d = (Q *)BDY(v);
202: rhs = (int *)MALLOC_ATOMIC(len*sizeof(int));
203: for ( i = 0; i < len; i++ ) rhs[i] = QTOS(d[i]);
204: solve_l(ll,l,len,rhs,modulus);
205: solve_u(ul,u,len,rhs,modulus);
206: NEWVECT(r); r->len = len;
207: r->body = (pointer *)MALLOC(len*sizeof(pointer));
208: p = (Q *)r->body;
209: for ( i = 0; i < len; i++ )
210: STOQ(rhs[i],p[i]);
211: *rp = r;
212: }
213:
214: ent *get_row(FILE *in,int *l)
215: {
216: int len,i;
217: ent *a;
218:
219: *l = len = getw(in);
220: a = (ent *)MALLOC_ATOMIC(len*sizeof(ent));
221: for ( i = 0; i < len; i++ ) {
222: a[i].j = getw(in);
223: a[i].e = getw(in);
224: }
225: return a;
226: }
227:
228: int lu_gauss(int *ul,ent **u,int *ll,ent **l,int n,int mod)
229: {
230: int i,j,k,s,mul;
231: unsigned int inv;
232: int *ll2;
233:
234: ll2 = (int *)MALLOC_ATOMIC(n*sizeof(int));
235: for ( i = 0; i < n; i++ ) ll2[i] = 0;
236: for ( i = 0; i < n; i++ ) {
237: fprintf(stderr,"i=%d\n",i);
238: inv = invm(u[i][0].e,mod);
239: for ( k = i+1; k < n; k++ )
240: if ( u[k][0].j == n-i ) {
241: s = u[k][0].e;
242: DMAR(s,inv,0,mod,mul);
243: lu_elim(ul,u,k,i,mul,mod);
244: lu_append(ll,l,ll2,k,i,mul);
245: }
246: }
247: }
248:
249: #define INITLEN 10
250:
251: lu_append(int *l,ent **a,int *l2,int k,int i,int mul)
252: {
253: int len;
254: ent *p;
255:
256: len = l[k];
257: if ( !len ) {
258: a[k] = p = (ent *)MALLOC_ATOMIC(INITLEN*sizeof(ent));
259: p[0].j = i; p[0].e = mul;
260: l[k] = 1; l2[k] = INITLEN;
261: } else {
262: if ( l2[k] == l[k] ) {
263: l2[k] *= 2;
264: a[k] = REALLOC(a[k],l2[k]*sizeof(ent));
265: }
266: p =a[k];
267: p[l[k]].j = i; p[l[k]].e = mul;
268: l[k]++;
269: }
270: }
271:
272: /* a[k] = a[k]-mul*a[i] */
273:
274: int lu_elim(int *l,ent **a,int k,int i,int mul,int mod)
275: {
276: ent *ak,*ai,*w;
277: int lk,li,j,m,p,q,r,s,t,j0;
278:
279: ak = a[k]; ai = a[i]; lk = l[k]; li = l[i];
280: w = (ent *)alloca((lk+li)*sizeof(ent));
281: p = 0; q = 0; j = 0;
282: mul = mod-mul;
283: while ( p < lk && q < li ) {
284: if ( ak[p].j > ai[q].j ) {
285: w[j] = ak[p]; j++; p++;
286: } else if ( ak[p].j < ai[q].j ) {
287: w[j].j = ai[q].j;
288: t = ai[q].e;
289: DMAR(t,mul,0,mod,r);
290: w[j].e = r;
291: j++; q++;
292: } else {
293: t = ai[q].e; s = ak[p].e;
294: DMAR(t,mul,s,mod,r);
295: if ( r ) {
296: w[j].j = ai[q].j; w[j].e = r; j++;
297: }
298: p++; q++;
299: }
300: }
301: if ( q == li )
302: while ( p < lk ) {
303: w[j] = ak[p]; j++; p++;
304: }
305: else if ( p == lk )
306: while ( q < li ) {
307: w[j].j = ai[q].j;
308: t = ai[q].e;
309: DMAR(t,mul,0,mod,r);
310: w[j].e = r;
311: j++; q++;
312: }
313: if ( j <= lk ) {
314: for ( m = 0; m < j; m++ ) ak[m] = w[m];
315: } else {
316: a[k] = ak = (ent *)MALLOC_ATOMIC(j*sizeof(ent));
317: for ( m = 0; m < j; m++ ) ak[m] = w[m];
318: }
319: l[k] = j;
320: }
321:
322: int solve_l(int *ll,ent **l,int n,int *rhs,int mod)
323: {
324: int j,k,s,len;
325: ent *p;
326:
327: for ( j = 0; j < n; j++ ) {
328: len = ll[j]; p = l[j];
329: for ( k = 0, s = 0; k < len; k++ )
330: s = dmar(p[k].e,rhs[p[k].j],s,mod);
331: rhs[j] -= s;
332: if ( rhs[j] < 0 ) rhs[j] += mod;
333: }
334: }
335:
336: int solve_u(int *ul,ent **u,int n,int *rhs,int mod)
337: {
338: int j,k,s,len,inv;
339: ent *p;
340:
341: for ( j = n-1; j >= 0; j-- ) {
342: len = ul[j]; p = u[j];
343: for ( k = 1, s = 0; k < len; k++ )
344: s = dmar(p[k].e,rhs[p[k].j],s,mod);
345: rhs[j] -= s;
346: if ( rhs[j] < 0 ) rhs[j] += mod;
347: inv = invm((unsigned int)p[0].e,mod);
348: rhs[j] = dmar(rhs[j],inv,0,mod);
349: }
350: }
351:
1.24 noro 352: int comp_obj(Obj *a,Obj *b)
1.1 noro 353: {
354: return arf_comp(CO,*a,*b);
355: }
356:
357: static FUNC generic_comp_obj_func;
358: static NODE generic_comp_obj_arg;
1.60 ohara 359: static NODE generic_comp_obj_option;
1.1 noro 360:
1.24 noro 361: int generic_comp_obj(Obj *a,Obj *b)
1.1 noro 362: {
363: Q r;
364:
365: BDY(generic_comp_obj_arg)=(pointer)(*a);
366: BDY(NEXT(generic_comp_obj_arg))=(pointer)(*b);
1.60 ohara 367: r = (Q)bevalf_with_opts(generic_comp_obj_func,generic_comp_obj_arg,generic_comp_obj_option);
1.1 noro 368: if ( !r )
369: return 0;
370: else
371: return SGN(r)>0?1:-1;
372: }
373:
374:
1.46 saito 375: void Pqsort(NODE arg,LIST *rp)
1.1 noro 376: {
377: VECT vect;
1.35 ohara 378: NODE n,n1;
1.1 noro 379: P p;
380: V v;
1.34 noro 381: FUNC func;
1.35 ohara 382: int len,i;
383: pointer *a;
384: Obj t;
385:
386: t = ARG0(arg);
387: if (OID(t) == O_LIST) {
388: n = (NODE)BDY((LIST)t);
389: len = length(n);
390: MKVECT(vect,len);
391: for ( i = 0; i < len; i++, n = NEXT(n) ) {
392: BDY(vect)[i] = BDY(n);
393: }
394:
395: }else if (OID(t) != O_VECT) {
396: error("qsort : invalid argument");
397: }else {
398: vect = (VECT)t;
399: }
1.1 noro 400: if ( argc(arg) == 1 )
401: qsort(BDY(vect),vect->len,sizeof(Obj),(int (*)(const void *,const void *))comp_obj);
402: else {
403: p = (P)ARG1(arg);
404: if ( !p || OID(p)!=2 )
405: error("qsort : invalid argument");
406: v = VR(p);
1.34 noro 407: gen_searchf(NAME(v),&func);
408: if ( !func ) {
409: if ( (int)v->attr != V_SR )
410: error("qsort : no such function");
411: func = (FUNC)v->priv;
412: }
413: generic_comp_obj_func = func;
1.60 ohara 414: MKNODE(n,0,0); MKNODE(generic_comp_obj_arg,0,n);
415: generic_comp_obj_option = current_option;
1.1 noro 416: qsort(BDY(vect),vect->len,sizeof(Obj),(int (*)(const void *,const void *))generic_comp_obj);
417: }
1.35 ohara 418: if (OID(t) == O_LIST) {
419: a = BDY(vect);
420: for ( i = len - 1, n = 0; i >= 0; i-- ) {
421: MKNODE(n1,a[i],n); n = n1;
422: }
1.46 saito 423: MKLIST(*rp,n);
1.35 ohara 424: }else {
1.46 saito 425: *rp = (LIST)vect;
1.35 ohara 426: }
1.1 noro 427: }
428:
1.24 noro 429: void PNBmul_gf2n(NODE arg,GF2N *rp)
1.1 noro 430: {
431: GF2N a,b;
432: GF2MAT mat;
433: int n,w;
434: unsigned int *ab,*bb;
435: UP2 r;
436:
437: a = (GF2N)ARG0(arg);
438: b = (GF2N)ARG1(arg);
439: mat = (GF2MAT)ARG2(arg);
440: if ( !a || !b )
441: *rp = 0;
442: else {
443: n = mat->row;
444: w = (n+BSH-1)/BSH;
445:
446: ab = (unsigned int *)ALLOCA(w*sizeof(unsigned int));
447: bzero((char *)ab,w*sizeof(unsigned int));
448: bcopy(a->body->b,ab,(a->body->w)*sizeof(unsigned int));
449:
450: bb = (unsigned int *)ALLOCA(w*sizeof(unsigned int));
451: bzero((char *)bb,w*sizeof(unsigned int));
452: bcopy(b->body->b,bb,(b->body->w)*sizeof(unsigned int));
453:
454: NEWUP2(r,w);
455: bzero((char *)r->b,w*sizeof(unsigned int));
456: mul_nb(mat,ab,bb,r->b);
457: r->w = w;
458: _adjup2(r);
459: if ( !r->w )
460: *rp = 0;
461: else
462: MKGF2N(r,*rp);
463: }
464: }
465:
1.24 noro 466: void Pmul_vect_mat_gf2n(NODE arg,GF2N *rp)
1.1 noro 467: {
468: GF2N a;
469: GF2MAT mat;
470: int n,w;
471: unsigned int *b;
472: UP2 r;
473:
474: a = (GF2N)ARG0(arg);
475: mat = (GF2MAT)ARG1(arg);
476: if ( !a )
477: *rp = 0;
478: else {
479: n = mat->row;
480: w = (n+BSH-1)/BSH;
481: b = (unsigned int *)ALLOCA(w*sizeof(unsigned int));
482: bzero((char *)b,w*sizeof(unsigned int));
483: bcopy(a->body->b,b,(a->body->w)*sizeof(unsigned int));
484: NEWUP2(r,w);
485: bzero((char *)r->b,w*sizeof(unsigned int));
486: mulgf2vectmat(mat->row,b,mat->body,r->b);
487: r->w = w;
488: _adjup2(r);
489: if ( !r->w )
490: *rp = 0;
491: else {
492: MKGF2N(r,*rp);
493: }
494: }
495: }
496:
1.24 noro 497: void Pbconvmat_gf2n(NODE arg,LIST *rp)
1.1 noro 498: {
499: P p0,p1;
500: int to;
501: GF2MAT p01,p10;
502: GF2N root;
503: NODE n0,n1;
504:
505: p0 = (P)ARG0(arg);
506: p1 = (P)ARG1(arg);
507: to = ARG2(arg)?1:0;
508: if ( argc(arg) == 4 ) {
509: root = (GF2N)ARG3(arg);
510: compute_change_of_basis_matrix_with_root(p0,p1,to,root,&p01,&p10);
511: } else
512: compute_change_of_basis_matrix(p0,p1,to,&p01,&p10);
513: MKNODE(n1,p10,0); MKNODE(n0,p01,n1);
514: MKLIST(*rp,n0);
515: }
516:
1.24 noro 517: void Pmulmat_gf2n(NODE arg,GF2MAT *rp)
1.1 noro 518: {
519: GF2MAT m;
520:
521: if ( !compute_multiplication_matrix((P)ARG0(arg),&m) )
522: error("mulmat_gf2n : input is not a normal polynomial");
523: *rp = m;
524: }
525:
1.24 noro 526: void Psepmat_destructive(NODE arg,LIST *rp)
1.1 noro 527: {
528: MAT mat,mat1;
529: int i,j,row,col;
530: Q **a,**a1;
531: Q ent;
532: N nm,mod,rem,quo;
533: int sgn;
534: NODE n0,n1;
535:
536: mat = (MAT)ARG0(arg); mod = NM((Q)ARG1(arg));
537: row = mat->row; col = mat->col;
538: MKMAT(mat1,row,col);
539: a = (Q **)mat->body; a1 = (Q **)mat1->body;
540: for ( i = 0; i < row; i++ )
541: for ( j = 0; j < col; j++ ) {
542: ent = a[i][j];
543: if ( !ent )
544: continue;
545: nm = NM(ent);
546: sgn = SGN(ent);
547: divn(nm,mod,&quo,&rem);
548: /* if ( quo != nm && rem != nm ) */
1.61 noro 549: /* GCFREE(nm); */
550: /* GCFREE(ent); */
1.1 noro 551: NTOQ(rem,sgn,a[i][j]); NTOQ(quo,sgn,a1[i][j]);
552: }
553: MKNODE(n1,mat1,0); MKNODE(n0,mat,n1);
554: MKLIST(*rp,n0);
555: }
556:
1.24 noro 557: void Psepvect(NODE arg,VECT *rp)
1.1 noro 558: {
559: sepvect((VECT)ARG0(arg),QTOS((Q)ARG1(arg)),rp);
560: }
561:
1.24 noro 562: void sepvect(VECT v,int d,VECT *rp)
1.1 noro 563: {
564: int i,j,k,n,q,q1,r;
565: pointer *pv,*pw,*pu;
566: VECT w,u;
567:
568: n = v->len;
569: if ( d > n )
570: d = n;
571: q = n/d; r = n%d; q1 = q+1;
572: MKVECT(w,d); *rp = w;
573: pv = BDY(v); pw = BDY(w); k = 0;
574: for ( i = 0; i < r; i++ ) {
575: MKVECT(u,q1); pw[i] = (pointer)u;
576: for ( pu = BDY(u), j = 0; j < q1; j++, k++ )
577: pu[j] = pv[k];
578: }
579: for ( ; i < d; i++ ) {
580: MKVECT(u,q); pw[i] = (pointer)u;
581: for ( pu = BDY(u), j = 0; j < q; j++, k++ )
582: pu[j] = pv[k];
583: }
584: }
585:
1.24 noro 586: void Pnewvect(NODE arg,VECT *rp)
1.1 noro 587: {
588: int len,i,r;
589: VECT vect;
590: pointer *vb;
591: LIST list;
592: NODE tn;
593:
594: asir_assert(ARG0(arg),O_N,"newvect");
595: len = QTOS((Q)ARG0(arg));
1.5 noro 596: if ( len < 0 )
1.1 noro 597: error("newvect : invalid size");
598: MKVECT(vect,len);
599: if ( argc(arg) == 2 ) {
600: list = (LIST)ARG1(arg);
601: asir_assert(list,O_LIST,"newvect");
1.56 ohara 602: #if 0
1.1 noro 603: for ( r = 0, tn = BDY(list); tn; r++, tn = NEXT(tn) );
604: if ( r > len ) {
605: *rp = vect;
606: return;
607: }
1.56 ohara 608: #endif
1.1 noro 609: for ( i = 0, tn = BDY(list), vb = BDY(vect); tn; i++, tn = NEXT(tn) )
610: vb[i] = (pointer)BDY(tn);
611: }
612: *rp = vect;
1.14 noro 613: }
614:
1.28 saito 615: void Pvect(NODE arg,VECT *rp) {
1.62 ohara 616: int len,i;
1.28 saito 617: VECT vect;
618: pointer *vb;
619: NODE tn;
620:
621: if ( !arg ) {
622: *rp =0;
623: return;
624: }
625:
626: for (len = 0, tn = arg; tn; tn = NEXT(tn), len++);
1.29 saito 627: if ( len == 1 ) {
628: if ( ARG0(arg) != 0 ) {
629: switch ( OID(ARG0(arg)) ) {
630: case O_VECT:
631: *rp = ARG0(arg);
632: return;
633: case O_LIST:
634: for ( len = 0, tn = ARG0(arg); tn; tn = NEXT(tn), len++ );
635: MKVECT(vect,len-1);
636: for ( i = 0, tn = BDY((LIST)ARG0(arg)), vb =BDY(vect);
637: tn; i++, tn = NEXT(tn) )
638: vb[i] = (pointer)BDY(tn);
639: *rp=vect;
640: return;
641: }
642: }
643: }
1.28 saito 644: MKVECT(vect,len);
645: for ( i = 0, tn = arg, vb = BDY(vect); tn; i++, tn = NEXT(tn) )
646: vb[i] = (pointer)BDY(tn);
647: *rp = vect;
648: }
649:
1.24 noro 650: void Pexponent_vector(NODE arg,DP *rp)
1.14 noro 651: {
652: nodetod(arg,rp);
1.9 noro 653: }
654:
1.24 noro 655: void Pnewbytearray(NODE arg,BYTEARRAY *rp)
1.9 noro 656: {
657: int len,i,r;
658: BYTEARRAY array;
659: unsigned char *vb;
1.10 noro 660: char *str;
1.9 noro 661: LIST list;
662: NODE tn;
1.51 noro 663: int ac;
664: struct stat sbuf;
665: char *fname;
666: FILE *fp;
667:
668: ac = argc(arg);
669: if ( ac == 1 ) {
1.57 noro 670: if ( !OID((Obj)ARG0(arg)) ) error("newbytearray : invalid argument");
671: switch ( OID((Obj)ARG0(arg)) ) {
672: case O_STR:
673: fname = BDY((STRING)ARG0(arg));
674: fp = fopen(fname,"rb");
675: if ( !fp ) error("newbytearray : fopen failed");
676: if ( stat(fname,&sbuf) < 0 )
677: error("newbytearray : stat failed");
678: len = sbuf.st_size;
679: MKBYTEARRAY(array,len);
680: fread(BDY(array),len,sizeof(char),fp);
681: break;
682: case O_N:
683: if ( !RATN(ARG0(arg)) )
684: error("newbytearray : invalid argument");
685: len = QTOS((Q)ARG0(arg));
686: if ( len < 0 )
687: error("newbytearray : invalid size");
688: MKBYTEARRAY(array,len);
689: break;
690: default:
691: error("newbytearray : invalid argument");
692: }
1.51 noro 693: } else if ( ac == 2 ) {
694: asir_assert(ARG0(arg),O_N,"newbytearray");
695: len = QTOS((Q)ARG0(arg));
696: if ( len < 0 )
697: error("newbytearray : invalid size");
698: MKBYTEARRAY(array,len);
1.10 noro 699: if ( !ARG1(arg) )
700: error("newbytearray : invalid initialization");
701: switch ( OID((Obj)ARG1(arg)) ) {
702: case O_LIST:
703: list = (LIST)ARG1(arg);
704: asir_assert(list,O_LIST,"newbytearray");
705: for ( r = 0, tn = BDY(list); tn; r++, tn = NEXT(tn) );
706: if ( r <= len ) {
707: for ( i = 0, tn = BDY(list), vb = BDY(array); tn;
708: i++, tn = NEXT(tn) )
709: vb[i] = (unsigned char)QTOS((Q)BDY(tn));
710: }
711: break;
712: case O_STR:
713: str = BDY((STRING)ARG1(arg));
714: r = strlen(str);
715: if ( r <= len )
716: bcopy(str,BDY(array),r);
717: break;
718: default:
719: if ( !ARG1(arg) )
720: error("newbytearray : invalid initialization");
1.9 noro 721: }
1.51 noro 722: } else
723: error("newbytearray : invalid argument");
1.9 noro 724: *rp = array;
1.49 noro 725: }
726:
727: #define MEMORY_GETPOINT(a,len,x,y) (((a)[(len)*(y)+((x)>>3)])&(1<<((x)&7)))
728:
729: void Pmemoryplot_to_coord(NODE arg,LIST *rp)
730: {
731: int len,blen,y,i,j;
1.51 noro 732: unsigned char *a;
1.49 noro 733: NODE r0,r,n;
734: LIST l;
735: BYTEARRAY ba;
736: Q iq,jq;
737:
738: asir_assert(ARG0(arg),O_LIST,"memoryplot_to_coord");
739: arg = BDY((LIST)ARG0(arg));
740: len = QTOS((Q)ARG0(arg));
741: blen = (len+7)/8;
742: y = QTOS((Q)ARG1(arg));
743: ba = (BYTEARRAY)ARG2(arg); a = ba->body;
744: r0 = 0;
745: for ( j = 0; j < y; j++ )
746: for ( i = 0; i < len; i++ )
747: if ( MEMORY_GETPOINT(a,blen,i,j) ) {
748: NEXTNODE(r0,r);
749: STOQ(i,iq); STOQ(j,jq);
750: n = mknode(2,iq,jq);
751: MKLIST(l,n);
752: BDY(r) = l;
753: }
754: if ( r0 ) NEXT(r) = 0;
755: MKLIST(*rp,r0);
1.1 noro 756: }
757:
1.24 noro 758: void Pnewmat(NODE arg,MAT *rp)
1.1 noro 759: {
760: int row,col;
761: int i,j,r,c;
762: NODE tn,sn;
763: MAT m;
764: pointer **mb;
765: LIST list;
766:
767: asir_assert(ARG0(arg),O_N,"newmat");
768: asir_assert(ARG1(arg),O_N,"newmat");
769: row = QTOS((Q)ARG0(arg)); col = QTOS((Q)ARG1(arg));
1.5 noro 770: if ( row < 0 || col < 0 )
1.1 noro 771: error("newmat : invalid size");
772: MKMAT(m,row,col);
773: if ( argc(arg) == 3 ) {
774: list = (LIST)ARG2(arg);
775: asir_assert(list,O_LIST,"newmat");
776: for ( r = 0, c = 0, tn = BDY(list); tn; r++, tn = NEXT(tn) ) {
777: for ( j = 0, sn = BDY((LIST)BDY(tn)); sn; j++, sn = NEXT(sn) );
778: c = MAX(c,j);
779: }
780: if ( (r > row) || (c > col) ) {
781: *rp = m;
782: return;
783: }
784: for ( i = 0, tn = BDY(list), mb = BDY(m); tn; i++, tn = NEXT(tn) ) {
785: asir_assert(BDY(tn),O_LIST,"newmat");
786: for ( j = 0, sn = BDY((LIST)BDY(tn)); sn; j++, sn = NEXT(sn) )
787: mb[i][j] = (pointer)BDY(sn);
788: }
789: }
1.28 saito 790: *rp = m;
791: }
792:
793: void Pmat(NODE arg, MAT *rp)
794: {
795: int row,col;
1.29 saito 796: int i;
1.28 saito 797: MAT m;
798: pointer **mb;
1.29 saito 799: pointer *ent;
1.28 saito 800: NODE tn, sn;
1.29 saito 801: VECT v;
1.28 saito 802:
803: if ( !arg ) {
804: *rp =0;
805: return;
806: }
807:
808: for (row = 0, tn = arg; tn; tn = NEXT(tn), row++);
1.30 saito 809: if ( row == 1 ) {
810: if ( OID(ARG0(arg)) == O_MAT ) {
811: *rp=ARG0(arg);
812: return;
813: } else if ( !(OID(ARG0(arg)) == O_LIST || OID(ARG0(arg)) == O_VECT)) {
814: error("mat : invalid argument");
815: }
816: }
1.29 saito 817: if ( OID(ARG0(arg)) == O_VECT ) {
818: v = ARG0(arg);
819: col = v->len;
820: } else if ( OID(ARG0(arg)) == O_LIST ) {
821: for (col = 0, tn = BDY((LIST)ARG0(arg)); tn ; tn = NEXT(tn), col++);
1.30 saito 822: } else {
823: error("mat : invalid argument");
1.29 saito 824: }
825:
1.28 saito 826: MKMAT(m,row,col);
1.29 saito 827: for (row = 0, tn = arg, mb = BDY(m); tn; tn = NEXT(tn), row++) {
828: if ( BDY(tn) == 0 ) {
829: error("mat : invalid argument");
830: } else if ( OID(BDY(tn)) == O_VECT ) {
831: v = tn->body;
832: ent = BDY(v);
833: for (i = 0; i < v->len; i++ ) mb[row][i] = (Obj)ent[i];
834: } else if ( OID(BDY(tn)) == O_LIST ) {
835: for (col = 0, sn = BDY((LIST)BDY(tn)); sn; col++, sn = NEXT(sn) )
836: mb[row][col] = (pointer)BDY(sn);
837: } else {
838: error("mat : invalid argument");
839: }
840: }
841: *rp = m;
842: }
843:
844: void Pmatc(NODE arg, MAT *rp)
845: {
846: int row,col;
847: int i;
848: MAT m;
849: pointer **mb;
850: pointer *ent;
851: NODE tn, sn;
852: VECT v;
853:
854: if ( !arg ) {
855: *rp =0;
856: return;
857: }
858:
859: for (col = 0, tn = arg; tn; tn = NEXT(tn), col++);
1.30 saito 860: if ( col == 1 ) {
861: if ( OID(ARG0(arg)) == O_MAT ) {
862: *rp=ARG0(arg);
863: return;
864: } else if ( !(OID(ARG0(arg)) == O_LIST || OID(ARG0(arg)) == O_VECT)) {
865: error("matc : invalid argument");
866: }
867: }
1.29 saito 868: if ( OID(ARG0(arg)) == O_VECT ) {
869: v = ARG0(arg);
870: row = v->len;
871: } else if ( OID(ARG0(arg)) == O_LIST ) {
872: for (row = 0, tn = BDY((LIST)ARG0(arg)); tn ; tn = NEXT(tn), row++);
1.30 saito 873: } else {
874: error("matc : invalid argument");
1.29 saito 875: }
876:
877: MKMAT(m,row,col);
878: for (col = 0, tn = arg, mb = BDY(m); tn; tn = NEXT(tn), col++) {
879: if ( BDY(tn) == 0 ) {
880: error("matc : invalid argument");
881: } else if ( OID(BDY(tn)) == O_VECT ) {
882: v = tn->body;
883: ent = BDY(v);
884: for (i = 0; i < v->len; i++ ) mb[i][col] = (Obj)ent[i];
885: } else if ( OID(BDY(tn)) == O_LIST ) {
886: for (row = 0, sn = BDY((LIST)BDY(tn)); sn; row++, sn = NEXT(sn) )
887: mb[row][col] = (pointer)BDY(sn);
888: } else {
889: error("matc : invalid argument");
890: }
891: }
1.1 noro 892: *rp = m;
893: }
894:
1.24 noro 895: void Pvtol(NODE arg,LIST *rp)
1.1 noro 896: {
897: NODE n,n1;
898: VECT v;
899: pointer *a;
900: int len,i;
901:
1.66 ohara 902: if ( OID(ARG0(arg)) == O_LIST ) {
903: *rp = ARG0(arg);
904: return;
905: }
1.1 noro 906: asir_assert(ARG0(arg),O_VECT,"vtol");
907: v = (VECT)ARG0(arg); len = v->len; a = BDY(v);
908: for ( i = len - 1, n = 0; i >= 0; i-- ) {
909: MKNODE(n1,a[i],n); n = n1;
910: }
911: MKLIST(*rp,n);
1.33 noro 912: }
913:
914: void Pltov(NODE arg,VECT *rp)
915: {
916: NODE n;
1.66 ohara 917: VECT v,v0;
1.33 noro 918: int len,i;
919:
1.66 ohara 920: if ( OID(ARG0(arg)) == O_VECT ) {
921: v0 = (VECT)ARG0(arg); len = v0->len;
922: MKVECT(v,len);
923: for ( i = 0; i < len; i++ ) {
924: BDY(v)[i] = BDY(v0)[i];
925: }
926: *rp = v;
927: return;
928: }
1.33 noro 929: asir_assert(ARG0(arg),O_LIST,"ltov");
930: n = (NODE)BDY((LIST)ARG0(arg));
931: len = length(n);
932: MKVECT(v,len);
933: for ( i = 0; i < len; i++, n = NEXT(n) )
934: BDY(v)[i] = BDY(n);
935: *rp = v;
1.1 noro 936: }
937:
1.24 noro 938: void Premainder(NODE arg,Obj *rp)
1.1 noro 939: {
940: Obj a;
941: VECT v,w;
942: MAT m,l;
943: pointer *vb,*wb;
944: pointer **mb,**lb;
945: int id,i,j,n,row,col,t,smd,sgn;
946: Q md,q;
947:
948: a = (Obj)ARG0(arg); md = (Q)ARG1(arg);
949: if ( !a )
950: *rp = 0;
951: else {
952: id = OID(a);
953: switch ( id ) {
954: case O_N:
955: case O_P:
956: cmp(md,(P)a,(P *)rp); break;
957: case O_VECT:
958: smd = QTOS(md);
959: v = (VECT)a; n = v->len; vb = v->body;
960: MKVECT(w,n); wb = w->body;
961: for ( i = 0; i < n; i++ ) {
962: if ( q = (Q)vb[i] ) {
963: sgn = SGN(q); t = rem(NM(q),smd);
964: STOQ(t,q);
965: if ( q )
966: SGN(q) = sgn;
967: }
968: wb[i] = (pointer)q;
969: }
970: *rp = (Obj)w;
971: break;
972: case O_MAT:
973: m = (MAT)a; row = m->row; col = m->col; mb = m->body;
974: MKMAT(l,row,col); lb = l->body;
975: for ( i = 0; i < row; i++ )
976: for ( j = 0, vb = mb[i], wb = lb[i]; j < col; j++ )
977: cmp(md,(P)vb[j],(P *)&wb[j]);
978: *rp = (Obj)l;
979: break;
980: default:
981: error("remainder : invalid argument");
982: }
983: }
984: }
985:
1.24 noro 986: void Psremainder(NODE arg,Obj *rp)
1.1 noro 987: {
988: Obj a;
989: VECT v,w;
990: MAT m,l;
991: pointer *vb,*wb;
992: pointer **mb,**lb;
993: unsigned int t,smd;
994: int id,i,j,n,row,col;
995: Q md,q;
996:
997: a = (Obj)ARG0(arg); md = (Q)ARG1(arg);
998: if ( !a )
999: *rp = 0;
1000: else {
1001: id = OID(a);
1002: switch ( id ) {
1003: case O_N:
1004: case O_P:
1005: cmp(md,(P)a,(P *)rp); break;
1006: case O_VECT:
1007: smd = QTOS(md);
1008: v = (VECT)a; n = v->len; vb = v->body;
1009: MKVECT(w,n); wb = w->body;
1010: for ( i = 0; i < n; i++ ) {
1011: if ( q = (Q)vb[i] ) {
1012: t = (unsigned int)rem(NM(q),smd);
1013: if ( SGN(q) < 0 )
1014: t = (smd - t) % smd;
1015: UTOQ(t,q);
1016: }
1017: wb[i] = (pointer)q;
1018: }
1019: *rp = (Obj)w;
1020: break;
1021: case O_MAT:
1022: m = (MAT)a; row = m->row; col = m->col; mb = m->body;
1023: MKMAT(l,row,col); lb = l->body;
1024: for ( i = 0; i < row; i++ )
1025: for ( j = 0, vb = mb[i], wb = lb[i]; j < col; j++ )
1026: cmp(md,(P)vb[j],(P *)&wb[j]);
1027: *rp = (Obj)l;
1028: break;
1029: default:
1030: error("remainder : invalid argument");
1031: }
1032: }
1033: }
1034:
1.24 noro 1035: void Psize(NODE arg,LIST *rp)
1.1 noro 1036: {
1037:
1038: int n,m;
1039: Q q;
1040: NODE t,s;
1041:
1042: if ( !ARG0(arg) )
1043: t = 0;
1044: else {
1045: switch (OID(ARG0(arg))) {
1046: case O_VECT:
1047: n = ((VECT)ARG0(arg))->len;
1048: STOQ(n,q); MKNODE(t,q,0);
1049: break;
1050: case O_MAT:
1051: n = ((MAT)ARG0(arg))->row; m = ((MAT)ARG0(arg))->col;
1.43 saito 1052: STOQ(m,q); MKNODE(s,q,0); STOQ(n,q); MKNODE(t,q,s);
1053: break;
1054: case O_IMAT:
1055: n = ((IMAT)ARG0(arg))->row; m = ((IMAT)ARG0(arg))->col;
1.1 noro 1056: STOQ(m,q); MKNODE(s,q,0); STOQ(n,q); MKNODE(t,q,s);
1057: break;
1058: default:
1059: error("size : invalid argument"); break;
1060: }
1061: }
1062: MKLIST(*rp,t);
1063: }
1064:
1.24 noro 1065: void Pdet(NODE arg,P *rp)
1.1 noro 1066: {
1067: MAT m;
1068: int n,i,j,mod;
1069: P d;
1070: P **mat,**w;
1071:
1072: m = (MAT)ARG0(arg);
1073: asir_assert(m,O_MAT,"det");
1074: if ( m->row != m->col )
1075: error("det : non-square matrix");
1076: else if ( argc(arg) == 1 )
1077: detp(CO,(P **)BDY(m),m->row,rp);
1078: else {
1079: n = m->row; mod = QTOS((Q)ARG1(arg)); mat = (P **)BDY(m);
1080: w = (P **)almat_pointer(n,n);
1081: for ( i = 0; i < n; i++ )
1082: for ( j = 0; j < n; j++ )
1083: ptomp(mod,mat[i][j],&w[i][j]);
1084: detmp(CO,mod,w,n,&d);
1085: mptop(d,rp);
1.23 noro 1086: }
1087: }
1088:
1.24 noro 1089: void Pinvmat(NODE arg,LIST *rp)
1.23 noro 1090: {
1091: MAT m,r;
1092: int n,i,j,mod;
1093: P dn;
1094: P **mat,**imat,**w;
1095: NODE nd;
1096:
1097: m = (MAT)ARG0(arg);
1098: asir_assert(m,O_MAT,"invmat");
1099: if ( m->row != m->col )
1100: error("invmat : non-square matrix");
1101: else if ( argc(arg) == 1 ) {
1102: n = m->row;
1103: invmatp(CO,(P **)BDY(m),n,&imat,&dn);
1104: NEWMAT(r); r->row = n; r->col = n; r->body = (pointer **)imat;
1105: nd = mknode(2,r,dn);
1106: MKLIST(*rp,nd);
1107: } else {
1108: n = m->row; mod = QTOS((Q)ARG1(arg)); mat = (P **)BDY(m);
1109: w = (P **)almat_pointer(n,n);
1110: for ( i = 0; i < n; i++ )
1111: for ( j = 0; j < n; j++ )
1112: ptomp(mod,mat[i][j],&w[i][j]);
1113: #if 0
1114: detmp(CO,mod,w,n,&d);
1115: mptop(d,rp);
1116: #else
1117: error("not implemented yet");
1118: #endif
1.1 noro 1119: }
1.25 noro 1120: }
1121:
1122: /*
1123: input : a row x col matrix A
1124: A[I] <-> A[I][0]*x_0+A[I][1]*x_1+...
1125:
1.45 noro 1126: output : [B,D,R,C]
1.25 noro 1127: B : a rank(A) x col-rank(A) matrix
1.45 noro 1128: D : the denominator
1.25 noro 1129: R : a vector of length rank(A)
1130: C : a vector of length col-rank(A)
1.45 noro 1131: B[I] <-> D*x_{R[I]}+B[I][0]x_{C[0]}+B[I][1]x_{C[1]}+...
1.25 noro 1132: */
1133:
1134: void Pgeneric_gauss_elim(NODE arg,LIST *rp)
1135: {
1.48 noro 1136: NODE n0,opt,p;
1.25 noro 1137: MAT m,nm;
1138: int *ri,*ci;
1139: VECT rind,cind;
1140: Q dn,q;
1.62 ohara 1141: int i,row,col,t,rank;
1.48 noro 1142: int is_hensel = 0;
1143: char *key;
1144: Obj value;
1145:
1146: if ( current_option ) {
1147: for ( opt = current_option; opt; opt = NEXT(opt) ) {
1148: p = BDY((LIST)BDY(opt));
1149: key = BDY((STRING)BDY(p));
1150: value = (Obj)BDY(NEXT(p));
1151: if ( !strcmp(key,"hensel") && value ) {
1152: is_hensel = value ? 1 : 0;
1153: break;
1154: }
1155: }
1156: }
1.25 noro 1157: asir_assert(ARG0(arg),O_MAT,"generic_gauss_elim");
1158: m = (MAT)ARG0(arg);
1159: row = m->row; col = m->col;
1.48 noro 1160: if ( is_hensel )
1161: rank = generic_gauss_elim_hensel(m,&nm,&dn,&ri,&ci);
1162: else
1163: rank = generic_gauss_elim(m,&nm,&dn,&ri,&ci);
1.25 noro 1164: t = col-rank;
1165: MKVECT(rind,rank);
1166: MKVECT(cind,t);
1167: for ( i = 0; i < rank; i++ ) {
1168: STOQ(ri[i],q);
1169: BDY(rind)[i] = (pointer)q;
1170: }
1171: for ( i = 0; i < t; i++ ) {
1172: STOQ(ci[i],q);
1173: BDY(cind)[i] = (pointer)q;
1174: }
1175: n0 = mknode(4,nm,dn,rind,cind);
1176: MKLIST(*rp,n0);
1.1 noro 1177: }
1178:
1.69 noro 1179: void Pindep_rows_mod(NODE arg,VECT *rp)
1180: {
1181: MAT m,mat;
1182: VECT rind;
1183: Q **tmat;
1184: int **wmat,**row0;
1185: Q *rib;
1186: int *rowstat,*p;
1187: Q q;
1188: int md,i,j,k,l,row,col,t,rank;
1189:
1190: asir_assert(ARG0(arg),O_MAT,"indep_rows_mod");
1191: asir_assert(ARG1(arg),O_N,"indep_rows_mod");
1192: m = (MAT)ARG0(arg); md = QTOS((Q)ARG1(arg));
1193: row = m->row; col = m->col; tmat = (Q **)m->body;
1194: wmat = (int **)almat(row,col);
1195:
1196: row0 = (int **)ALLOCA(row*sizeof(int *));
1197: for ( i = 0; i < row; i++ ) row0[i] = wmat[i];
1198:
1199: rowstat = (int *)MALLOC_ATOMIC(row*sizeof(int));
1200: for ( i = 0; i < row; i++ )
1201: for ( j = 0; j < col; j++ )
1202: if ( q = (Q)tmat[i][j] ) {
1203: t = rem(NM(q),md);
1204: if ( t && SGN(q) < 0 )
1205: t = (md - t) % md;
1206: wmat[i][j] = t;
1207: } else
1208: wmat[i][j] = 0;
1209: rank = indep_rows_mod(wmat,row,col,md,rowstat);
1210:
1211: MKVECT(rind,rank);
1212: rib = (Q *)rind->body;
1213: for ( j = 0; j < rank; j++ ) {
1214: STOQ(rowstat[j],rib[j]);
1215: }
1216: *rp = rind;
1217: }
1218:
1.1 noro 1219: /*
1220: input : a row x col matrix A
1221: A[I] <-> A[I][0]*x_0+A[I][1]*x_1+...
1222:
1223: output : [B,R,C]
1224: B : a rank(A) x col-rank(A) matrix
1225: R : a vector of length rank(A)
1226: C : a vector of length col-rank(A)
1.47 noro 1227: RN : a vector of length rank(A) indicating useful rows
1228:
1.1 noro 1229: B[I] <-> x_{R[I]}+B[I][0]x_{C[0]}+B[I][1]x_{C[1]}+...
1230: */
1231:
1.24 noro 1232: void Pgeneric_gauss_elim_mod(NODE arg,LIST *rp)
1.1 noro 1233: {
1234: NODE n0;
1235: MAT m,mat;
1.47 noro 1236: VECT rind,cind,rnum;
1.1 noro 1237: Q **tmat;
1.47 noro 1238: int **wmat,**row0;
1239: Q *rib,*cib,*rnb;
1240: int *colstat,*p;
1.1 noro 1241: Q q;
1.24 noro 1242: int md,i,j,k,l,row,col,t,rank;
1.1 noro 1243:
1244: asir_assert(ARG0(arg),O_MAT,"generic_gauss_elim_mod");
1245: asir_assert(ARG1(arg),O_N,"generic_gauss_elim_mod");
1246: m = (MAT)ARG0(arg); md = QTOS((Q)ARG1(arg));
1247: row = m->row; col = m->col; tmat = (Q **)m->body;
1248: wmat = (int **)almat(row,col);
1.47 noro 1249:
1250: row0 = (int **)ALLOCA(row*sizeof(int *));
1251: for ( i = 0; i < row; i++ ) row0[i] = wmat[i];
1252:
1.1 noro 1253: colstat = (int *)MALLOC_ATOMIC(col*sizeof(int));
1254: for ( i = 0; i < row; i++ )
1255: for ( j = 0; j < col; j++ )
1256: if ( q = (Q)tmat[i][j] ) {
1257: t = rem(NM(q),md);
1258: if ( t && SGN(q) < 0 )
1259: t = (md - t) % md;
1260: wmat[i][j] = t;
1261: } else
1262: wmat[i][j] = 0;
1263: rank = generic_gauss_elim_mod(wmat,row,col,md,colstat);
1264:
1.47 noro 1265: MKVECT(rnum,rank);
1266: rnb = (Q *)rnum->body;
1267: for ( i = 0; i < rank; i++ )
1268: for ( j = 0, p = wmat[i]; j < row; j++ )
1269: if ( p == row0[j] )
1270: STOQ(j,rnb[i]);
1271:
1.1 noro 1272: MKMAT(mat,rank,col-rank);
1273: tmat = (Q **)mat->body;
1274: for ( i = 0; i < rank; i++ )
1275: for ( j = k = 0; j < col; j++ )
1276: if ( !colstat[j] ) {
1277: UTOQ(wmat[i][j],tmat[i][k]); k++;
1278: }
1279:
1280: MKVECT(rind,rank);
1281: MKVECT(cind,col-rank);
1282: rib = (Q *)rind->body; cib = (Q *)cind->body;
1283: for ( j = k = l = 0; j < col; j++ )
1284: if ( colstat[j] ) {
1285: STOQ(j,rib[k]); k++;
1286: } else {
1287: STOQ(j,cib[l]); l++;
1288: }
1.47 noro 1289: n0 = mknode(4,mat,rind,cind,rnum);
1.1 noro 1290: MKLIST(*rp,n0);
1291: }
1292:
1.24 noro 1293: void Pleqm(NODE arg,VECT *rp)
1.1 noro 1294: {
1295: MAT m;
1296: VECT vect;
1297: pointer **mat;
1298: Q *v;
1299: Q q;
1300: int **wmat;
1301: int md,i,j,row,col,t,n,status;
1302:
1303: asir_assert(ARG0(arg),O_MAT,"leqm");
1304: asir_assert(ARG1(arg),O_N,"leqm");
1305: m = (MAT)ARG0(arg); md = QTOS((Q)ARG1(arg));
1306: row = m->row; col = m->col; mat = m->body;
1307: wmat = (int **)almat(row,col);
1308: for ( i = 0; i < row; i++ )
1309: for ( j = 0; j < col; j++ )
1310: if ( q = (Q)mat[i][j] ) {
1311: t = rem(NM(q),md);
1312: if ( SGN(q) < 0 )
1313: t = (md - t) % md;
1314: wmat[i][j] = t;
1315: } else
1316: wmat[i][j] = 0;
1317: status = gauss_elim_mod(wmat,row,col,md);
1318: if ( status < 0 )
1319: *rp = 0;
1320: else if ( status > 0 )
1321: *rp = (VECT)ONE;
1322: else {
1323: n = col - 1;
1324: MKVECT(vect,n);
1325: for ( i = 0, v = (Q *)vect->body; i < n; i++ ) {
1326: t = (md-wmat[i][n])%md; STOQ(t,v[i]);
1327: }
1328: *rp = vect;
1329: }
1330: }
1331:
1.24 noro 1332: int gauss_elim_mod(int **mat,int row,int col,int md)
1.1 noro 1333: {
1334: int i,j,k,inv,a,n;
1335: int *t,*pivot;
1336:
1337: n = col - 1;
1338: for ( j = 0; j < n; j++ ) {
1339: for ( i = j; i < row && !mat[i][j]; i++ );
1340: if ( i == row )
1341: return 1;
1342: if ( i != j ) {
1343: t = mat[i]; mat[i] = mat[j]; mat[j] = t;
1344: }
1345: pivot = mat[j];
1346: inv = invm(pivot[j],md);
1347: for ( k = j; k <= n; k++ ) {
1348: /* pivot[k] = dmar(pivot[k],inv,0,md); */
1349: DMAR(pivot[k],inv,0,md,pivot[k])
1350: }
1351: for ( i = 0; i < row; i++ ) {
1352: t = mat[i];
1353: if ( i != j && (a = t[j]) )
1354: for ( k = j, a = md - a; k <= n; k++ ) {
1.8 noro 1355: unsigned int tk;
1.1 noro 1356: /* t[k] = dmar(pivot[k],a,t[k],md); */
1.8 noro 1357: DMAR(pivot[k],a,t[k],md,tk)
1358: t[k] = tk;
1.1 noro 1359: }
1360: }
1361: }
1362: for ( i = n; i < row && !mat[i][n]; i++ );
1363: if ( i == row )
1364: return 0;
1365: else
1366: return -1;
1367: }
1368:
1.4 noro 1369: struct oEGT eg_mod,eg_elim,eg_elim1,eg_elim2,eg_chrem,eg_gschk,eg_intrat,eg_symb;
1.31 noro 1370: struct oEGT eg_conv;
1.1 noro 1371:
1.24 noro 1372: int generic_gauss_elim(MAT mat,MAT *nm,Q *dn,int **rindp,int **cindp)
1.1 noro 1373: {
1374: int **wmat;
1375: Q **bmat;
1376: N **tmat;
1377: Q *bmi;
1378: N *tmi;
1379: Q q;
1380: int *wmi;
1381: int *colstat,*wcolstat,*rind,*cind;
1382: int row,col,ind,md,i,j,k,l,t,t1,rank,rank0,inv;
1383: N m1,m2,m3,s,u;
1384: MAT r,crmat;
1385: struct oEGT tmp0,tmp1;
1386: struct oEGT eg_mod_split,eg_elim_split,eg_chrem_split;
1387: struct oEGT eg_intrat_split,eg_gschk_split;
1388: int ret;
1389:
1390: init_eg(&eg_mod_split); init_eg(&eg_chrem_split);
1391: init_eg(&eg_elim_split); init_eg(&eg_intrat_split);
1392: init_eg(&eg_gschk_split);
1393: bmat = (Q **)mat->body;
1394: row = mat->row; col = mat->col;
1395: wmat = (int **)almat(row,col);
1396: colstat = (int *)MALLOC_ATOMIC(col*sizeof(int));
1397: wcolstat = (int *)MALLOC_ATOMIC(col*sizeof(int));
1398: for ( ind = 0; ; ind++ ) {
1.11 noro 1399: if ( DP_Print ) {
1.2 noro 1400: fprintf(asir_out,"."); fflush(asir_out);
1401: }
1.12 noro 1402: md = get_lprime(ind);
1.1 noro 1403: get_eg(&tmp0);
1404: for ( i = 0; i < row; i++ )
1405: for ( j = 0, bmi = bmat[i], wmi = wmat[i]; j < col; j++ )
1406: if ( q = (Q)bmi[j] ) {
1407: t = rem(NM(q),md);
1408: if ( t && SGN(q) < 0 )
1409: t = (md - t) % md;
1410: wmi[j] = t;
1411: } else
1412: wmi[j] = 0;
1413: get_eg(&tmp1);
1414: add_eg(&eg_mod,&tmp0,&tmp1);
1415: add_eg(&eg_mod_split,&tmp0,&tmp1);
1416: get_eg(&tmp0);
1417: rank = generic_gauss_elim_mod(wmat,row,col,md,wcolstat);
1418: get_eg(&tmp1);
1419: add_eg(&eg_elim,&tmp0,&tmp1);
1420: add_eg(&eg_elim_split,&tmp0,&tmp1);
1421: if ( !ind ) {
1422: RESET:
1423: UTON(md,m1);
1424: rank0 = rank;
1425: bcopy(wcolstat,colstat,col*sizeof(int));
1426: MKMAT(crmat,rank,col-rank);
1427: MKMAT(r,rank,col-rank); *nm = r;
1428: tmat = (N **)crmat->body;
1429: for ( i = 0; i < rank; i++ )
1430: for ( j = k = 0, tmi = tmat[i], wmi = wmat[i]; j < col; j++ )
1431: if ( !colstat[j] ) {
1432: UTON(wmi[j],tmi[k]); k++;
1433: }
1434: } else {
1435: if ( rank < rank0 ) {
1.11 noro 1436: if ( DP_Print ) {
1.1 noro 1437: fprintf(asir_out,"lower rank matrix; continuing...\n");
1.2 noro 1438: fflush(asir_out);
1439: }
1.1 noro 1440: continue;
1441: } else if ( rank > rank0 ) {
1.11 noro 1442: if ( DP_Print ) {
1.1 noro 1443: fprintf(asir_out,"higher rank matrix; resetting...\n");
1.2 noro 1444: fflush(asir_out);
1445: }
1.1 noro 1446: goto RESET;
1447: } else {
1448: for ( j = 0; (j<col) && (colstat[j]==wcolstat[j]); j++ );
1449: if ( j < col ) {
1.11 noro 1450: if ( DP_Print ) {
1.1 noro 1451: fprintf(asir_out,"inconsitent colstat; resetting...\n");
1.2 noro 1452: fflush(asir_out);
1453: }
1.1 noro 1454: goto RESET;
1455: }
1456: }
1457:
1458: get_eg(&tmp0);
1459: inv = invm(rem(m1,md),md);
1460: UTON(md,m2); muln(m1,m2,&m3);
1461: for ( i = 0; i < rank; i++ )
1462: for ( j = k = 0, tmi = tmat[i], wmi = wmat[i]; j < col; j++ )
1463: if ( !colstat[j] ) {
1464: if ( tmi[k] ) {
1465: /* f3 = f1+m1*(m1 mod m2)^(-1)*(f2 - f1 mod m2) */
1466: t = rem(tmi[k],md);
1467: if ( wmi[j] >= t )
1468: t = wmi[j]-t;
1469: else
1470: t = md-(t-wmi[j]);
1471: DMAR(t,inv,0,md,t1)
1472: UTON(t1,u);
1473: muln(m1,u,&s);
1474: addn(tmi[k],s,&u); tmi[k] = u;
1475: } else if ( wmi[j] ) {
1476: /* f3 = m1*(m1 mod m2)^(-1)*f2 */
1477: DMAR(wmi[j],inv,0,md,t)
1478: UTON(t,u);
1479: muln(m1,u,&s); tmi[k] = s;
1480: }
1481: k++;
1482: }
1483: m1 = m3;
1484: get_eg(&tmp1);
1485: add_eg(&eg_chrem,&tmp0,&tmp1);
1486: add_eg(&eg_chrem_split,&tmp0,&tmp1);
1487:
1488: get_eg(&tmp0);
1.38 noro 1489: if ( ind % F4_INTRAT_PERIOD )
1.13 noro 1490: ret = 0;
1491: else
1492: ret = intmtoratm(crmat,m1,*nm,dn);
1.1 noro 1493: get_eg(&tmp1);
1494: add_eg(&eg_intrat,&tmp0,&tmp1);
1495: add_eg(&eg_intrat_split,&tmp0,&tmp1);
1496: if ( ret ) {
1497: *rindp = rind = (int *)MALLOC_ATOMIC(rank*sizeof(int));
1498: *cindp = cind = (int *)MALLOC_ATOMIC((col-rank)*sizeof(int));
1499: for ( j = k = l = 0; j < col; j++ )
1500: if ( colstat[j] )
1501: rind[k++] = j;
1502: else
1503: cind[l++] = j;
1504: get_eg(&tmp0);
1.3 noro 1505: if ( gensolve_check(mat,*nm,*dn,rind,cind) ) {
1506: get_eg(&tmp1);
1507: add_eg(&eg_gschk,&tmp0,&tmp1);
1508: add_eg(&eg_gschk_split,&tmp0,&tmp1);
1.11 noro 1509: if ( DP_Print ) {
1.3 noro 1510: print_eg("Mod",&eg_mod_split);
1511: print_eg("Elim",&eg_elim_split);
1512: print_eg("ChRem",&eg_chrem_split);
1513: print_eg("IntRat",&eg_intrat_split);
1514: print_eg("Check",&eg_gschk_split);
1515: fflush(asir_out);
1516: }
1517: return rank;
1518: }
1519: }
1520: }
1521: }
1522: }
1523:
1.64 noro 1524: void lu_dec_cr(MAT mat,MAT lu,Q *dn,int **perm);
1525:
1.53 noro 1526: /* XXX broken */
1.64 noro 1527: void lu_dec_cr(MAT mat,MAT lu,Q *dn,int **perm)
1.53 noro 1528: {
1529: Q **a0,**b;
1530: Q *aiq;
1531: N **a;
1532: N *ai;
1533: Q q,q1,dn2,a1,q0,bik;
1534: MAT m;
1535: unsigned int md;
1536: int n,ind,i,j,rank,t,inv,t1,ret,min,k;
1537: int **w;
1538: int *wi,*rinfo0,*rinfo;
1539: N m1,m2,m3,u,s;
1540:
1541: a0 = (Q **)mat->body;
1542: n = mat->row;
1543: if ( n != mat->col )
1544: error("lu_dec_cr : non-square matrix");
1545: w = (int **)almat(n,n);
1546: MKMAT(m,n,n);
1547: a = (N **)m->body;
1548: UTON(1,m1);
1549: rinfo0 = 0;
1550: ind = 0;
1551: while ( 1 ) {
1552: md = get_lprime(ind);
1553: /* mat mod md */
1554: for ( i = 0; i < n; i++ )
1555: for ( j = 0, aiq = a0[i], wi = w[i]; j < n; j++ )
1556: if ( q = aiq[j] ) {
1557: t = rem(NM(q),md);
1558: if ( t && SGN(q) < 0 )
1559: t = (md - t) % md;
1560: wi[j] = t;
1561: } else
1562: wi[j] = 0;
1563:
1564: if ( !lu_mod((unsigned int **)w,n,md,&rinfo) ) continue;
1565: printf("."); fflush(stdout);
1566: if ( !rinfo0 )
1567: *perm = rinfo0 = rinfo;
1568: else {
1569: for ( i = 0; i < n; i++ )
1570: if ( rinfo[i] != rinfo0[i] ) break;
1571: if ( i < n ) continue;
1572: }
1573: if ( UNIN(m1) ) {
1574: for ( i = 0; i < n; i++ )
1575: for ( j = 0, ai = a[i], wi = w[i]; j < n; j++ ) {
1576: UTON(wi[j],u); ai[j] = u;
1577: }
1578: UTON(md,m1);
1579: } else {
1580: inv = invm(rem(m1,md),md);
1581: UTON(md,m2); muln(m1,m2,&m3);
1582: for ( i = 0; i < n; i++ )
1583: for ( j = 0, ai = a[i], wi = w[i]; j < n; j++ )
1584: if ( ai[i] ) {
1585: /* f3 = f1+m1*(m1 mod m2)^(-1)*(f2 - f1 mod m2) */
1586: t = rem(ai[j],md);
1587: if ( wi[j] >= t )
1588: t = wi[j]-t;
1589: else
1590: t = md-(t-wi[j]);
1591: DMAR(t,inv,0,md,t1)
1592: UTON(t1,u);
1593: muln(m1,u,&s);
1594: addn(ai[j],s,&u); ai[j] = u;
1595: } else if ( wi[j] ) {
1596: /* f3 = m1*(m1 mod m2)^(-1)*f2 */
1597: DMAR(wi[j],inv,0,md,t)
1598: UTON(t,u);
1599: muln(m1,u,&s); ai[j] = s;
1600: }
1601: m1 = m3;
1602: }
1603: if ( (++ind%8) == 0 ) {
1604: ret = intmtoratm(m,m1,lu,dn);
1605: if ( ret ) {
1606: b = (Q **)lu->body;
1607: mulq(*dn,*dn,&dn2);
1608: for ( i = 0; i < n; i++ ) {
1609: for ( j = 0; j < n; j++ ) {
1610: q = 0;
1611: min = MIN(i,j);
1612: for ( k = 0; k <= min; k++ ) {
1613: bik = k==i ? *dn : b[i][k];
1614: mulq(bik,b[k][j],&q0);
1615: addq(q,q0,&q1); q = q1;
1616: }
1617: mulq(a0[rinfo0[i]][j],dn2,&q1);
1618: if ( cmpq(q,q1) ) break;
1619: }
1620: if ( j < n ) break;
1621: }
1622: if ( i == n )
1623: return;
1624: }
1625: }
1626: }
1627: }
1628:
1.64 noro 1629: void nmat(N **m,int n)
1.53 noro 1630: {
1631: int i,j;
1632:
1633: for ( i = 0; i < n; i++ ) {
1634: for ( j = 0; j < n; j++ ) {
1635: printn(m[i][j]); printf(" ");
1636: }
1637: printf("\n");
1638: }
1639: }
1640:
1.24 noro 1641: int generic_gauss_elim_hensel(MAT mat,MAT *nmmat,Q *dn,int **rindp,int **cindp)
1.3 noro 1642: {
1643: MAT bmat,xmat;
1644: Q **a0,**a,**b,**x,**nm;
1645: Q *ai,*bi,*xi;
1646: int row,col;
1647: int **w;
1648: int *wi;
1649: int **wc;
1650: Q mdq,q,s,u;
1651: N tn;
1652: int ind,md,i,j,k,l,li,ri,rank;
1653: unsigned int t;
1654: int *cinfo,*rinfo;
1655: int *rind,*cind;
1656: int count;
1.41 noro 1657: int ret;
1658: struct oEGT eg_mul,eg_inv,eg_intrat,eg_check,tmp0,tmp1;
1.39 noro 1659: int period;
1.44 noro 1660: int *wx,*ptr;
1661: int wxsize,nsize;
1662: N wn;
1663: Q wq;
1.3 noro 1664:
1665: a0 = (Q **)mat->body;
1666: row = mat->row; col = mat->col;
1667: w = (int **)almat(row,col);
1668: for ( ind = 0; ; ind++ ) {
1.12 noro 1669: md = get_lprime(ind);
1.3 noro 1670: STOQ(md,mdq);
1671: for ( i = 0; i < row; i++ )
1672: for ( j = 0, ai = a0[i], wi = w[i]; j < col; j++ )
1673: if ( q = (Q)ai[j] ) {
1674: t = rem(NM(q),md);
1675: if ( t && SGN(q) < 0 )
1676: t = (md - t) % md;
1677: wi[j] = t;
1678: } else
1679: wi[j] = 0;
1680:
1.52 noro 1681: if ( DP_Print > 3 ) {
1.48 noro 1682: fprintf(asir_out,"LU decomposition.."); fflush(asir_out);
1683: }
1.27 noro 1684: rank = find_lhs_and_lu_mod((unsigned int **)w,row,col,md,&rinfo,&cinfo);
1.52 noro 1685: if ( DP_Print > 3 ) {
1.48 noro 1686: fprintf(asir_out,"done.\n"); fflush(asir_out);
1687: }
1.3 noro 1688: a = (Q **)almat_pointer(rank,rank); /* lhs mat */
1689: MKMAT(bmat,rank,col-rank); b = (Q **)bmat->body; /* lhs mat */
1690: for ( j = li = ri = 0; j < col; j++ )
1691: if ( cinfo[j] ) {
1692: /* the column is in lhs */
1693: for ( i = 0; i < rank; i++ ) {
1694: w[i][li] = w[i][j];
1695: a[i][li] = a0[rinfo[i]][j];
1696: }
1697: li++;
1698: } else {
1699: /* the column is in rhs */
1700: for ( i = 0; i < rank; i++ )
1701: b[i][ri] = a0[rinfo[i]][j];
1702: ri++;
1703: }
1704:
1705: /* solve Ax+B=0; A: rank x rank, B: rank x ri */
1706: MKMAT(xmat,rank,ri); x = (Q **)(xmat)->body;
1707: MKMAT(*nmmat,rank,ri); nm = (Q **)(*nmmat)->body;
1708: /* use the right part of w as work area */
1709: /* ri = col - rank */
1710: wc = (int **)almat(rank,ri);
1711: for ( i = 0; i < rank; i++ )
1712: wc[i] = w[i]+rank;
1713: *rindp = rind = (int *)MALLOC_ATOMIC(rank*sizeof(int));
1714: *cindp = cind = (int *)MALLOC_ATOMIC((ri)*sizeof(int));
1715:
1716: init_eg(&eg_mul); init_eg(&eg_inv);
1.41 noro 1717: init_eg(&eg_check); init_eg(&eg_intrat);
1.39 noro 1718: period = F4_INTRAT_PERIOD;
1.44 noro 1719: nsize = period;
1720: wxsize = rank*ri*nsize;
1721: wx = (int *)MALLOC_ATOMIC(wxsize*sizeof(int));
1722: for ( i = 0; i < wxsize; i++ ) wx[i] = 0;
1723: for ( q = ONE, count = 0; ; ) {
1.52 noro 1724: if ( DP_Print > 3 )
1.41 noro 1725: fprintf(stderr,"o");
1.3 noro 1726: /* wc = -b mod md */
1.44 noro 1727: get_eg(&tmp0);
1.3 noro 1728: for ( i = 0; i < rank; i++ )
1729: for ( j = 0, bi = b[i], wi = wc[i]; j < ri; j++ )
1730: if ( u = (Q)bi[j] ) {
1731: t = rem(NM(u),md);
1732: if ( t && SGN(u) > 0 )
1733: t = (md - t) % md;
1734: wi[j] = t;
1735: } else
1736: wi[j] = 0;
1.44 noro 1737: /* wc = A^(-1)wc; wc is not normalized */
1738: solve_by_lu_mod(w,rank,md,wc,ri,0);
1739: /* wx += q*wc */
1740: ptr = wx;
1741: for ( i = 0; i < rank; i++ )
1742: for ( j = 0, wi = wc[i]; j < ri; j++ ) {
1743: if ( wi[j] )
1744: muln_1(BD(NM(q)),PL(NM(q)),wi[j],ptr);
1745: ptr += nsize;
1746: }
1747: count++;
1.1 noro 1748: get_eg(&tmp1);
1.3 noro 1749: add_eg(&eg_inv,&tmp0,&tmp1);
1750: get_eg(&tmp0);
1751: for ( i = 0; i < rank; i++ )
1752: for ( j = 0; j < ri; j++ ) {
1753: inner_product_mat_int_mod(a,wc,rank,i,j,&u);
1754: addq(b[i][j],u,&s);
1755: if ( s ) {
1756: t = divin(NM(s),md,&tn);
1757: if ( t )
1758: error("generic_gauss_elim_hensel:incosistent");
1759: NTOQ(tn,SGN(s),b[i][j]);
1760: } else
1761: b[i][j] = 0;
1762: }
1763: get_eg(&tmp1);
1764: add_eg(&eg_mul,&tmp0,&tmp1);
1765: /* q = q*md */
1766: mulq(q,mdq,&u); q = u;
1.44 noro 1767: if ( count == period ) {
1.41 noro 1768: get_eg(&tmp0);
1.44 noro 1769: ptr = wx;
1770: for ( i = 0; i < rank; i++ )
1771: for ( j = 0, xi = x[i]; j < ri;
1772: j++, ptr += nsize ) {
1773: for ( k = nsize-1; k >= 0 && !ptr[k]; k-- );
1774: if ( k >= 0 ) {
1775: wn = NALLOC(k+1);
1776: PL(wn) = k+1;
1777: for ( l = 0; l <= k; l++ ) BD(wn)[l] = (unsigned int)ptr[l];
1778: NTOQ(wn,1,wq);
1779: subq(xi[j],wq,&u); xi[j] = u;
1780: }
1781: }
1.41 noro 1782: ret = intmtoratm_q(xmat,NM(q),*nmmat,dn);
1783: get_eg(&tmp1); add_eg(&eg_intrat,&tmp0,&tmp1);
1784: if ( ret ) {
1.50 noro 1785: rind = (int *)MALLOC_ATOMIC(rank*sizeof(int));
1786: cind = (int *)MALLOC_ATOMIC((col-rank)*sizeof(int));
1.39 noro 1787: for ( j = k = l = 0; j < col; j++ )
1788: if ( cinfo[j] )
1789: rind[k++] = j;
1790: else
1.50 noro 1791: cind[l++] = j;
1792: get_eg(&tmp0);
1793: ret = gensolve_check(mat,*nmmat,*dn,rind,cind);
1794: get_eg(&tmp1); add_eg(&eg_check,&tmp0,&tmp1);
1795: if ( ret ) {
1796: if ( DP_Print > 3 ) {
1797: fprintf(stderr,"\n");
1798: print_eg("INV",&eg_inv);
1799: print_eg("MUL",&eg_mul);
1800: print_eg("INTRAT",&eg_intrat);
1801: print_eg("CHECK",&eg_check);
1802: fflush(asir_out);
1803: }
1804: *rindp = rind;
1805: *cindp = cind;
1806: for ( j = k = 0; j < col; j++ )
1807: if ( !cinfo[j] )
1808: cind[k++] = j;
1809: return rank;
1810: }
1811: } else {
1812: period = period*3/2;
1813: count = 0;
1814: nsize += period;
1815: wxsize += rank*ri*nsize;
1816: wx = (int *)REALLOC(wx,wxsize*sizeof(int));
1817: for ( i = 0; i < wxsize; i++ ) wx[i] = 0;
1818: }
1819: }
1820: }
1821: }
1822: }
1823:
1.55 noro 1824: int generic_gauss_elim_hensel_dalg(MAT mat,DP *mb,MAT *nmmat,Q *dn,int **rindp,int **cindp)
1.50 noro 1825: {
1826: MAT bmat,xmat;
1827: Q **a0,**a,**b,**x,**nm;
1828: Q *ai,*bi,*xi;
1829: int row,col;
1830: int **w;
1831: int *wi;
1832: int **wc;
1833: Q mdq,q,s,u;
1834: N tn;
1835: int ind,md,i,j,k,l,li,ri,rank;
1836: unsigned int t;
1837: int *cinfo,*rinfo;
1838: int *rind,*cind;
1839: int count;
1840: int ret;
1841: struct oEGT eg_mul,eg_inv,eg_intrat,eg_check,tmp0,tmp1;
1842: int period;
1843: int *wx,*ptr;
1844: int wxsize,nsize;
1845: N wn;
1846: Q wq;
1847: NumberField nf;
1848: DP m;
1849: int col1;
1850:
1851: a0 = (Q **)mat->body;
1852: row = mat->row; col = mat->col;
1853: w = (int **)almat(row,col);
1854: for ( ind = 0; ; ind++ ) {
1855: md = get_lprime(ind);
1856: STOQ(md,mdq);
1857: for ( i = 0; i < row; i++ )
1858: for ( j = 0, ai = a0[i], wi = w[i]; j < col; j++ )
1859: if ( q = (Q)ai[j] ) {
1860: t = rem(NM(q),md);
1861: if ( t && SGN(q) < 0 )
1862: t = (md - t) % md;
1863: wi[j] = t;
1864: } else
1865: wi[j] = 0;
1866:
1867: if ( DP_Print ) {
1868: fprintf(asir_out,"LU decomposition.."); fflush(asir_out);
1869: }
1870: rank = find_lhs_and_lu_mod((unsigned int **)w,row,col,md,&rinfo,&cinfo);
1871: if ( DP_Print ) {
1872: fprintf(asir_out,"done.\n"); fflush(asir_out);
1873: }
1874: for ( i = 0; i < col-1; i++ ) {
1875: if ( !cinfo[i] ) {
1876: m = mb[i];
1877: for ( j = i+1; j < col-1; j++ )
1878: if ( dp_redble(mb[j],m) )
1879: cinfo[j] = -1;
1880: }
1881: }
1882: a = (Q **)almat_pointer(rank,rank); /* lhs mat */
1883: MKMAT(bmat,rank,col-rank); b = (Q **)bmat->body; /* lhs mat */
1884: for ( j = li = ri = 0; j < col; j++ )
1885: if ( cinfo[j] > 0 ) {
1886: /* the column is in lhs */
1887: for ( i = 0; i < rank; i++ ) {
1888: w[i][li] = w[i][j];
1889: a[i][li] = a0[rinfo[i]][j];
1890: }
1891: li++;
1892: } else if ( !cinfo[j] ) {
1893: /* the column is in rhs */
1894: for ( i = 0; i < rank; i++ )
1895: b[i][ri] = a0[rinfo[i]][j];
1896: ri++;
1897: }
1898:
1899: /* solve Ax+B=0; A: rank x rank, B: rank x ri */
1900: MKMAT(xmat,rank,ri); x = (Q **)(xmat)->body;
1901: MKMAT(*nmmat,rank,ri); nm = (Q **)(*nmmat)->body;
1902: /* use the right part of w as work area */
1903: wc = (int **)almat(rank,ri);
1904: for ( i = 0; i < rank; i++ )
1905: wc[i] = w[i]+rank;
1906: *rindp = rind = (int *)MALLOC_ATOMIC(rank*sizeof(int));
1907: *cindp = cind = (int *)MALLOC_ATOMIC((col-rank)*sizeof(int));
1908: init_eg(&eg_mul); init_eg(&eg_inv);
1909: init_eg(&eg_check); init_eg(&eg_intrat);
1910: period = F4_INTRAT_PERIOD;
1911: nsize = period;
1912: wxsize = rank*ri*nsize;
1913: wx = (int *)MALLOC_ATOMIC(wxsize*sizeof(int));
1914: for ( i = 0; i < wxsize; i++ ) wx[i] = 0;
1915: for ( q = ONE, count = 0; ; ) {
1916: if ( DP_Print )
1917: fprintf(stderr,"o");
1918: /* wc = -b mod md */
1919: get_eg(&tmp0);
1920: for ( i = 0; i < rank; i++ )
1921: for ( j = 0, bi = b[i], wi = wc[i]; j < ri; j++ )
1922: if ( u = (Q)bi[j] ) {
1923: t = rem(NM(u),md);
1924: if ( t && SGN(u) > 0 )
1925: t = (md - t) % md;
1926: wi[j] = t;
1927: } else
1928: wi[j] = 0;
1929: /* wc = A^(-1)wc; wc is not normalized */
1930: solve_by_lu_mod(w,rank,md,wc,ri,0);
1931: /* wx += q*wc */
1932: ptr = wx;
1933: for ( i = 0; i < rank; i++ )
1934: for ( j = 0, wi = wc[i]; j < ri; j++ ) {
1935: if ( wi[j] )
1936: muln_1(BD(NM(q)),PL(NM(q)),wi[j],ptr);
1937: ptr += nsize;
1938: }
1939: count++;
1940: get_eg(&tmp1);
1941: add_eg(&eg_inv,&tmp0,&tmp1);
1942: get_eg(&tmp0);
1943: for ( i = 0; i < rank; i++ )
1944: for ( j = 0; j < ri; j++ ) {
1945: inner_product_mat_int_mod(a,wc,rank,i,j,&u);
1946: addq(b[i][j],u,&s);
1947: if ( s ) {
1948: t = divin(NM(s),md,&tn);
1949: if ( t )
1950: error("generic_gauss_elim_hensel:incosistent");
1951: NTOQ(tn,SGN(s),b[i][j]);
1952: } else
1953: b[i][j] = 0;
1954: }
1955: get_eg(&tmp1);
1956: add_eg(&eg_mul,&tmp0,&tmp1);
1957: /* q = q*md */
1958: mulq(q,mdq,&u); q = u;
1959: if ( count == period ) {
1960: get_eg(&tmp0);
1961: ptr = wx;
1962: for ( i = 0; i < rank; i++ )
1963: for ( j = 0, xi = x[i]; j < ri;
1964: j++, ptr += nsize ) {
1965: for ( k = nsize-1; k >= 0 && !ptr[k]; k-- );
1966: if ( k >= 0 ) {
1967: wn = NALLOC(k+1);
1968: PL(wn) = k+1;
1969: for ( l = 0; l <= k; l++ ) BD(wn)[l] = (unsigned int)ptr[l];
1970: NTOQ(wn,1,wq);
1971: subq(xi[j],wq,&u); xi[j] = u;
1972: }
1973: }
1974: ret = intmtoratm_q(xmat,NM(q),*nmmat,dn);
1975: get_eg(&tmp1); add_eg(&eg_intrat,&tmp0,&tmp1);
1976: if ( ret ) {
1977: for ( j = k = l = 0; j < col; j++ )
1978: if ( cinfo[j] > 0 )
1979: rind[k++] = j;
1980: else if ( !cinfo[j] )
1.39 noro 1981: cind[l++] = j;
1.41 noro 1982: get_eg(&tmp0);
1983: ret = gensolve_check(mat,*nmmat,*dn,rind,cind);
1984: get_eg(&tmp1); add_eg(&eg_check,&tmp0,&tmp1);
1985: if ( ret ) {
1.42 noro 1986: if ( DP_Print > 3 ) {
1.40 noro 1987: fprintf(stderr,"\n");
1988: print_eg("INV",&eg_inv);
1989: print_eg("MUL",&eg_mul);
1.41 noro 1990: print_eg("INTRAT",&eg_intrat);
1991: print_eg("CHECK",&eg_check);
1.40 noro 1992: fflush(asir_out);
1993: }
1.39 noro 1994: return rank;
1995: }
1.44 noro 1996: } else {
1997: period = period*3/2;
1998: count = 0;
1999: nsize += period;
2000: wxsize += rank*ri*nsize;
2001: wx = (int *)REALLOC(wx,wxsize*sizeof(int));
2002: for ( i = 0; i < wxsize; i++ ) wx[i] = 0;
2003: }
1.41 noro 2004: }
1.1 noro 2005: }
2006: }
2007: }
2008:
2009: int f4_nocheck;
2010:
1.24 noro 2011: int gensolve_check(MAT mat,MAT nm,Q dn,int *rind,int *cind)
1.1 noro 2012: {
2013: int row,col,rank,clen,i,j,k,l;
1.24 noro 2014: Q s,t;
1.1 noro 2015: Q *w;
2016: Q *mati,*nmk;
2017:
2018: if ( f4_nocheck )
2019: return 1;
2020: row = mat->row; col = mat->col;
2021: rank = nm->row; clen = nm->col;
2022: w = (Q *)MALLOC(clen*sizeof(Q));
2023: for ( i = 0; i < row; i++ ) {
2024: mati = (Q *)mat->body[i];
2025: #if 1
2026: bzero(w,clen*sizeof(Q));
2027: for ( k = 0; k < rank; k++ )
2028: for ( l = 0, nmk = (Q *)nm->body[k]; l < clen; l++ ) {
2029: mulq(mati[rind[k]],nmk[l],&t);
2030: addq(w[l],t,&s); w[l] = s;
2031: }
2032: for ( j = 0; j < clen; j++ ) {
2033: mulq(dn,mati[cind[j]],&t);
2034: if ( cmpq(w[j],t) )
2035: break;
2036: }
2037: #else
2038: for ( j = 0; j < clen; j++ ) {
2039: for ( k = 0, s = 0; k < rank; k++ ) {
2040: mulq(mati[rind[k]],nm->body[k][j],&t);
2041: addq(s,t,&u); s = u;
2042: }
2043: mulq(dn,mati[cind[j]],&t);
2044: if ( cmpq(s,t) )
2045: break;
2046: }
2047: #endif
2048: if ( j != clen )
2049: break;
2050: }
2051: if ( i != row )
2052: return 0;
2053: else
2054: return 1;
2055: }
2056:
2057: /* assuming 0 < c < m */
2058:
1.24 noro 2059: int inttorat(N c,N m,N b,int *sgnp,N *nmp,N *dnp)
1.1 noro 2060: {
1.24 noro 2061: Q qq,t,u1,v1,r1;
2062: N q,u2,v2,r2;
1.1 noro 2063:
2064: u1 = 0; v1 = ONE; u2 = m; v2 = c;
2065: while ( cmpn(v2,b) >= 0 ) {
2066: divn(u2,v2,&q,&r2); u2 = v2; v2 = r2;
2067: NTOQ(q,1,qq); mulq(qq,v1,&t); subq(u1,t,&r1); u1 = v1; v1 = r1;
2068: }
2069: if ( cmpn(NM(v1),b) >= 0 )
2070: return 0;
2071: else {
2072: *nmp = v2;
2073: *dnp = NM(v1);
2074: *sgnp = SGN(v1);
2075: return 1;
2076: }
2077: }
2078:
2079: /* mat->body = N ** */
2080:
1.24 noro 2081: int intmtoratm(MAT mat,N md,MAT nm,Q *dn)
1.1 noro 2082: {
2083: N t,s,b;
1.24 noro 2084: Q dn0,dn1,nm1,q;
1.1 noro 2085: int i,j,k,l,row,col;
2086: Q **rmat;
2087: N **tmat;
2088: N *tmi;
2089: Q *nmk;
2090: N u,unm,udn;
2091: int sgn,ret;
2092:
1.3 noro 2093: if ( UNIN(md) )
2094: return 0;
1.1 noro 2095: row = mat->row; col = mat->col;
2096: bshiftn(md,1,&t);
2097: isqrt(t,&s);
2098: bshiftn(s,64,&b);
2099: if ( !b )
2100: b = ONEN;
2101: dn0 = ONE;
2102: tmat = (N **)mat->body;
2103: rmat = (Q **)nm->body;
2104: for ( i = 0; i < row; i++ )
2105: for ( j = 0, tmi = tmat[i]; j < col; j++ )
2106: if ( tmi[j] ) {
2107: muln(tmi[j],NM(dn0),&s);
2108: remn(s,md,&u);
2109: ret = inttorat(u,md,b,&sgn,&unm,&udn);
2110: if ( !ret )
2111: return 0;
2112: else {
2113: NTOQ(unm,sgn,nm1);
2114: NTOQ(udn,1,dn1);
2115: if ( !UNIQ(dn1) ) {
2116: for ( k = 0; k < i; k++ )
2117: for ( l = 0, nmk = rmat[k]; l < col; l++ ) {
2118: mulq(nmk[l],dn1,&q); nmk[l] = q;
2119: }
2120: for ( l = 0, nmk = rmat[i]; l < j; l++ ) {
2121: mulq(nmk[l],dn1,&q); nmk[l] = q;
2122: }
2123: }
2124: rmat[i][j] = nm1;
2125: mulq(dn0,dn1,&q); dn0 = q;
2126: }
2127: }
2128: *dn = dn0;
2129: return 1;
2130: }
2131:
1.3 noro 2132: /* mat->body = Q ** */
2133:
1.24 noro 2134: int intmtoratm_q(MAT mat,N md,MAT nm,Q *dn)
1.3 noro 2135: {
2136: N t,s,b;
1.24 noro 2137: Q dn0,dn1,nm1,q;
1.3 noro 2138: int i,j,k,l,row,col;
2139: Q **rmat;
2140: Q **tmat;
2141: Q *tmi;
2142: Q *nmk;
2143: N u,unm,udn;
2144: int sgn,ret;
2145:
2146: if ( UNIN(md) )
2147: return 0;
2148: row = mat->row; col = mat->col;
2149: bshiftn(md,1,&t);
2150: isqrt(t,&s);
2151: bshiftn(s,64,&b);
2152: if ( !b )
2153: b = ONEN;
2154: dn0 = ONE;
2155: tmat = (Q **)mat->body;
2156: rmat = (Q **)nm->body;
2157: for ( i = 0; i < row; i++ )
2158: for ( j = 0, tmi = tmat[i]; j < col; j++ )
2159: if ( tmi[j] ) {
2160: muln(NM(tmi[j]),NM(dn0),&s);
2161: remn(s,md,&u);
2162: ret = inttorat(u,md,b,&sgn,&unm,&udn);
2163: if ( !ret )
2164: return 0;
2165: else {
2166: if ( SGN(tmi[j])<0 )
2167: sgn = -sgn;
2168: NTOQ(unm,sgn,nm1);
2169: NTOQ(udn,1,dn1);
2170: if ( !UNIQ(dn1) ) {
2171: for ( k = 0; k < i; k++ )
2172: for ( l = 0, nmk = rmat[k]; l < col; l++ ) {
2173: mulq(nmk[l],dn1,&q); nmk[l] = q;
2174: }
2175: for ( l = 0, nmk = rmat[i]; l < j; l++ ) {
2176: mulq(nmk[l],dn1,&q); nmk[l] = q;
2177: }
2178: }
2179: rmat[i][j] = nm1;
2180: mulq(dn0,dn1,&q); dn0 = q;
2181: }
2182: }
2183: *dn = dn0;
2184: return 1;
2185: }
2186:
1.4 noro 2187: #define ONE_STEP1 if ( zzz = *s ) { DMAR(zzz,hc,*tj,md,*tj) } tj++; s++;
2188:
1.24 noro 2189: void reduce_reducers_mod(int **mat,int row,int col,int md)
1.4 noro 2190: {
2191: int i,j,k,l,hc,zzz;
2192: int *t,*s,*tj,*ind;
2193:
2194: /* reduce the reducers */
2195: ind = (int *)ALLOCA(row*sizeof(int));
2196: for ( i = 0; i < row; i++ ) {
2197: t = mat[i];
2198: for ( j = 0; j < col && !t[j]; j++ );
2199: /* register the position of the head term */
2200: ind[i] = j;
2201: for ( l = i-1; l >= 0; l-- ) {
2202: /* reduce mat[i] by mat[l] */
2203: if ( hc = t[ind[l]] ) {
2204: /* mat[i] = mat[i]-hc*mat[l] */
2205: j = ind[l];
2206: s = mat[l]+j;
2207: tj = t+j;
2208: hc = md-hc;
2209: k = col-j;
2210: for ( ; k >= 64; k -= 64 ) {
2211: ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
2212: ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
2213: ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
2214: ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
2215: ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
2216: ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
2217: ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
2218: ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
2219: ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
2220: ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
2221: ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
2222: ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
2223: ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
2224: ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
2225: ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
2226: ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
2227: }
1.16 noro 2228: for ( ; k > 0; k-- ) {
1.4 noro 2229: if ( zzz = *s ) { DMAR(zzz,hc,*tj,md,*tj) } tj++; s++;
2230: }
2231: }
2232: }
2233: }
2234: }
2235:
2236: /*
2237: mat[i] : reducers (i=0,...,nred-1)
2238: spolys (i=nred,...,row-1)
2239: mat[0] < mat[1] < ... < mat[nred-1] w.r.t the term order
2240: 1. reduce the reducers
2241: 2. reduce spolys by the reduced reducers
2242: */
2243:
1.24 noro 2244: void pre_reduce_mod(int **mat,int row,int col,int nred,int md)
1.4 noro 2245: {
2246: int i,j,k,l,hc,inv;
2247: int *t,*s,*tk,*ind;
2248:
2249: #if 1
2250: /* reduce the reducers */
2251: ind = (int *)ALLOCA(row*sizeof(int));
2252: for ( i = 0; i < nred; i++ ) {
2253: /* make mat[i] monic and mat[i] by mat[0],...,mat[i-1] */
2254: t = mat[i];
2255: for ( j = 0; j < col && !t[j]; j++ );
2256: /* register the position of the head term */
2257: ind[i] = j;
2258: inv = invm(t[j],md);
2259: for ( k = j; k < col; k++ )
2260: if ( t[k] )
2261: DMAR(t[k],inv,0,md,t[k])
2262: for ( l = i-1; l >= 0; l-- ) {
2263: /* reduce mat[i] by mat[l] */
2264: if ( hc = t[ind[l]] ) {
2265: /* mat[i] = mat[i]-hc*mat[l] */
2266: for ( k = ind[l], hc = md-hc, s = mat[l]+k, tk = t+k;
2267: k < col; k++, tk++, s++ )
2268: if ( *s )
2269: DMAR(*s,hc,*tk,md,*tk)
2270: }
2271: }
2272: }
2273: /* reduce the spolys */
2274: for ( i = nred; i < row; i++ ) {
2275: t = mat[i];
2276: for ( l = nred-1; l >= 0; l-- ) {
2277: /* reduce mat[i] by mat[l] */
2278: if ( hc = t[ind[l]] ) {
2279: /* mat[i] = mat[i]-hc*mat[l] */
2280: for ( k = ind[l], hc = md-hc, s = mat[l]+k, tk = t+k;
2281: k < col; k++, tk++, s++ )
2282: if ( *s )
2283: DMAR(*s,hc,*tk,md,*tk)
2284: }
2285: }
2286: }
2287: #endif
2288: }
2289: /*
2290: mat[i] : reducers (i=0,...,nred-1)
2291: mat[0] < mat[1] < ... < mat[nred-1] w.r.t the term order
2292: */
2293:
1.24 noro 2294: void reduce_sp_by_red_mod(int *sp,int **redmat,int *ind,int nred,int col,int md)
1.4 noro 2295: {
2296: int i,j,k,hc,zzz;
1.24 noro 2297: int *s,*tj;
1.4 noro 2298:
2299: /* reduce the spolys by redmat */
2300: for ( i = nred-1; i >= 0; i-- ) {
2301: /* reduce sp by redmat[i] */
2302: if ( hc = sp[ind[i]] ) {
2303: /* sp = sp-hc*redmat[i] */
2304: j = ind[i];
2305: hc = md-hc;
2306: s = redmat[i]+j;
2307: tj = sp+j;
1.16 noro 2308: for ( k = col-j; k > 0; k-- ) {
1.4 noro 2309: if ( zzz = *s ) { DMAR(zzz,hc,*tj,md,*tj) } tj++; s++;
1.15 noro 2310: }
2311: }
1.17 noro 2312: }
2313: }
2314:
2315: /*
1.15 noro 2316: mat[i] : compressed reducers (i=0,...,nred-1)
2317: mat[0] < mat[1] < ... < mat[nred-1] w.r.t the term order
2318: */
2319:
1.24 noro 2320: void red_by_compress(int m,unsigned int *p,unsigned int *r,
2321: unsigned int *ri,unsigned int hc,int len)
1.18 noro 2322: {
1.19 noro 2323: unsigned int up,lo;
1.18 noro 2324: unsigned int dmy;
2325: unsigned int *pj;
2326:
1.21 noro 2327: p[*ri] = 0; r++; ri++;
2328: for ( len--; len; len--, r++, ri++ ) {
2329: pj = p+ *ri;
2330: DMA(*r,hc,*pj,up,lo);
1.18 noro 2331: if ( up ) {
2332: DSAB(m,up,lo,dmy,*pj);
2333: } else
2334: *pj = lo;
2335: }
2336: }
2337:
2338: /* p -= hc*r */
2339:
1.24 noro 2340: void red_by_vect(int m,unsigned int *p,unsigned int *r,unsigned int hc,int len)
1.18 noro 2341: {
1.58 ohara 2342: unsigned int up,lo,dmy;
1.18 noro 2343:
2344: *p++ = 0; r++; len--;
2345: for ( ; len; len--, r++, p++ )
2346: if ( *r ) {
1.20 noro 2347: DMA(*r,hc,*p,up,lo);
1.18 noro 2348: if ( up ) {
2349: DSAB(m,up,lo,dmy,*p);
2350: } else
2351: *p = lo;
2352: }
2353: }
2354:
1.32 noro 2355: void red_by_vect_sf(int m,unsigned int *p,unsigned int *r,unsigned int hc,int len)
2356: {
2357: *p++ = 0; r++; len--;
2358: for ( ; len; len--, r++, p++ )
2359: if ( *r )
2360: *p = _addsf(_mulsf(*r,hc),*p);
2361: }
2362:
1.71 ! noro 2363: extern GZ current_mod_lf;
! 2364: extern int current_mod_lf_size;
! 2365:
1.70 noro 2366: void red_by_vect_lf(mpz_t *p,mpz_t *r,mpz_t hc,int len)
2367: {
2368: mpz_set_ui(*p++,0); r++; len--;
1.71 ! noro 2369: for ( ; len; len--, r++, p++ ) {
1.70 noro 2370: mpz_addmul(*p,*r,hc);
1.71 ! noro 2371: #if 0
! 2372: if ( mpz_size(*p) > current_mod_lf_size )
! 2373: mpz_mod(*p,*p,BDY(current_mod_lf));
! 2374: #endif
! 2375: }
1.70 noro 2376: }
2377:
2378:
1.21 noro 2379: extern unsigned int **psca;
2380:
1.24 noro 2381: void reduce_sp_by_red_mod_compress (int *sp,CDP *redmat,int *ind,
2382: int nred,int col,int md)
1.15 noro 2383: {
1.24 noro 2384: int i,len;
1.15 noro 2385: CDP ri;
1.24 noro 2386: unsigned int hc;
1.18 noro 2387: unsigned int *usp;
1.15 noro 2388:
1.18 noro 2389: usp = (unsigned int *)sp;
1.15 noro 2390: /* reduce the spolys by redmat */
2391: for ( i = nred-1; i >= 0; i-- ) {
2392: /* reduce sp by redmat[i] */
1.18 noro 2393: usp[ind[i]] %= md;
2394: if ( hc = usp[ind[i]] ) {
1.15 noro 2395: /* sp = sp-hc*redmat[i] */
2396: hc = md-hc;
2397: ri = redmat[i];
2398: len = ri->len;
1.21 noro 2399: red_by_compress(md,usp,psca[ri->psindex],ri->body,hc,len);
1.4 noro 2400: }
2401: }
1.18 noro 2402: for ( i = 0; i < col; i++ )
1.24 noro 2403: if ( usp[i] >= (unsigned int)md )
1.18 noro 2404: usp[i] %= md;
1.4 noro 2405: }
2406:
2407: #define ONE_STEP2 if ( zzz = *pk ) { DMAR(zzz,a,*tk,md,*tk) } pk++; tk++;
2408:
1.24 noro 2409: int generic_gauss_elim_mod(int **mat0,int row,int col,int md,int *colstat)
1.1 noro 2410: {
1.24 noro 2411: int i,j,k,l,inv,a,rank;
2412: unsigned int *t,*pivot,*pk;
1.18 noro 2413: unsigned int **mat;
1.1 noro 2414:
1.18 noro 2415: mat = (unsigned int **)mat0;
1.1 noro 2416: for ( rank = 0, j = 0; j < col; j++ ) {
1.18 noro 2417: for ( i = rank; i < row; i++ )
2418: mat[i][j] %= md;
2419: for ( i = rank; i < row; i++ )
2420: if ( mat[i][j] )
2421: break;
1.1 noro 2422: if ( i == row ) {
2423: colstat[j] = 0;
2424: continue;
2425: } else
2426: colstat[j] = 1;
2427: if ( i != rank ) {
2428: t = mat[i]; mat[i] = mat[rank]; mat[rank] = t;
2429: }
2430: pivot = mat[rank];
2431: inv = invm(pivot[j],md);
1.4 noro 2432: for ( k = j, pk = pivot+k; k < col; k++, pk++ )
2433: if ( *pk ) {
1.24 noro 2434: if ( *pk >= (unsigned int)md )
1.18 noro 2435: *pk %= md;
1.4 noro 2436: DMAR(*pk,inv,0,md,*pk)
1.1 noro 2437: }
2438: for ( i = rank+1; i < row; i++ ) {
2439: t = mat[i];
1.18 noro 2440: if ( a = t[j] )
2441: red_by_vect(md,t+j,pivot+j,md-a,col-j);
1.1 noro 2442: }
2443: rank++;
2444: }
2445: for ( j = col-1, l = rank-1; j >= 0; j-- )
2446: if ( colstat[j] ) {
2447: pivot = mat[l];
2448: for ( i = 0; i < l; i++ ) {
2449: t = mat[i];
1.18 noro 2450: t[j] %= md;
2451: if ( a = t[j] )
2452: red_by_vect(md,t+j,pivot+j,md-a,col-j);
1.1 noro 2453: }
2454: l--;
1.18 noro 2455: }
2456: for ( j = 0, l = 0; l < rank; j++ )
2457: if ( colstat[j] ) {
2458: t = mat[l];
2459: for ( k = j; k < col; k++ )
1.24 noro 2460: if ( t[k] >= (unsigned int)md )
1.18 noro 2461: t[k] %= md;
2462: l++;
1.32 noro 2463: }
2464: return rank;
2465: }
2466:
1.65 noro 2467: int generic_gauss_elim_mod2(int **mat0,int row,int col,int md,int *colstat,int *rowstat)
2468: {
2469: int i,j,k,l,inv,a,rank;
2470: unsigned int *t,*pivot,*pk;
2471: unsigned int **mat;
2472:
2473: for ( i = 0; i < row; i++ ) rowstat[i] = i;
2474: mat = (unsigned int **)mat0;
2475: for ( rank = 0, j = 0; j < col; j++ ) {
2476: for ( i = rank; i < row; i++ )
2477: mat[i][j] %= md;
2478: for ( i = rank; i < row; i++ )
2479: if ( mat[i][j] )
2480: break;
2481: if ( i == row ) {
2482: colstat[j] = 0;
2483: continue;
2484: } else
2485: colstat[j] = 1;
2486: if ( i != rank ) {
2487: t = mat[i]; mat[i] = mat[rank]; mat[rank] = t;
2488: k = rowstat[i]; rowstat[i] = rowstat[rank]; rowstat[rank] = k;
2489: }
2490: pivot = mat[rank];
2491: inv = invm(pivot[j],md);
2492: for ( k = j, pk = pivot+k; k < col; k++, pk++ )
2493: if ( *pk ) {
2494: if ( *pk >= (unsigned int)md )
2495: *pk %= md;
2496: DMAR(*pk,inv,0,md,*pk)
2497: }
2498: for ( i = rank+1; i < row; i++ ) {
2499: t = mat[i];
2500: if ( a = t[j] )
2501: red_by_vect(md,t+j,pivot+j,md-a,col-j);
2502: }
2503: rank++;
2504: }
2505: for ( j = col-1, l = rank-1; j >= 0; j-- )
2506: if ( colstat[j] ) {
2507: pivot = mat[l];
2508: for ( i = 0; i < l; i++ ) {
2509: t = mat[i];
2510: t[j] %= md;
2511: if ( a = t[j] )
2512: red_by_vect(md,t+j,pivot+j,md-a,col-j);
2513: }
2514: l--;
2515: }
2516: for ( j = 0, l = 0; l < rank; j++ )
2517: if ( colstat[j] ) {
2518: t = mat[l];
2519: for ( k = j; k < col; k++ )
2520: if ( t[k] >= (unsigned int)md )
2521: t[k] %= md;
2522: l++;
2523: }
2524: return rank;
2525: }
2526:
1.69 noro 2527: int indep_rows_mod(int **mat0,int row,int col,int md,int *rowstat)
2528: {
2529: int i,j,k,l,inv,a,rank;
2530: unsigned int *t,*pivot,*pk;
2531: unsigned int **mat;
2532:
2533: for ( i = 0; i < row; i++ ) rowstat[i] = i;
2534: mat = (unsigned int **)mat0;
2535: for ( rank = 0, j = 0; j < col; j++ ) {
2536: for ( i = rank; i < row; i++ )
2537: mat[i][j] %= md;
2538: for ( i = rank; i < row; i++ )
2539: if ( mat[i][j] )
2540: break;
2541: if ( i == row ) continue;
2542: if ( i != rank ) {
2543: t = mat[i]; mat[i] = mat[rank]; mat[rank] = t;
2544: k = rowstat[i]; rowstat[i] = rowstat[rank]; rowstat[rank] = k;
2545: }
2546: pivot = mat[rank];
2547: inv = invm(pivot[j],md);
2548: for ( k = j, pk = pivot+k; k < col; k++, pk++ )
2549: if ( *pk ) {
2550: if ( *pk >= (unsigned int)md )
2551: *pk %= md;
2552: DMAR(*pk,inv,0,md,*pk)
2553: }
2554: for ( i = rank+1; i < row; i++ ) {
2555: t = mat[i];
2556: if ( a = t[j] )
2557: red_by_vect(md,t+j,pivot+j,md-a,col-j);
2558: }
2559: rank++;
2560: }
2561: return rank;
2562: }
2563:
1.32 noro 2564: int generic_gauss_elim_sf(int **mat0,int row,int col,int md,int *colstat)
2565: {
2566: int i,j,k,l,inv,a,rank;
2567: unsigned int *t,*pivot,*pk;
2568: unsigned int **mat;
2569:
2570: mat = (unsigned int **)mat0;
2571: for ( rank = 0, j = 0; j < col; j++ ) {
2572: for ( i = rank; i < row; i++ )
2573: if ( mat[i][j] )
2574: break;
2575: if ( i == row ) {
2576: colstat[j] = 0;
2577: continue;
2578: } else
2579: colstat[j] = 1;
2580: if ( i != rank ) {
2581: t = mat[i]; mat[i] = mat[rank]; mat[rank] = t;
2582: }
2583: pivot = mat[rank];
2584: inv = _invsf(pivot[j]);
2585: for ( k = j, pk = pivot+k; k < col; k++, pk++ )
2586: if ( *pk )
2587: *pk = _mulsf(*pk,inv);
2588: for ( i = rank+1; i < row; i++ ) {
2589: t = mat[i];
2590: if ( a = t[j] )
2591: red_by_vect_sf(md,t+j,pivot+j,_chsgnsf(a),col-j);
2592: }
2593: rank++;
2594: }
2595: for ( j = col-1, l = rank-1; j >= 0; j-- )
2596: if ( colstat[j] ) {
2597: pivot = mat[l];
2598: for ( i = 0; i < l; i++ ) {
2599: t = mat[i];
2600: if ( a = t[j] )
2601: red_by_vect_sf(md,t+j,pivot+j,_chsgnsf(a),col-j);
2602: }
2603: l--;
1.1 noro 2604: }
2605: return rank;
2606: }
2607:
2608: /* LU decomposition; a[i][i] = 1/U[i][i] */
2609:
1.24 noro 2610: int lu_gfmmat(GFMMAT mat,unsigned int md,int *perm)
1.1 noro 2611: {
2612: int row,col;
1.24 noro 2613: int i,j,k;
1.1 noro 2614: unsigned int *t,*pivot;
2615: unsigned int **a;
2616: unsigned int inv,m;
2617:
2618: row = mat->row; col = mat->col;
2619: a = mat->body;
2620: bzero(perm,row*sizeof(int));
2621:
2622: for ( i = 0; i < row; i++ )
2623: perm[i] = i;
2624: for ( k = 0; k < col; k++ ) {
2625: for ( i = k; i < row && !a[i][k]; i++ );
2626: if ( i == row )
2627: return 0;
2628: if ( i != k ) {
2629: j = perm[i]; perm[i] = perm[k]; perm[k] = j;
2630: t = a[i]; a[i] = a[k]; a[k] = t;
2631: }
2632: pivot = a[k];
2633: pivot[k] = inv = invm(pivot[k],md);
2634: for ( i = k+1; i < row; i++ ) {
2635: t = a[i];
2636: if ( m = t[k] ) {
2637: DMAR(inv,m,0,md,t[k])
2638: for ( j = k+1, m = md - t[k]; j < col; j++ )
2639: if ( pivot[j] ) {
1.8 noro 2640: unsigned int tj;
2641:
2642: DMAR(m,pivot[j],t[j],md,tj)
2643: t[j] = tj;
1.1 noro 2644: }
2645: }
2646: }
2647: }
2648: return 1;
2649: }
2650:
1.3 noro 2651: /*
2652: Input
2653: a: a row x col matrix
2654: md : a modulus
2655:
2656: Output:
2657: return : d = the rank of mat
2658: a[0..(d-1)][0..(d-1)] : LU decomposition (a[i][i] = 1/U[i][i])
2659: rinfo: array of length row
2660: cinfo: array of length col
2661: i-th row in new a <-> rinfo[i]-th row in old a
2662: cinfo[j]=1 <=> j-th column is contained in the LU decomp.
2663: */
2664:
1.24 noro 2665: int find_lhs_and_lu_mod(unsigned int **a,int row,int col,
2666: unsigned int md,int **rinfo,int **cinfo)
1.3 noro 2667: {
1.24 noro 2668: int i,j,k,d;
1.3 noro 2669: int *rp,*cp;
2670: unsigned int *t,*pivot;
2671: unsigned int inv,m;
2672:
2673: *rinfo = rp = (int *)MALLOC_ATOMIC(row*sizeof(int));
2674: *cinfo = cp = (int *)MALLOC_ATOMIC(col*sizeof(int));
2675: for ( i = 0; i < row; i++ )
2676: rp[i] = i;
2677: for ( k = 0, d = 0; k < col; k++ ) {
2678: for ( i = d; i < row && !a[i][k]; i++ );
2679: if ( i == row ) {
2680: cp[k] = 0;
2681: continue;
2682: } else
2683: cp[k] = 1;
2684: if ( i != d ) {
2685: j = rp[i]; rp[i] = rp[d]; rp[d] = j;
2686: t = a[i]; a[i] = a[d]; a[d] = t;
2687: }
2688: pivot = a[d];
2689: pivot[k] = inv = invm(pivot[k],md);
2690: for ( i = d+1; i < row; i++ ) {
2691: t = a[i];
2692: if ( m = t[k] ) {
2693: DMAR(inv,m,0,md,t[k])
2694: for ( j = k+1, m = md - t[k]; j < col; j++ )
2695: if ( pivot[j] ) {
1.8 noro 2696: unsigned int tj;
2697: DMAR(m,pivot[j],t[j],md,tj)
2698: t[j] = tj;
1.3 noro 2699: }
2700: }
2701: }
2702: d++;
2703: }
2704: return d;
2705: }
2706:
1.53 noro 2707: int lu_mod(unsigned int **a,int n,unsigned int md,int **rinfo)
2708: {
2709: int i,j,k;
2710: int *rp;
2711: unsigned int *t,*pivot;
2712: unsigned int inv,m;
2713:
2714: *rinfo = rp = (int *)MALLOC_ATOMIC(n*sizeof(int));
2715: for ( i = 0; i < n; i++ ) rp[i] = i;
2716: for ( k = 0; k < n; k++ ) {
2717: for ( i = k; i < n && !a[i][k]; i++ );
2718: if ( i == n ) return 0;
2719: if ( i != k ) {
2720: j = rp[i]; rp[i] = rp[k]; rp[k] = j;
2721: t = a[i]; a[i] = a[k]; a[k] = t;
2722: }
2723: pivot = a[k];
2724: inv = invm(pivot[k],md);
2725: for ( i = k+1; i < n; i++ ) {
2726: t = a[i];
2727: if ( m = t[k] ) {
2728: DMAR(inv,m,0,md,t[k])
2729: for ( j = k+1, m = md - t[k]; j < n; j++ )
2730: if ( pivot[j] ) {
2731: unsigned int tj;
2732: DMAR(m,pivot[j],t[j],md,tj)
2733: t[j] = tj;
2734: }
2735: }
2736: }
2737: }
2738: return 1;
2739: }
2740:
1.3 noro 2741: /*
2742: Input
2743: a : n x n matrix; a result of LU-decomposition
2744: md : modulus
2745: b : n x l matrix
2746: Output
2747: b = a^(-1)b
2748: */
2749:
1.44 noro 2750: void solve_by_lu_mod(int **a,int n,int md,int **b,int l,int normalize)
1.3 noro 2751: {
2752: unsigned int *y,*c;
2753: int i,j,k;
2754: unsigned int t,m,m2;
2755:
2756: y = (int *)MALLOC_ATOMIC(n*sizeof(int));
2757: c = (int *)MALLOC_ATOMIC(n*sizeof(int));
2758: m2 = md>>1;
2759: for ( k = 0; k < l; k++ ) {
2760: /* copy b[.][k] to c */
2761: for ( i = 0; i < n; i++ )
2762: c[i] = (unsigned int)b[i][k];
2763: /* solve Ly=c */
2764: for ( i = 0; i < n; i++ ) {
2765: for ( t = c[i], j = 0; j < i; j++ )
2766: if ( a[i][j] ) {
2767: m = md - a[i][j];
2768: DMAR(m,y[j],t,md,t)
2769: }
2770: y[i] = t;
2771: }
2772: /* solve Uc=y */
2773: for ( i = n-1; i >= 0; i-- ) {
2774: for ( t = y[i], j =i+1; j < n; j++ )
2775: if ( a[i][j] ) {
2776: m = md - a[i][j];
2777: DMAR(m,c[j],t,md,t)
2778: }
2779: /* a[i][i] = 1/U[i][i] */
2780: DMAR(t,a[i][i],0,md,c[i])
2781: }
2782: /* copy c to b[.][k] with normalization */
1.44 noro 2783: if ( normalize )
2784: for ( i = 0; i < n; i++ )
2785: b[i][k] = (int)(c[i]>m2 ? c[i]-md : c[i]);
2786: else
2787: for ( i = 0; i < n; i++ )
2788: b[i][k] = c[i];
1.3 noro 2789: }
2790: }
2791:
1.24 noro 2792: void Pleqm1(NODE arg,VECT *rp)
1.1 noro 2793: {
2794: MAT m;
2795: VECT vect;
2796: pointer **mat;
2797: Q *v;
2798: Q q;
2799: int **wmat;
2800: int md,i,j,row,col,t,n,status;
2801:
2802: asir_assert(ARG0(arg),O_MAT,"leqm1");
2803: asir_assert(ARG1(arg),O_N,"leqm1");
2804: m = (MAT)ARG0(arg); md = QTOS((Q)ARG1(arg));
2805: row = m->row; col = m->col; mat = m->body;
2806: wmat = (int **)almat(row,col);
2807: for ( i = 0; i < row; i++ )
2808: for ( j = 0; j < col; j++ )
2809: if ( q = (Q)mat[i][j] ) {
2810: t = rem(NM(q),md);
2811: if ( SGN(q) < 0 )
2812: t = (md - t) % md;
2813: wmat[i][j] = t;
2814: } else
2815: wmat[i][j] = 0;
2816: status = gauss_elim_mod1(wmat,row,col,md);
2817: if ( status < 0 )
2818: *rp = 0;
2819: else if ( status > 0 )
2820: *rp = (VECT)ONE;
2821: else {
2822: n = col - 1;
2823: MKVECT(vect,n);
2824: for ( i = 0, v = (Q *)vect->body; i < n; i++ ) {
2825: t = (md-wmat[i][n])%md; STOQ(t,v[i]);
2826: }
2827: *rp = vect;
2828: }
2829: }
2830:
1.24 noro 2831: int gauss_elim_mod1(int **mat,int row,int col,int md)
1.1 noro 2832: {
2833: int i,j,k,inv,a,n;
2834: int *t,*pivot;
2835:
2836: n = col - 1;
2837: for ( j = 0; j < n; j++ ) {
2838: for ( i = j; i < row && !mat[i][j]; i++ );
2839: if ( i == row )
2840: return 1;
2841: if ( i != j ) {
2842: t = mat[i]; mat[i] = mat[j]; mat[j] = t;
2843: }
2844: pivot = mat[j];
2845: inv = invm(pivot[j],md);
2846: for ( k = j; k <= n; k++ )
2847: pivot[k] = dmar(pivot[k],inv,0,md);
2848: for ( i = j+1; i < row; i++ ) {
2849: t = mat[i];
2850: if ( i != j && (a = t[j]) )
2851: for ( k = j, a = md - a; k <= n; k++ )
2852: t[k] = dmar(pivot[k],a,t[k],md);
2853: }
2854: }
2855: for ( i = n; i < row && !mat[i][n]; i++ );
2856: if ( i == row ) {
2857: for ( j = n-1; j >= 0; j-- ) {
2858: for ( i = j-1, a = (md-mat[j][n])%md; i >= 0; i-- ) {
2859: mat[i][n] = dmar(mat[i][j],a,mat[i][n],md);
2860: mat[i][j] = 0;
2861: }
2862: }
2863: return 0;
2864: } else
2865: return -1;
2866: }
2867:
1.24 noro 2868: void Pgeninvm(NODE arg,LIST *rp)
1.1 noro 2869: {
2870: MAT m;
2871: pointer **mat;
2872: Q **tmat;
2873: Q q;
2874: unsigned int **wmat;
2875: int md,i,j,row,col,t,status;
2876: MAT mat1,mat2;
2877: NODE node1,node2;
2878:
2879: asir_assert(ARG0(arg),O_MAT,"leqm1");
2880: asir_assert(ARG1(arg),O_N,"leqm1");
2881: m = (MAT)ARG0(arg); md = QTOS((Q)ARG1(arg));
2882: row = m->row; col = m->col; mat = m->body;
2883: wmat = (unsigned int **)almat(row,col+row);
2884: for ( i = 0; i < row; i++ ) {
2885: bzero((char *)wmat[i],(col+row)*sizeof(int));
2886: for ( j = 0; j < col; j++ )
2887: if ( q = (Q)mat[i][j] ) {
2888: t = rem(NM(q),md);
2889: if ( SGN(q) < 0 )
2890: t = (md - t) % md;
2891: wmat[i][j] = t;
2892: }
2893: wmat[i][col+i] = 1;
2894: }
2895: status = gauss_elim_geninv_mod(wmat,row,col,md);
2896: if ( status > 0 )
2897: *rp = 0;
2898: else {
2899: MKMAT(mat1,col,row); MKMAT(mat2,row-col,row);
2900: for ( i = 0, tmat = (Q **)mat1->body; i < col; i++ )
2901: for ( j = 0; j < row; j++ )
1.24 noro 2902: UTOQ(wmat[i][j+col],tmat[i][j]);
1.1 noro 2903: for ( tmat = (Q **)mat2->body; i < row; i++ )
2904: for ( j = 0; j < row; j++ )
1.24 noro 2905: UTOQ(wmat[i][j+col],tmat[i-col][j]);
1.1 noro 2906: MKNODE(node2,mat2,0); MKNODE(node1,mat1,node2); MKLIST(*rp,node1);
2907: }
2908: }
2909:
1.24 noro 2910: int gauss_elim_geninv_mod(unsigned int **mat,int row,int col,int md)
1.1 noro 2911: {
2912: int i,j,k,inv,a,n,m;
2913: unsigned int *t,*pivot;
2914:
2915: n = col; m = row+col;
2916: for ( j = 0; j < n; j++ ) {
2917: for ( i = j; i < row && !mat[i][j]; i++ );
2918: if ( i == row )
2919: return 1;
2920: if ( i != j ) {
2921: t = mat[i]; mat[i] = mat[j]; mat[j] = t;
2922: }
2923: pivot = mat[j];
2924: inv = invm(pivot[j],md);
2925: for ( k = j; k < m; k++ )
2926: pivot[k] = dmar(pivot[k],inv,0,md);
2927: for ( i = j+1; i < row; i++ ) {
2928: t = mat[i];
2929: if ( a = t[j] )
2930: for ( k = j, a = md - a; k < m; k++ )
2931: t[k] = dmar(pivot[k],a,t[k],md);
2932: }
2933: }
2934: for ( j = n-1; j >= 0; j-- ) {
2935: pivot = mat[j];
2936: for ( i = j-1; i >= 0; i-- ) {
2937: t = mat[i];
2938: if ( a = t[j] )
2939: for ( k = j, a = md - a; k < m; k++ )
2940: t[k] = dmar(pivot[k],a,t[k],md);
2941: }
2942: }
2943: return 0;
2944: }
2945:
1.24 noro 2946: void Psolve_by_lu_gfmmat(NODE arg,VECT *rp)
1.1 noro 2947: {
2948: GFMMAT lu;
2949: Q *perm,*rhs,*v;
2950: int n,i;
2951: unsigned int md;
2952: unsigned int *b,*sol;
2953: VECT r;
2954:
2955: lu = (GFMMAT)ARG0(arg);
2956: perm = (Q *)BDY((VECT)ARG1(arg));
2957: rhs = (Q *)BDY((VECT)ARG2(arg));
2958: md = (unsigned int)QTOS((Q)ARG3(arg));
2959: n = lu->col;
2960: b = (unsigned int *)MALLOC_ATOMIC(n*sizeof(int));
2961: sol = (unsigned int *)MALLOC_ATOMIC(n*sizeof(int));
2962: for ( i = 0; i < n; i++ )
2963: b[i] = QTOS(rhs[QTOS(perm[i])]);
2964: solve_by_lu_gfmmat(lu,md,b,sol);
2965: MKVECT(r,n);
2966: for ( i = 0, v = (Q *)r->body; i < n; i++ )
1.24 noro 2967: UTOQ(sol[i],v[i]);
1.1 noro 2968: *rp = r;
2969: }
2970:
1.24 noro 2971: void solve_by_lu_gfmmat(GFMMAT lu,unsigned int md,
2972: unsigned int *b,unsigned int *x)
1.1 noro 2973: {
2974: int n;
2975: unsigned int **a;
2976: unsigned int *y;
2977: int i,j;
2978: unsigned int t,m;
2979:
2980: n = lu->col;
2981: a = lu->body;
2982: y = (unsigned int *)MALLOC_ATOMIC(n*sizeof(int));
2983: /* solve Ly=b */
2984: for ( i = 0; i < n; i++ ) {
2985: for ( t = b[i], j = 0; j < i; j++ )
2986: if ( a[i][j] ) {
2987: m = md - a[i][j];
2988: DMAR(m,y[j],t,md,t)
2989: }
2990: y[i] = t;
2991: }
2992: /* solve Ux=y */
2993: for ( i = n-1; i >= 0; i-- ) {
2994: for ( t = y[i], j =i+1; j < n; j++ )
2995: if ( a[i][j] ) {
2996: m = md - a[i][j];
2997: DMAR(m,x[j],t,md,t)
2998: }
2999: /* a[i][i] = 1/U[i][i] */
3000: DMAR(t,a[i][i],0,md,x[i])
3001: }
3002: }
3003:
1.53 noro 3004: void Plu_mat(NODE arg,LIST *rp)
3005: {
3006: MAT m,lu;
3007: Q dn;
3008: Q *v;
3009: int n,i;
3010: int *iperm;
3011: VECT perm;
3012: NODE n0;
3013:
3014: asir_assert(ARG0(arg),O_MAT,"lu_mat");
3015: m = (MAT)ARG0(arg);
3016: n = m->row;
3017: MKMAT(lu,n,n);
3018: lu_dec_cr(m,lu,&dn,&iperm);
3019: MKVECT(perm,n);
3020: for ( i = 0, v = (Q *)perm->body; i < n; i++ )
3021: STOQ(iperm[i],v[i]);
3022: n0 = mknode(3,lu,dn,perm);
3023: MKLIST(*rp,n0);
3024: }
3025:
1.24 noro 3026: void Plu_gfmmat(NODE arg,LIST *rp)
1.1 noro 3027: {
3028: MAT m;
3029: GFMMAT mm;
3030: unsigned int md;
3031: int i,row,col,status;
3032: int *iperm;
3033: Q *v;
3034: VECT perm;
3035: NODE n0;
3036:
1.53 noro 3037: asir_assert(ARG0(arg),O_MAT,"lu_gfmmat");
3038: asir_assert(ARG1(arg),O_N,"lu_gfmmat");
1.1 noro 3039: m = (MAT)ARG0(arg); md = (unsigned int)QTOS((Q)ARG1(arg));
3040: mat_to_gfmmat(m,md,&mm);
3041: row = m->row;
3042: col = m->col;
3043: iperm = (int *)MALLOC_ATOMIC(row*sizeof(int));
3044: status = lu_gfmmat(mm,md,iperm);
3045: if ( !status )
3046: n0 = 0;
3047: else {
3048: MKVECT(perm,row);
3049: for ( i = 0, v = (Q *)perm->body; i < row; i++ )
3050: STOQ(iperm[i],v[i]);
3051: n0 = mknode(2,mm,perm);
3052: }
3053: MKLIST(*rp,n0);
3054: }
3055:
1.24 noro 3056: void Pmat_to_gfmmat(NODE arg,GFMMAT *rp)
1.1 noro 3057: {
3058: MAT m;
3059: unsigned int md;
3060:
3061: asir_assert(ARG0(arg),O_MAT,"mat_to_gfmmat");
3062: asir_assert(ARG1(arg),O_N,"mat_to_gfmmat");
3063: m = (MAT)ARG0(arg); md = (unsigned int)QTOS((Q)ARG1(arg));
3064: mat_to_gfmmat(m,md,rp);
3065: }
3066:
1.24 noro 3067: void mat_to_gfmmat(MAT m,unsigned int md,GFMMAT *rp)
1.1 noro 3068: {
3069: unsigned int **wmat;
3070: unsigned int t;
3071: Q **mat;
3072: Q q;
3073: int i,j,row,col;
3074:
3075: row = m->row; col = m->col; mat = (Q **)m->body;
3076: wmat = (unsigned int **)almat(row,col);
3077: for ( i = 0; i < row; i++ ) {
3078: bzero((char *)wmat[i],col*sizeof(unsigned int));
3079: for ( j = 0; j < col; j++ )
3080: if ( q = mat[i][j] ) {
3081: t = (unsigned int)rem(NM(q),md);
3082: if ( SGN(q) < 0 )
3083: t = (md - t) % md;
3084: wmat[i][j] = t;
3085: }
3086: }
3087: TOGFMMAT(row,col,wmat,*rp);
3088: }
3089:
1.27 noro 3090: void Pgeninvm_swap(arg,rp)
3091: NODE arg;
3092: LIST *rp;
1.1 noro 3093: {
3094: MAT m;
3095: pointer **mat;
3096: Q **tmat;
3097: Q *tvect;
3098: Q q;
3099: unsigned int **wmat,**invmat;
3100: int *index;
3101: unsigned int t,md;
3102: int i,j,row,col,status;
3103: MAT mat1;
3104: VECT vect1;
3105: NODE node1,node2;
3106:
3107: asir_assert(ARG0(arg),O_MAT,"geninvm_swap");
3108: asir_assert(ARG1(arg),O_N,"geninvm_swap");
3109: m = (MAT)ARG0(arg); md = QTOS((Q)ARG1(arg));
3110: row = m->row; col = m->col; mat = m->body;
3111: wmat = (unsigned int **)almat(row,col+row);
3112: for ( i = 0; i < row; i++ ) {
3113: bzero((char *)wmat[i],(col+row)*sizeof(int));
3114: for ( j = 0; j < col; j++ )
3115: if ( q = (Q)mat[i][j] ) {
3116: t = (unsigned int)rem(NM(q),md);
3117: if ( SGN(q) < 0 )
3118: t = (md - t) % md;
3119: wmat[i][j] = t;
3120: }
3121: wmat[i][col+i] = 1;
3122: }
3123: status = gauss_elim_geninv_mod_swap(wmat,row,col,md,&invmat,&index);
3124: if ( status > 0 )
3125: *rp = 0;
3126: else {
3127: MKMAT(mat1,col,col);
3128: for ( i = 0, tmat = (Q **)mat1->body; i < col; i++ )
3129: for ( j = 0; j < col; j++ )
3130: UTOQ(invmat[i][j],tmat[i][j]);
3131: MKVECT(vect1,row);
3132: for ( i = 0, tvect = (Q *)vect1->body; i < row; i++ )
3133: STOQ(index[i],tvect[i]);
3134: MKNODE(node2,vect1,0); MKNODE(node1,mat1,node2); MKLIST(*rp,node1);
3135: }
3136: }
3137:
1.27 noro 3138: gauss_elim_geninv_mod_swap(mat,row,col,md,invmatp,indexp)
3139: unsigned int **mat;
3140: int row,col;
3141: unsigned int md;
3142: unsigned int ***invmatp;
3143: int **indexp;
1.1 noro 3144: {
3145: int i,j,k,inv,a,n,m;
3146: unsigned int *t,*pivot,*s;
3147: int *index;
3148: unsigned int **invmat;
3149:
3150: n = col; m = row+col;
3151: *indexp = index = (int *)MALLOC_ATOMIC(row*sizeof(int));
3152: for ( i = 0; i < row; i++ )
3153: index[i] = i;
3154: for ( j = 0; j < n; j++ ) {
3155: for ( i = j; i < row && !mat[i][j]; i++ );
3156: if ( i == row ) {
3157: *indexp = 0; *invmatp = 0; return 1;
3158: }
3159: if ( i != j ) {
3160: t = mat[i]; mat[i] = mat[j]; mat[j] = t;
3161: k = index[i]; index[i] = index[j]; index[j] = k;
3162: }
3163: pivot = mat[j];
3164: inv = (unsigned int)invm(pivot[j],md);
3165: for ( k = j; k < m; k++ )
3166: if ( pivot[k] )
3167: pivot[k] = (unsigned int)dmar(pivot[k],inv,0,md);
3168: for ( i = j+1; i < row; i++ ) {
3169: t = mat[i];
3170: if ( a = t[j] )
3171: for ( k = j, a = md - a; k < m; k++ )
3172: if ( pivot[k] )
3173: t[k] = dmar(pivot[k],a,t[k],md);
3174: }
3175: }
3176: for ( j = n-1; j >= 0; j-- ) {
3177: pivot = mat[j];
3178: for ( i = j-1; i >= 0; i-- ) {
3179: t = mat[i];
3180: if ( a = t[j] )
3181: for ( k = j, a = md - a; k < m; k++ )
3182: if ( pivot[k] )
3183: t[k] = dmar(pivot[k],a,t[k],md);
3184: }
3185: }
3186: *invmatp = invmat = (unsigned int **)almat(col,col);
1.27 noro 3187: for ( i = 0; i < col; i++ )
3188: for ( j = 0, s = invmat[i], t = mat[i]; j < col; j++ )
3189: s[j] = t[col+index[j]];
3190: return 0;
3191: }
3192:
3193: void Pgeninv_sf_swap(NODE arg,LIST *rp)
3194: {
3195: MAT m;
3196: GFS **mat,**tmat;
3197: Q *tvect;
3198: GFS q;
3199: int **wmat,**invmat;
3200: int *index;
3201: unsigned int t;
3202: int i,j,row,col,status;
3203: MAT mat1;
3204: VECT vect1;
3205: NODE node1,node2;
3206:
3207: asir_assert(ARG0(arg),O_MAT,"geninv_sf_swap");
3208: m = (MAT)ARG0(arg);
3209: row = m->row; col = m->col; mat = (GFS **)m->body;
3210: wmat = (int **)almat(row,col+row);
3211: for ( i = 0; i < row; i++ ) {
3212: bzero((char *)wmat[i],(col+row)*sizeof(int));
3213: for ( j = 0; j < col; j++ )
3214: if ( q = (GFS)mat[i][j] )
3215: wmat[i][j] = FTOIF(CONT(q));
3216: wmat[i][col+i] = _onesf();
3217: }
3218: status = gauss_elim_geninv_sf_swap(wmat,row,col,&invmat,&index);
3219: if ( status > 0 )
3220: *rp = 0;
3221: else {
3222: MKMAT(mat1,col,col);
3223: for ( i = 0, tmat = (GFS **)mat1->body; i < col; i++ )
3224: for ( j = 0; j < col; j++ )
3225: if ( t = invmat[i][j] ) {
3226: MKGFS(IFTOF(t),tmat[i][j]);
3227: }
3228: MKVECT(vect1,row);
3229: for ( i = 0, tvect = (Q *)vect1->body; i < row; i++ )
3230: STOQ(index[i],tvect[i]);
3231: MKNODE(node2,vect1,0); MKNODE(node1,mat1,node2); MKLIST(*rp,node1);
3232: }
3233: }
3234:
3235: int gauss_elim_geninv_sf_swap(int **mat,int row,int col,
3236: int ***invmatp,int **indexp)
3237: {
3238: int i,j,k,inv,a,n,m,u;
3239: int *t,*pivot,*s;
3240: int *index;
3241: int **invmat;
3242:
3243: n = col; m = row+col;
3244: *indexp = index = (int *)MALLOC_ATOMIC(row*sizeof(int));
3245: for ( i = 0; i < row; i++ )
3246: index[i] = i;
3247: for ( j = 0; j < n; j++ ) {
3248: for ( i = j; i < row && !mat[i][j]; i++ );
3249: if ( i == row ) {
3250: *indexp = 0; *invmatp = 0; return 1;
3251: }
3252: if ( i != j ) {
3253: t = mat[i]; mat[i] = mat[j]; mat[j] = t;
3254: k = index[i]; index[i] = index[j]; index[j] = k;
3255: }
3256: pivot = mat[j];
3257: inv = _invsf(pivot[j]);
3258: for ( k = j; k < m; k++ )
3259: if ( pivot[k] )
3260: pivot[k] = _mulsf(pivot[k],inv);
3261: for ( i = j+1; i < row; i++ ) {
3262: t = mat[i];
3263: if ( a = t[j] )
3264: for ( k = j, a = _chsgnsf(a); k < m; k++ )
3265: if ( pivot[k] ) {
3266: u = _mulsf(pivot[k],a);
3267: t[k] = _addsf(u,t[k]);
3268: }
3269: }
3270: }
3271: for ( j = n-1; j >= 0; j-- ) {
3272: pivot = mat[j];
3273: for ( i = j-1; i >= 0; i-- ) {
3274: t = mat[i];
3275: if ( a = t[j] )
3276: for ( k = j, a = _chsgnsf(a); k < m; k++ )
3277: if ( pivot[k] ) {
3278: u = _mulsf(pivot[k],a);
3279: t[k] = _addsf(u,t[k]);
3280: }
3281: }
3282: }
3283: *invmatp = invmat = (int **)almat(col,col);
1.1 noro 3284: for ( i = 0; i < col; i++ )
3285: for ( j = 0, s = invmat[i], t = mat[i]; j < col; j++ )
3286: s[j] = t[col+index[j]];
3287: return 0;
3288: }
3289:
3290: void _addn(N,N,N);
3291: int _subn(N,N,N);
3292: void _muln(N,N,N);
3293:
1.24 noro 3294: void inner_product_int(Q *a,Q *b,int n,Q *r)
1.1 noro 3295: {
3296: int la,lb,i;
3297: int sgn,sgn1;
3298: N wm,wma,sum,t;
3299:
3300: for ( la = lb = 0, i = 0; i < n; i++ ) {
3301: if ( a[i] )
3302: if ( DN(a[i]) )
3303: error("inner_product_int : invalid argument");
3304: else
3305: la = MAX(PL(NM(a[i])),la);
3306: if ( b[i] )
3307: if ( DN(b[i]) )
3308: error("inner_product_int : invalid argument");
3309: else
3310: lb = MAX(PL(NM(b[i])),lb);
3311: }
3312: sgn = 0;
3313: sum= NALLOC(la+lb+2);
3314: bzero((char *)sum,(la+lb+3)*sizeof(unsigned int));
3315: wm = NALLOC(la+lb+2);
3316: wma = NALLOC(la+lb+2);
3317: for ( i = 0; i < n; i++ ) {
3318: if ( !a[i] || !b[i] )
3319: continue;
3320: _muln(NM(a[i]),NM(b[i]),wm);
3321: sgn1 = SGN(a[i])*SGN(b[i]);
3322: if ( !sgn ) {
3323: sgn = sgn1;
3324: t = wm; wm = sum; sum = t;
3325: } else if ( sgn == sgn1 ) {
3326: _addn(sum,wm,wma);
3327: if ( !PL(wma) )
3328: sgn = 0;
3329: t = wma; wma = sum; sum = t;
3330: } else {
3331: /* sgn*sum+sgn1*wm = sgn*(sum-wm) */
3332: sgn *= _subn(sum,wm,wma);
3333: t = wma; wma = sum; sum = t;
3334: }
3335: }
1.61 noro 3336: GCFREE(wm);
3337: GCFREE(wma);
1.1 noro 3338: if ( !sgn ) {
1.61 noro 3339: GCFREE(sum);
1.1 noro 3340: *r = 0;
3341: } else
3342: NTOQ(sum,sgn,*r);
3343: }
3344:
1.3 noro 3345: /* (k,l) element of a*b where a: .x n matrix, b: n x . integer matrix */
3346:
1.24 noro 3347: void inner_product_mat_int_mod(Q **a,int **b,int n,int k,int l,Q *r)
1.3 noro 3348: {
3349: int la,lb,i;
3350: int sgn,sgn1;
3351: N wm,wma,sum,t;
3352: Q aki;
3353: int bil,bilsgn;
3354: struct oN tn;
3355:
3356: for ( la = 0, i = 0; i < n; i++ ) {
3357: if ( aki = a[k][i] )
3358: if ( DN(aki) )
3359: error("inner_product_int : invalid argument");
3360: else
3361: la = MAX(PL(NM(aki)),la);
3362: }
3363: lb = 1;
3364: sgn = 0;
3365: sum= NALLOC(la+lb+2);
3366: bzero((char *)sum,(la+lb+3)*sizeof(unsigned int));
3367: wm = NALLOC(la+lb+2);
3368: wma = NALLOC(la+lb+2);
3369: for ( i = 0; i < n; i++ ) {
3370: if ( !(aki = a[k][i]) || !(bil = b[i][l]) )
3371: continue;
3372: tn.p = 1;
3373: if ( bil > 0 ) {
3374: tn.b[0] = bil; bilsgn = 1;
3375: } else {
3376: tn.b[0] = -bil; bilsgn = -1;
3377: }
3378: _muln(NM(aki),&tn,wm);
3379: sgn1 = SGN(aki)*bilsgn;
3380: if ( !sgn ) {
3381: sgn = sgn1;
3382: t = wm; wm = sum; sum = t;
3383: } else if ( sgn == sgn1 ) {
3384: _addn(sum,wm,wma);
3385: if ( !PL(wma) )
3386: sgn = 0;
3387: t = wma; wma = sum; sum = t;
3388: } else {
3389: /* sgn*sum+sgn1*wm = sgn*(sum-wm) */
3390: sgn *= _subn(sum,wm,wma);
3391: t = wma; wma = sum; sum = t;
3392: }
3393: }
1.61 noro 3394: GCFREE(wm);
3395: GCFREE(wma);
1.3 noro 3396: if ( !sgn ) {
1.61 noro 3397: GCFREE(sum);
1.3 noro 3398: *r = 0;
3399: } else
3400: NTOQ(sum,sgn,*r);
3401: }
3402:
1.24 noro 3403: void Pmul_mat_vect_int(NODE arg,VECT *rp)
1.1 noro 3404: {
3405: MAT mat;
3406: VECT vect,r;
3407: int row,col,i;
3408:
3409: mat = (MAT)ARG0(arg);
3410: vect = (VECT)ARG1(arg);
3411: row = mat->row;
3412: col = mat->col;
3413: MKVECT(r,row);
1.24 noro 3414: for ( i = 0; i < row; i++ ) {
3415: inner_product_int((Q *)mat->body[i],(Q *)vect->body,col,(Q *)&r->body[i]);
3416: }
1.1 noro 3417: *rp = r;
3418: }
3419:
1.24 noro 3420: void Pnbpoly_up2(NODE arg,GF2N *rp)
1.1 noro 3421: {
3422: int m,type,ret;
3423: UP2 r;
3424:
3425: m = QTOS((Q)ARG0(arg));
3426: type = QTOS((Q)ARG1(arg));
3427: ret = generate_ONB_polynomial(&r,m,type);
3428: if ( ret == 0 )
3429: MKGF2N(r,*rp);
3430: else
3431: *rp = 0;
3432: }
3433:
1.24 noro 3434: void Px962_irredpoly_up2(NODE arg,GF2N *rp)
1.1 noro 3435: {
1.24 noro 3436: int m,ret,w;
1.1 noro 3437: GF2N prev;
3438: UP2 r;
3439:
3440: m = QTOS((Q)ARG0(arg));
3441: prev = (GF2N)ARG1(arg);
3442: if ( !prev ) {
3443: w = (m>>5)+1; NEWUP2(r,w); r->w = 0;
3444: bzero((char *)r->b,w*sizeof(unsigned int));
3445: } else {
3446: r = prev->body;
3447: if ( degup2(r) != m ) {
3448: w = (m>>5)+1; NEWUP2(r,w); r->w = 0;
3449: bzero((char *)r->b,w*sizeof(unsigned int));
3450: }
3451: }
1.24 noro 3452: ret = _generate_irreducible_polynomial(r,m);
1.1 noro 3453: if ( ret == 0 )
3454: MKGF2N(r,*rp);
3455: else
3456: *rp = 0;
3457: }
3458:
1.24 noro 3459: void Pirredpoly_up2(NODE arg,GF2N *rp)
1.1 noro 3460: {
1.24 noro 3461: int m,ret,w;
1.1 noro 3462: GF2N prev;
3463: UP2 r;
3464:
3465: m = QTOS((Q)ARG0(arg));
3466: prev = (GF2N)ARG1(arg);
3467: if ( !prev ) {
3468: w = (m>>5)+1; NEWUP2(r,w); r->w = 0;
3469: bzero((char *)r->b,w*sizeof(unsigned int));
3470: } else {
3471: r = prev->body;
3472: if ( degup2(r) != m ) {
3473: w = (m>>5)+1; NEWUP2(r,w); r->w = 0;
3474: bzero((char *)r->b,w*sizeof(unsigned int));
3475: }
3476: }
1.24 noro 3477: ret = _generate_good_irreducible_polynomial(r,m);
1.1 noro 3478: if ( ret == 0 )
3479: MKGF2N(r,*rp);
3480: else
3481: *rp = 0;
3482: }
3483:
1.26 noro 3484: void Pmat_swap_row_destructive(NODE arg, MAT *m)
3485: {
3486: int i1,i2;
3487: pointer *t;
3488: MAT mat;
3489:
3490: asir_assert(ARG0(arg),O_MAT,"mat_swap_row_destructive");
3491: asir_assert(ARG1(arg),O_N,"mat_swap_row_destructive");
3492: asir_assert(ARG2(arg),O_N,"mat_swap_row_destructive");
3493: mat = (MAT)ARG0(arg);
3494: i1 = QTOS((Q)ARG1(arg));
3495: i2 = QTOS((Q)ARG2(arg));
3496: if ( i1 < 0 || i2 < 0 || i1 >= mat->row || i2 >= mat->row )
3497: error("mat_swap_row_destructive : Out of range");
3498: t = mat->body[i1];
3499: mat->body[i1] = mat->body[i2];
3500: mat->body[i2] = t;
3501: *m = mat;
3502: }
3503:
3504: void Pmat_swap_col_destructive(NODE arg, MAT *m)
3505: {
3506: int j1,j2,i,n;
3507: pointer *mi;
3508: pointer t;
3509: MAT mat;
3510:
3511: asir_assert(ARG0(arg),O_MAT,"mat_swap_col_destructive");
3512: asir_assert(ARG1(arg),O_N,"mat_swap_col_destructive");
3513: asir_assert(ARG2(arg),O_N,"mat_swap_col_destructive");
3514: mat = (MAT)ARG0(arg);
3515: j1 = QTOS((Q)ARG1(arg));
3516: j2 = QTOS((Q)ARG2(arg));
3517: if ( j1 < 0 || j2 < 0 || j1 >= mat->col || j2 >= mat->col )
3518: error("mat_swap_col_destructive : Out of range");
3519: n = mat->row;
3520: for ( i = 0; i < n; i++ ) {
3521: mi = mat->body[i];
3522: t = mi[j1]; mi[j1] = mi[j2]; mi[j2] = t;
3523: }
3524: *m = mat;
3525: }
1.1 noro 3526: /*
3527: * f = type 'type' normal polynomial of degree m if exists
3528: * IEEE P1363 A.7.2
3529: *
3530: * return value : 0 --- exists
3531: * 1 --- does not exist
3532: * -1 --- failure (memory allocation error)
3533: */
3534:
3535: int generate_ONB_polynomial(UP2 *rp,int m,int type)
3536: {
3537: int i,r;
3538: int w;
3539: UP2 f,f0,f1,f2,t;
3540:
3541: w = (m>>5)+1;
3542: switch ( type ) {
3543: case 1:
3544: if ( !TypeT_NB_check(m,1) ) return 1;
3545: NEWUP2(f,w); *rp = f; f->w = w;
3546: /* set all the bits */
3547: for ( i = 0; i < w; i++ )
3548: f->b[i] = 0xffffffff;
3549: /* mask the top word if necessary */
3550: if ( r = (m+1)&31 )
3551: f->b[w-1] &= (1<<r)-1;
3552: return 0;
3553: break;
3554: case 2:
3555: if ( !TypeT_NB_check(m,2) ) return 1;
3556: NEWUP2(f,w); *rp = f;
3557: W_NEWUP2(f0,w);
3558: W_NEWUP2(f1,w);
3559: W_NEWUP2(f2,w);
3560:
3561: /* recursion for genrating Type II normal polynomial */
3562:
3563: /* f0 = 1, f1 = t+1 */
3564: f0->w = 1; f0->b[0] = 1;
3565: f1->w = 1; f1->b[0] = 3;
3566: for ( i = 2; i <= m; i++ ) {
3567: /* f2 = t*f1+f0 */
3568: _bshiftup2(f1,-1,f2);
3569: _addup2_destructive(f2,f0);
3570: /* cyclic change of the variables */
3571: t = f0; f0 = f1; f1 = f2; f2 = t;
3572: }
3573: _copyup2(f1,f);
3574: return 0;
3575: break;
3576: default:
3577: return -1;
3578: break;
3579: }
3580: }
3581:
3582: /*
3583: * f = an irreducible trinomial or pentanomial of degree d 'after' f
3584: * return value : 0 --- exists
3585: * 1 --- does not exist (exhaustion)
3586: */
3587:
3588: int _generate_irreducible_polynomial(UP2 f,int d)
3589: {
3590: int ret,i,j,k,nz,i0,j0,k0;
3591: int w;
3592: unsigned int *fd;
3593:
3594: /*
3595: * if f = x^d+x^i+1 then i0 <- i, j0 <- 0, k0 <-0.
3596: * if f = x^d+x^k+x^j+x^i+1 (k>j>i) then i0 <- i, j0 <- j, k0 <-k.
3597: * otherwise i0,j0,k0 is set to 0.
3598: */
3599:
3600: fd = f->b;
3601: w = (d>>5)+1;
3602: if ( f->w && (d==degup2(f)) ) {
3603: for ( nz = 0, i = d; i >= 0; i-- )
3604: if ( fd[i>>5]&(1<<(i&31)) ) nz++;
3605: switch ( nz ) {
3606: case 3:
3607: for ( i0 = 1; !(fd[i0>>5]&(1<<(i0&31))) ; i0++ );
3608: /* reset i0-th bit */
3609: fd[i0>>5] &= ~(1<<(i0&31));
3610: j0 = k0 = 0;
3611: break;
3612: case 5:
3613: for ( i0 = 1; !(fd[i0>>5]&(1<<(i0&31))) ; i0++ );
3614: /* reset i0-th bit */
3615: fd[i0>>5] &= ~(1<<(i0&31));
3616: for ( j0 = i0+1; !(fd[j0>>5]&(1<<(j0&31))) ; j0++ );
3617: /* reset j0-th bit */
3618: fd[j0>>5] &= ~(1<<(j0&31));
3619: for ( k0 = j0+1; !(fd[k0>>5]&(1<<(k0&31))) ; k0++ );
3620: /* reset k0-th bit */
3621: fd[k0>>5] &= ~(1<<(k0&31));
3622: break;
3623: default:
3624: f->w = 0; break;
3625: }
3626: } else
3627: f->w = 0;
3628:
3629: if ( !f->w ) {
3630: fd = f->b;
3631: f->w = w; fd[0] |= 1; fd[d>>5] |= (1<<(d&31));
3632: i0 = j0 = k0 = 0;
3633: }
3634: /* if j0 > 0 then f is already a pentanomial */
3635: if ( j0 > 0 ) goto PENTA;
3636:
3637: /* searching for an irreducible trinomial */
3638:
3639: for ( i = 1; 2*i <= d; i++ ) {
3640: /* skip the polynomials 'before' f */
3641: if ( i < i0 ) continue;
3642: if ( i == i0 ) { i0 = 0; continue; }
3643: /* set i-th bit */
3644: fd[i>>5] |= (1<<(i&31));
3645: ret = irredcheck_dddup2(f);
3646: if ( ret == 1 ) return 0;
3647: /* reset i-th bit */
3648: fd[i>>5] &= ~(1<<(i&31));
3649: }
3650:
3651: /* searching for an irreducible pentanomial */
3652: PENTA:
3653: for ( i = 1; i < d; i++ ) {
3654: /* skip the polynomials 'before' f */
3655: if ( i < i0 ) continue;
3656: if ( i == i0 ) i0 = 0;
3657: /* set i-th bit */
3658: fd[i>>5] |= (1<<(i&31));
3659: for ( j = i+1; j < d; j++ ) {
3660: /* skip the polynomials 'before' f */
3661: if ( j < j0 ) continue;
3662: if ( j == j0 ) j0 = 0;
3663: /* set j-th bit */
3664: fd[j>>5] |= (1<<(j&31));
3665: for ( k = j+1; k < d; k++ ) {
3666: /* skip the polynomials 'before' f */
3667: if ( k < k0 ) continue;
3668: else if ( k == k0 ) { k0 = 0; continue; }
3669: /* set k-th bit */
3670: fd[k>>5] |= (1<<(k&31));
3671: ret = irredcheck_dddup2(f);
3672: if ( ret == 1 ) return 0;
3673: /* reset k-th bit */
3674: fd[k>>5] &= ~(1<<(k&31));
3675: }
3676: /* reset j-th bit */
3677: fd[j>>5] &= ~(1<<(j&31));
3678: }
3679: /* reset i-th bit */
3680: fd[i>>5] &= ~(1<<(i&31));
3681: }
3682: /* exhausted */
3683: return 1;
3684: }
3685:
3686: /*
3687: * f = an irreducible trinomial or pentanomial of degree d 'after' f
3688: *
3689: * searching strategy:
3690: * trinomial x^d+x^i+1:
3691: * i is as small as possible.
3692: * trinomial x^d+x^i+x^j+x^k+1:
3693: * i is as small as possible.
3694: * For such i, j is as small as possible.
3695: * For such i and j, 'k' is as small as possible.
3696: *
3697: * return value : 0 --- exists
3698: * 1 --- does not exist (exhaustion)
3699: */
3700:
3701: int _generate_good_irreducible_polynomial(UP2 f,int d)
3702: {
3703: int ret,i,j,k,nz,i0,j0,k0;
3704: int w;
3705: unsigned int *fd;
3706:
3707: /*
3708: * if f = x^d+x^i+1 then i0 <- i, j0 <- 0, k0 <-0.
3709: * if f = x^d+x^k+x^j+x^i+1 (k>j>i) then i0 <- i, j0 <- j, k0 <-k.
3710: * otherwise i0,j0,k0 is set to 0.
3711: */
3712:
3713: fd = f->b;
3714: w = (d>>5)+1;
3715: if ( f->w && (d==degup2(f)) ) {
3716: for ( nz = 0, i = d; i >= 0; i-- )
3717: if ( fd[i>>5]&(1<<(i&31)) ) nz++;
3718: switch ( nz ) {
3719: case 3:
3720: for ( i0 = 1; !(fd[i0>>5]&(1<<(i0&31))) ; i0++ );
3721: /* reset i0-th bit */
3722: fd[i0>>5] &= ~(1<<(i0&31));
3723: j0 = k0 = 0;
3724: break;
3725: case 5:
3726: for ( i0 = 1; !(fd[i0>>5]&(1<<(i0&31))) ; i0++ );
3727: /* reset i0-th bit */
3728: fd[i0>>5] &= ~(1<<(i0&31));
3729: for ( j0 = i0+1; !(fd[j0>>5]&(1<<(j0&31))) ; j0++ );
3730: /* reset j0-th bit */
3731: fd[j0>>5] &= ~(1<<(j0&31));
3732: for ( k0 = j0+1; !(fd[k0>>5]&(1<<(k0&31))) ; k0++ );
3733: /* reset k0-th bit */
3734: fd[k0>>5] &= ~(1<<(k0&31));
3735: break;
3736: default:
3737: f->w = 0; break;
3738: }
3739: } else
3740: f->w = 0;
3741:
3742: if ( !f->w ) {
3743: fd = f->b;
3744: f->w = w; fd[0] |= 1; fd[d>>5] |= (1<<(d&31));
3745: i0 = j0 = k0 = 0;
3746: }
3747: /* if j0 > 0 then f is already a pentanomial */
3748: if ( j0 > 0 ) goto PENTA;
3749:
3750: /* searching for an irreducible trinomial */
3751:
3752: for ( i = 1; 2*i <= d; i++ ) {
3753: /* skip the polynomials 'before' f */
3754: if ( i < i0 ) continue;
3755: if ( i == i0 ) { i0 = 0; continue; }
3756: /* set i-th bit */
3757: fd[i>>5] |= (1<<(i&31));
3758: ret = irredcheck_dddup2(f);
3759: if ( ret == 1 ) return 0;
3760: /* reset i-th bit */
3761: fd[i>>5] &= ~(1<<(i&31));
3762: }
3763:
3764: /* searching for an irreducible pentanomial */
3765: PENTA:
3766: for ( i = 3; i < d; i++ ) {
3767: /* skip the polynomials 'before' f */
3768: if ( i < i0 ) continue;
3769: if ( i == i0 ) i0 = 0;
3770: /* set i-th bit */
3771: fd[i>>5] |= (1<<(i&31));
3772: for ( j = 2; j < i; j++ ) {
3773: /* skip the polynomials 'before' f */
3774: if ( j < j0 ) continue;
3775: if ( j == j0 ) j0 = 0;
3776: /* set j-th bit */
3777: fd[j>>5] |= (1<<(j&31));
3778: for ( k = 1; k < j; k++ ) {
3779: /* skip the polynomials 'before' f */
3780: if ( k < k0 ) continue;
3781: else if ( k == k0 ) { k0 = 0; continue; }
3782: /* set k-th bit */
3783: fd[k>>5] |= (1<<(k&31));
3784: ret = irredcheck_dddup2(f);
3785: if ( ret == 1 ) return 0;
3786: /* reset k-th bit */
3787: fd[k>>5] &= ~(1<<(k&31));
3788: }
3789: /* reset j-th bit */
3790: fd[j>>5] &= ~(1<<(j&31));
3791: }
3792: /* reset i-th bit */
3793: fd[i>>5] &= ~(1<<(i&31));
3794: }
3795: /* exhausted */
3796: return 1;
1.3 noro 3797: }
3798:
1.24 noro 3799: void printqmat(Q **mat,int row,int col)
1.3 noro 3800: {
3801: int i,j;
3802:
3803: for ( i = 0; i < row; i++ ) {
3804: for ( j = 0; j < col; j++ ) {
1.8 noro 3805: printnum((Num)mat[i][j]); printf(" ");
1.3 noro 3806: }
3807: printf("\n");
3808: }
3809: }
3810:
1.24 noro 3811: void printimat(int **mat,int row,int col)
1.3 noro 3812: {
3813: int i,j;
3814:
3815: for ( i = 0; i < row; i++ ) {
3816: for ( j = 0; j < col; j++ ) {
3817: printf("%d ",mat[i][j]);
3818: }
3819: printf("\n");
3820: }
1.36 noro 3821: }
3822:
3823: void Pnd_det(NODE arg,P *rp)
3824: {
1.37 noro 3825: if ( argc(arg) == 1 )
3826: nd_det(0,ARG0(arg),rp);
3827: else
3828: nd_det(QTOS((Q)ARG1(arg)),ARG0(arg),rp);
1.1 noro 3829: }
1.59 ohara 3830:
1.62 ohara 3831: void Pmat_col(NODE arg,VECT *rp)
1.59 ohara 3832: {
3833: int i,j,n;
3834: MAT mat;
3835: VECT vect;
3836:
3837: asir_assert(ARG0(arg),O_MAT,"mat_col");
3838: asir_assert(ARG1(arg),O_N,"mat_col");
3839: mat = (MAT)ARG0(arg);
3840: j = QTOS((Q)ARG1(arg));
3841: if ( j < 0 || j >= mat->col) {
3842: error("mat_col : Out of range");
3843: }
3844: n = mat->row;
3845: MKVECT(vect,n);
3846: for(i=0; i<n; i++) {
3847: BDY(vect)[i] = BDY(mat)[i][j];
3848: }
3849: *rp = vect;
3850: }
1.71 ! noro 3851:
! 3852: NODE triangleq(NODE e)
! 3853: {
! 3854: int n,i,k;
! 3855: V v;
! 3856: VL vl;
! 3857: P *p;
! 3858: NODE r,r1;
! 3859:
! 3860: n = length(e);
! 3861: p = (P *)MALLOC(n*sizeof(P));
! 3862: for ( i = 0; i < n; i++, e = NEXT(e) ) p[i] = (P)BDY(e);
! 3863: i = 0;
! 3864: while ( 1 ) {
! 3865: for ( ; i < n && !p[i]; i++ );
! 3866: if ( i == n ) break;
! 3867: if ( OID(p[i]) == O_N ) return 0;
! 3868: v = p[i]->v;
! 3869: for ( k = i+1; k < n; k++ )
! 3870: if ( p[k] ) {
! 3871: if ( OID(p[k]) == O_N ) return 0;
! 3872: if ( p[k]->v == v ) p[k] = 0;
! 3873: }
! 3874: i++;
! 3875: }
! 3876: for ( r = 0, i = 0; i < n; i++ ) {
! 3877: if ( p[i] ) {
! 3878: MKNODE(r1,p[i],r); r = r1;
! 3879: }
! 3880: }
! 3881: return r;
! 3882: }
! 3883:
! 3884: void Ptriangleq(NODE arg,LIST *rp)
! 3885: {
! 3886: NODE ret;
! 3887:
! 3888: asir_assert(ARG0(arg),O_LIST,"sparseleq");
! 3889: ret = triangleq(BDY((LIST)ARG0(arg)));
! 3890: MKLIST(*rp,ret);
! 3891: }