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