![[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.65 ! noro 48: * $OpenXM: OpenXM_contrib2/asir2000/builtin/array.c,v 1.64 2013/11/05 02:55:02 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:
898: asir_assert(ARG0(arg),O_VECT,"vtol");
899: v = (VECT)ARG0(arg); len = v->len; a = BDY(v);
900: for ( i = len - 1, n = 0; i >= 0; i-- ) {
901: MKNODE(n1,a[i],n); n = n1;
902: }
903: MKLIST(*rp,n);
1.33 noro 904: }
905:
906: void Pltov(NODE arg,VECT *rp)
907: {
908: NODE n;
909: VECT v;
910: int len,i;
911:
912: asir_assert(ARG0(arg),O_LIST,"ltov");
913: n = (NODE)BDY((LIST)ARG0(arg));
914: len = length(n);
915: MKVECT(v,len);
916: for ( i = 0; i < len; i++, n = NEXT(n) )
917: BDY(v)[i] = BDY(n);
918: *rp = v;
1.1 noro 919: }
920:
1.24 noro 921: void Premainder(NODE arg,Obj *rp)
1.1 noro 922: {
923: Obj a;
924: VECT v,w;
925: MAT m,l;
926: pointer *vb,*wb;
927: pointer **mb,**lb;
928: int id,i,j,n,row,col,t,smd,sgn;
929: Q md,q;
930:
931: a = (Obj)ARG0(arg); md = (Q)ARG1(arg);
932: if ( !a )
933: *rp = 0;
934: else {
935: id = OID(a);
936: switch ( id ) {
937: case O_N:
938: case O_P:
939: cmp(md,(P)a,(P *)rp); break;
940: case O_VECT:
941: smd = QTOS(md);
942: v = (VECT)a; n = v->len; vb = v->body;
943: MKVECT(w,n); wb = w->body;
944: for ( i = 0; i < n; i++ ) {
945: if ( q = (Q)vb[i] ) {
946: sgn = SGN(q); t = rem(NM(q),smd);
947: STOQ(t,q);
948: if ( q )
949: SGN(q) = sgn;
950: }
951: wb[i] = (pointer)q;
952: }
953: *rp = (Obj)w;
954: break;
955: case O_MAT:
956: m = (MAT)a; row = m->row; col = m->col; mb = m->body;
957: MKMAT(l,row,col); lb = l->body;
958: for ( i = 0; i < row; i++ )
959: for ( j = 0, vb = mb[i], wb = lb[i]; j < col; j++ )
960: cmp(md,(P)vb[j],(P *)&wb[j]);
961: *rp = (Obj)l;
962: break;
963: default:
964: error("remainder : invalid argument");
965: }
966: }
967: }
968:
1.24 noro 969: void Psremainder(NODE arg,Obj *rp)
1.1 noro 970: {
971: Obj a;
972: VECT v,w;
973: MAT m,l;
974: pointer *vb,*wb;
975: pointer **mb,**lb;
976: unsigned int t,smd;
977: int id,i,j,n,row,col;
978: Q md,q;
979:
980: a = (Obj)ARG0(arg); md = (Q)ARG1(arg);
981: if ( !a )
982: *rp = 0;
983: else {
984: id = OID(a);
985: switch ( id ) {
986: case O_N:
987: case O_P:
988: cmp(md,(P)a,(P *)rp); break;
989: case O_VECT:
990: smd = QTOS(md);
991: v = (VECT)a; n = v->len; vb = v->body;
992: MKVECT(w,n); wb = w->body;
993: for ( i = 0; i < n; i++ ) {
994: if ( q = (Q)vb[i] ) {
995: t = (unsigned int)rem(NM(q),smd);
996: if ( SGN(q) < 0 )
997: t = (smd - t) % smd;
998: UTOQ(t,q);
999: }
1000: wb[i] = (pointer)q;
1001: }
1002: *rp = (Obj)w;
1003: break;
1004: case O_MAT:
1005: m = (MAT)a; row = m->row; col = m->col; mb = m->body;
1006: MKMAT(l,row,col); lb = l->body;
1007: for ( i = 0; i < row; i++ )
1008: for ( j = 0, vb = mb[i], wb = lb[i]; j < col; j++ )
1009: cmp(md,(P)vb[j],(P *)&wb[j]);
1010: *rp = (Obj)l;
1011: break;
1012: default:
1013: error("remainder : invalid argument");
1014: }
1015: }
1016: }
1017:
1.24 noro 1018: void Psize(NODE arg,LIST *rp)
1.1 noro 1019: {
1020:
1021: int n,m;
1022: Q q;
1023: NODE t,s;
1024:
1025: if ( !ARG0(arg) )
1026: t = 0;
1027: else {
1028: switch (OID(ARG0(arg))) {
1029: case O_VECT:
1030: n = ((VECT)ARG0(arg))->len;
1031: STOQ(n,q); MKNODE(t,q,0);
1032: break;
1033: case O_MAT:
1034: n = ((MAT)ARG0(arg))->row; m = ((MAT)ARG0(arg))->col;
1.43 saito 1035: STOQ(m,q); MKNODE(s,q,0); STOQ(n,q); MKNODE(t,q,s);
1036: break;
1037: case O_IMAT:
1038: n = ((IMAT)ARG0(arg))->row; m = ((IMAT)ARG0(arg))->col;
1.1 noro 1039: STOQ(m,q); MKNODE(s,q,0); STOQ(n,q); MKNODE(t,q,s);
1040: break;
1041: default:
1042: error("size : invalid argument"); break;
1043: }
1044: }
1045: MKLIST(*rp,t);
1046: }
1047:
1.24 noro 1048: void Pdet(NODE arg,P *rp)
1.1 noro 1049: {
1050: MAT m;
1051: int n,i,j,mod;
1052: P d;
1053: P **mat,**w;
1054:
1055: m = (MAT)ARG0(arg);
1056: asir_assert(m,O_MAT,"det");
1057: if ( m->row != m->col )
1058: error("det : non-square matrix");
1059: else if ( argc(arg) == 1 )
1060: detp(CO,(P **)BDY(m),m->row,rp);
1061: else {
1062: n = m->row; mod = QTOS((Q)ARG1(arg)); mat = (P **)BDY(m);
1063: w = (P **)almat_pointer(n,n);
1064: for ( i = 0; i < n; i++ )
1065: for ( j = 0; j < n; j++ )
1066: ptomp(mod,mat[i][j],&w[i][j]);
1067: detmp(CO,mod,w,n,&d);
1068: mptop(d,rp);
1.23 noro 1069: }
1070: }
1071:
1.24 noro 1072: void Pinvmat(NODE arg,LIST *rp)
1.23 noro 1073: {
1074: MAT m,r;
1075: int n,i,j,mod;
1076: P dn;
1077: P **mat,**imat,**w;
1078: NODE nd;
1079:
1080: m = (MAT)ARG0(arg);
1081: asir_assert(m,O_MAT,"invmat");
1082: if ( m->row != m->col )
1083: error("invmat : non-square matrix");
1084: else if ( argc(arg) == 1 ) {
1085: n = m->row;
1086: invmatp(CO,(P **)BDY(m),n,&imat,&dn);
1087: NEWMAT(r); r->row = n; r->col = n; r->body = (pointer **)imat;
1088: nd = mknode(2,r,dn);
1089: MKLIST(*rp,nd);
1090: } else {
1091: n = m->row; mod = QTOS((Q)ARG1(arg)); mat = (P **)BDY(m);
1092: w = (P **)almat_pointer(n,n);
1093: for ( i = 0; i < n; i++ )
1094: for ( j = 0; j < n; j++ )
1095: ptomp(mod,mat[i][j],&w[i][j]);
1096: #if 0
1097: detmp(CO,mod,w,n,&d);
1098: mptop(d,rp);
1099: #else
1100: error("not implemented yet");
1101: #endif
1.1 noro 1102: }
1.25 noro 1103: }
1104:
1105: /*
1106: input : a row x col matrix A
1107: A[I] <-> A[I][0]*x_0+A[I][1]*x_1+...
1108:
1.45 noro 1109: output : [B,D,R,C]
1.25 noro 1110: B : a rank(A) x col-rank(A) matrix
1.45 noro 1111: D : the denominator
1.25 noro 1112: R : a vector of length rank(A)
1113: C : a vector of length col-rank(A)
1.45 noro 1114: B[I] <-> D*x_{R[I]}+B[I][0]x_{C[0]}+B[I][1]x_{C[1]}+...
1.25 noro 1115: */
1116:
1117: void Pgeneric_gauss_elim(NODE arg,LIST *rp)
1118: {
1.48 noro 1119: NODE n0,opt,p;
1.25 noro 1120: MAT m,nm;
1121: int *ri,*ci;
1122: VECT rind,cind;
1123: Q dn,q;
1.62 ohara 1124: int i,row,col,t,rank;
1.48 noro 1125: int is_hensel = 0;
1126: char *key;
1127: Obj value;
1128:
1129: if ( current_option ) {
1130: for ( opt = current_option; opt; opt = NEXT(opt) ) {
1131: p = BDY((LIST)BDY(opt));
1132: key = BDY((STRING)BDY(p));
1133: value = (Obj)BDY(NEXT(p));
1134: if ( !strcmp(key,"hensel") && value ) {
1135: is_hensel = value ? 1 : 0;
1136: break;
1137: }
1138: }
1139: }
1.25 noro 1140: asir_assert(ARG0(arg),O_MAT,"generic_gauss_elim");
1141: m = (MAT)ARG0(arg);
1142: row = m->row; col = m->col;
1.48 noro 1143: if ( is_hensel )
1144: rank = generic_gauss_elim_hensel(m,&nm,&dn,&ri,&ci);
1145: else
1146: rank = generic_gauss_elim(m,&nm,&dn,&ri,&ci);
1.25 noro 1147: t = col-rank;
1148: MKVECT(rind,rank);
1149: MKVECT(cind,t);
1150: for ( i = 0; i < rank; i++ ) {
1151: STOQ(ri[i],q);
1152: BDY(rind)[i] = (pointer)q;
1153: }
1154: for ( i = 0; i < t; i++ ) {
1155: STOQ(ci[i],q);
1156: BDY(cind)[i] = (pointer)q;
1157: }
1158: n0 = mknode(4,nm,dn,rind,cind);
1159: MKLIST(*rp,n0);
1.1 noro 1160: }
1161:
1162: /*
1163: input : a row x col matrix A
1164: A[I] <-> A[I][0]*x_0+A[I][1]*x_1+...
1165:
1166: output : [B,R,C]
1167: B : a rank(A) x col-rank(A) matrix
1168: R : a vector of length rank(A)
1169: C : a vector of length col-rank(A)
1.47 noro 1170: RN : a vector of length rank(A) indicating useful rows
1171:
1.1 noro 1172: B[I] <-> x_{R[I]}+B[I][0]x_{C[0]}+B[I][1]x_{C[1]}+...
1173: */
1174:
1.24 noro 1175: void Pgeneric_gauss_elim_mod(NODE arg,LIST *rp)
1.1 noro 1176: {
1177: NODE n0;
1178: MAT m,mat;
1.47 noro 1179: VECT rind,cind,rnum;
1.1 noro 1180: Q **tmat;
1.47 noro 1181: int **wmat,**row0;
1182: Q *rib,*cib,*rnb;
1183: int *colstat,*p;
1.1 noro 1184: Q q;
1.24 noro 1185: int md,i,j,k,l,row,col,t,rank;
1.1 noro 1186:
1187: asir_assert(ARG0(arg),O_MAT,"generic_gauss_elim_mod");
1188: asir_assert(ARG1(arg),O_N,"generic_gauss_elim_mod");
1189: m = (MAT)ARG0(arg); md = QTOS((Q)ARG1(arg));
1190: row = m->row; col = m->col; tmat = (Q **)m->body;
1191: wmat = (int **)almat(row,col);
1.47 noro 1192:
1193: row0 = (int **)ALLOCA(row*sizeof(int *));
1194: for ( i = 0; i < row; i++ ) row0[i] = wmat[i];
1195:
1.1 noro 1196: colstat = (int *)MALLOC_ATOMIC(col*sizeof(int));
1197: for ( i = 0; i < row; i++ )
1198: for ( j = 0; j < col; j++ )
1199: if ( q = (Q)tmat[i][j] ) {
1200: t = rem(NM(q),md);
1201: if ( t && SGN(q) < 0 )
1202: t = (md - t) % md;
1203: wmat[i][j] = t;
1204: } else
1205: wmat[i][j] = 0;
1206: rank = generic_gauss_elim_mod(wmat,row,col,md,colstat);
1207:
1.47 noro 1208: MKVECT(rnum,rank);
1209: rnb = (Q *)rnum->body;
1210: for ( i = 0; i < rank; i++ )
1211: for ( j = 0, p = wmat[i]; j < row; j++ )
1212: if ( p == row0[j] )
1213: STOQ(j,rnb[i]);
1214:
1.1 noro 1215: MKMAT(mat,rank,col-rank);
1216: tmat = (Q **)mat->body;
1217: for ( i = 0; i < rank; i++ )
1218: for ( j = k = 0; j < col; j++ )
1219: if ( !colstat[j] ) {
1220: UTOQ(wmat[i][j],tmat[i][k]); k++;
1221: }
1222:
1223: MKVECT(rind,rank);
1224: MKVECT(cind,col-rank);
1225: rib = (Q *)rind->body; cib = (Q *)cind->body;
1226: for ( j = k = l = 0; j < col; j++ )
1227: if ( colstat[j] ) {
1228: STOQ(j,rib[k]); k++;
1229: } else {
1230: STOQ(j,cib[l]); l++;
1231: }
1.47 noro 1232: n0 = mknode(4,mat,rind,cind,rnum);
1.1 noro 1233: MKLIST(*rp,n0);
1234: }
1235:
1.24 noro 1236: void Pleqm(NODE arg,VECT *rp)
1.1 noro 1237: {
1238: MAT m;
1239: VECT vect;
1240: pointer **mat;
1241: Q *v;
1242: Q q;
1243: int **wmat;
1244: int md,i,j,row,col,t,n,status;
1245:
1246: asir_assert(ARG0(arg),O_MAT,"leqm");
1247: asir_assert(ARG1(arg),O_N,"leqm");
1248: m = (MAT)ARG0(arg); md = QTOS((Q)ARG1(arg));
1249: row = m->row; col = m->col; mat = m->body;
1250: wmat = (int **)almat(row,col);
1251: for ( i = 0; i < row; i++ )
1252: for ( j = 0; j < col; j++ )
1253: if ( q = (Q)mat[i][j] ) {
1254: t = rem(NM(q),md);
1255: if ( SGN(q) < 0 )
1256: t = (md - t) % md;
1257: wmat[i][j] = t;
1258: } else
1259: wmat[i][j] = 0;
1260: status = gauss_elim_mod(wmat,row,col,md);
1261: if ( status < 0 )
1262: *rp = 0;
1263: else if ( status > 0 )
1264: *rp = (VECT)ONE;
1265: else {
1266: n = col - 1;
1267: MKVECT(vect,n);
1268: for ( i = 0, v = (Q *)vect->body; i < n; i++ ) {
1269: t = (md-wmat[i][n])%md; STOQ(t,v[i]);
1270: }
1271: *rp = vect;
1272: }
1273: }
1274:
1.24 noro 1275: int gauss_elim_mod(int **mat,int row,int col,int md)
1.1 noro 1276: {
1277: int i,j,k,inv,a,n;
1278: int *t,*pivot;
1279:
1280: n = col - 1;
1281: for ( j = 0; j < n; j++ ) {
1282: for ( i = j; i < row && !mat[i][j]; i++ );
1283: if ( i == row )
1284: return 1;
1285: if ( i != j ) {
1286: t = mat[i]; mat[i] = mat[j]; mat[j] = t;
1287: }
1288: pivot = mat[j];
1289: inv = invm(pivot[j],md);
1290: for ( k = j; k <= n; k++ ) {
1291: /* pivot[k] = dmar(pivot[k],inv,0,md); */
1292: DMAR(pivot[k],inv,0,md,pivot[k])
1293: }
1294: for ( i = 0; i < row; i++ ) {
1295: t = mat[i];
1296: if ( i != j && (a = t[j]) )
1297: for ( k = j, a = md - a; k <= n; k++ ) {
1.8 noro 1298: unsigned int tk;
1.1 noro 1299: /* t[k] = dmar(pivot[k],a,t[k],md); */
1.8 noro 1300: DMAR(pivot[k],a,t[k],md,tk)
1301: t[k] = tk;
1.1 noro 1302: }
1303: }
1304: }
1305: for ( i = n; i < row && !mat[i][n]; i++ );
1306: if ( i == row )
1307: return 0;
1308: else
1309: return -1;
1310: }
1311:
1.4 noro 1312: struct oEGT eg_mod,eg_elim,eg_elim1,eg_elim2,eg_chrem,eg_gschk,eg_intrat,eg_symb;
1.31 noro 1313: struct oEGT eg_conv;
1.1 noro 1314:
1.24 noro 1315: int generic_gauss_elim(MAT mat,MAT *nm,Q *dn,int **rindp,int **cindp)
1.1 noro 1316: {
1317: int **wmat;
1318: Q **bmat;
1319: N **tmat;
1320: Q *bmi;
1321: N *tmi;
1322: Q q;
1323: int *wmi;
1324: int *colstat,*wcolstat,*rind,*cind;
1325: int row,col,ind,md,i,j,k,l,t,t1,rank,rank0,inv;
1326: N m1,m2,m3,s,u;
1327: MAT r,crmat;
1328: struct oEGT tmp0,tmp1;
1329: struct oEGT eg_mod_split,eg_elim_split,eg_chrem_split;
1330: struct oEGT eg_intrat_split,eg_gschk_split;
1331: int ret;
1332:
1333: init_eg(&eg_mod_split); init_eg(&eg_chrem_split);
1334: init_eg(&eg_elim_split); init_eg(&eg_intrat_split);
1335: init_eg(&eg_gschk_split);
1336: bmat = (Q **)mat->body;
1337: row = mat->row; col = mat->col;
1338: wmat = (int **)almat(row,col);
1339: colstat = (int *)MALLOC_ATOMIC(col*sizeof(int));
1340: wcolstat = (int *)MALLOC_ATOMIC(col*sizeof(int));
1341: for ( ind = 0; ; ind++ ) {
1.11 noro 1342: if ( DP_Print ) {
1.2 noro 1343: fprintf(asir_out,"."); fflush(asir_out);
1344: }
1.12 noro 1345: md = get_lprime(ind);
1.1 noro 1346: get_eg(&tmp0);
1347: for ( i = 0; i < row; i++ )
1348: for ( j = 0, bmi = bmat[i], wmi = wmat[i]; j < col; j++ )
1349: if ( q = (Q)bmi[j] ) {
1350: t = rem(NM(q),md);
1351: if ( t && SGN(q) < 0 )
1352: t = (md - t) % md;
1353: wmi[j] = t;
1354: } else
1355: wmi[j] = 0;
1356: get_eg(&tmp1);
1357: add_eg(&eg_mod,&tmp0,&tmp1);
1358: add_eg(&eg_mod_split,&tmp0,&tmp1);
1359: get_eg(&tmp0);
1360: rank = generic_gauss_elim_mod(wmat,row,col,md,wcolstat);
1361: get_eg(&tmp1);
1362: add_eg(&eg_elim,&tmp0,&tmp1);
1363: add_eg(&eg_elim_split,&tmp0,&tmp1);
1364: if ( !ind ) {
1365: RESET:
1366: UTON(md,m1);
1367: rank0 = rank;
1368: bcopy(wcolstat,colstat,col*sizeof(int));
1369: MKMAT(crmat,rank,col-rank);
1370: MKMAT(r,rank,col-rank); *nm = r;
1371: tmat = (N **)crmat->body;
1372: for ( i = 0; i < rank; i++ )
1373: for ( j = k = 0, tmi = tmat[i], wmi = wmat[i]; j < col; j++ )
1374: if ( !colstat[j] ) {
1375: UTON(wmi[j],tmi[k]); k++;
1376: }
1377: } else {
1378: if ( rank < rank0 ) {
1.11 noro 1379: if ( DP_Print ) {
1.1 noro 1380: fprintf(asir_out,"lower rank matrix; continuing...\n");
1.2 noro 1381: fflush(asir_out);
1382: }
1.1 noro 1383: continue;
1384: } else if ( rank > rank0 ) {
1.11 noro 1385: if ( DP_Print ) {
1.1 noro 1386: fprintf(asir_out,"higher rank matrix; resetting...\n");
1.2 noro 1387: fflush(asir_out);
1388: }
1.1 noro 1389: goto RESET;
1390: } else {
1391: for ( j = 0; (j<col) && (colstat[j]==wcolstat[j]); j++ );
1392: if ( j < col ) {
1.11 noro 1393: if ( DP_Print ) {
1.1 noro 1394: fprintf(asir_out,"inconsitent colstat; resetting...\n");
1.2 noro 1395: fflush(asir_out);
1396: }
1.1 noro 1397: goto RESET;
1398: }
1399: }
1400:
1401: get_eg(&tmp0);
1402: inv = invm(rem(m1,md),md);
1403: UTON(md,m2); muln(m1,m2,&m3);
1404: for ( i = 0; i < rank; i++ )
1405: for ( j = k = 0, tmi = tmat[i], wmi = wmat[i]; j < col; j++ )
1406: if ( !colstat[j] ) {
1407: if ( tmi[k] ) {
1408: /* f3 = f1+m1*(m1 mod m2)^(-1)*(f2 - f1 mod m2) */
1409: t = rem(tmi[k],md);
1410: if ( wmi[j] >= t )
1411: t = wmi[j]-t;
1412: else
1413: t = md-(t-wmi[j]);
1414: DMAR(t,inv,0,md,t1)
1415: UTON(t1,u);
1416: muln(m1,u,&s);
1417: addn(tmi[k],s,&u); tmi[k] = u;
1418: } else if ( wmi[j] ) {
1419: /* f3 = m1*(m1 mod m2)^(-1)*f2 */
1420: DMAR(wmi[j],inv,0,md,t)
1421: UTON(t,u);
1422: muln(m1,u,&s); tmi[k] = s;
1423: }
1424: k++;
1425: }
1426: m1 = m3;
1427: get_eg(&tmp1);
1428: add_eg(&eg_chrem,&tmp0,&tmp1);
1429: add_eg(&eg_chrem_split,&tmp0,&tmp1);
1430:
1431: get_eg(&tmp0);
1.38 noro 1432: if ( ind % F4_INTRAT_PERIOD )
1.13 noro 1433: ret = 0;
1434: else
1435: ret = intmtoratm(crmat,m1,*nm,dn);
1.1 noro 1436: get_eg(&tmp1);
1437: add_eg(&eg_intrat,&tmp0,&tmp1);
1438: add_eg(&eg_intrat_split,&tmp0,&tmp1);
1439: if ( ret ) {
1440: *rindp = rind = (int *)MALLOC_ATOMIC(rank*sizeof(int));
1441: *cindp = cind = (int *)MALLOC_ATOMIC((col-rank)*sizeof(int));
1442: for ( j = k = l = 0; j < col; j++ )
1443: if ( colstat[j] )
1444: rind[k++] = j;
1445: else
1446: cind[l++] = j;
1447: get_eg(&tmp0);
1.3 noro 1448: if ( gensolve_check(mat,*nm,*dn,rind,cind) ) {
1449: get_eg(&tmp1);
1450: add_eg(&eg_gschk,&tmp0,&tmp1);
1451: add_eg(&eg_gschk_split,&tmp0,&tmp1);
1.11 noro 1452: if ( DP_Print ) {
1.3 noro 1453: print_eg("Mod",&eg_mod_split);
1454: print_eg("Elim",&eg_elim_split);
1455: print_eg("ChRem",&eg_chrem_split);
1456: print_eg("IntRat",&eg_intrat_split);
1457: print_eg("Check",&eg_gschk_split);
1458: fflush(asir_out);
1459: }
1460: return rank;
1461: }
1462: }
1463: }
1464: }
1465: }
1466:
1.64 noro 1467: void lu_dec_cr(MAT mat,MAT lu,Q *dn,int **perm);
1468:
1.53 noro 1469: /* XXX broken */
1.64 noro 1470: void lu_dec_cr(MAT mat,MAT lu,Q *dn,int **perm)
1.53 noro 1471: {
1472: Q **a0,**b;
1473: Q *aiq;
1474: N **a;
1475: N *ai;
1476: Q q,q1,dn2,a1,q0,bik;
1477: MAT m;
1478: unsigned int md;
1479: int n,ind,i,j,rank,t,inv,t1,ret,min,k;
1480: int **w;
1481: int *wi,*rinfo0,*rinfo;
1482: N m1,m2,m3,u,s;
1483:
1484: a0 = (Q **)mat->body;
1485: n = mat->row;
1486: if ( n != mat->col )
1487: error("lu_dec_cr : non-square matrix");
1488: w = (int **)almat(n,n);
1489: MKMAT(m,n,n);
1490: a = (N **)m->body;
1491: UTON(1,m1);
1492: rinfo0 = 0;
1493: ind = 0;
1494: while ( 1 ) {
1495: md = get_lprime(ind);
1496: /* mat mod md */
1497: for ( i = 0; i < n; i++ )
1498: for ( j = 0, aiq = a0[i], wi = w[i]; j < n; j++ )
1499: if ( q = aiq[j] ) {
1500: t = rem(NM(q),md);
1501: if ( t && SGN(q) < 0 )
1502: t = (md - t) % md;
1503: wi[j] = t;
1504: } else
1505: wi[j] = 0;
1506:
1507: if ( !lu_mod((unsigned int **)w,n,md,&rinfo) ) continue;
1508: printf("."); fflush(stdout);
1509: if ( !rinfo0 )
1510: *perm = rinfo0 = rinfo;
1511: else {
1512: for ( i = 0; i < n; i++ )
1513: if ( rinfo[i] != rinfo0[i] ) break;
1514: if ( i < n ) continue;
1515: }
1516: if ( UNIN(m1) ) {
1517: for ( i = 0; i < n; i++ )
1518: for ( j = 0, ai = a[i], wi = w[i]; j < n; j++ ) {
1519: UTON(wi[j],u); ai[j] = u;
1520: }
1521: UTON(md,m1);
1522: } else {
1523: inv = invm(rem(m1,md),md);
1524: UTON(md,m2); muln(m1,m2,&m3);
1525: for ( i = 0; i < n; i++ )
1526: for ( j = 0, ai = a[i], wi = w[i]; j < n; j++ )
1527: if ( ai[i] ) {
1528: /* f3 = f1+m1*(m1 mod m2)^(-1)*(f2 - f1 mod m2) */
1529: t = rem(ai[j],md);
1530: if ( wi[j] >= t )
1531: t = wi[j]-t;
1532: else
1533: t = md-(t-wi[j]);
1534: DMAR(t,inv,0,md,t1)
1535: UTON(t1,u);
1536: muln(m1,u,&s);
1537: addn(ai[j],s,&u); ai[j] = u;
1538: } else if ( wi[j] ) {
1539: /* f3 = m1*(m1 mod m2)^(-1)*f2 */
1540: DMAR(wi[j],inv,0,md,t)
1541: UTON(t,u);
1542: muln(m1,u,&s); ai[j] = s;
1543: }
1544: m1 = m3;
1545: }
1546: if ( (++ind%8) == 0 ) {
1547: ret = intmtoratm(m,m1,lu,dn);
1548: if ( ret ) {
1549: b = (Q **)lu->body;
1550: mulq(*dn,*dn,&dn2);
1551: for ( i = 0; i < n; i++ ) {
1552: for ( j = 0; j < n; j++ ) {
1553: q = 0;
1554: min = MIN(i,j);
1555: for ( k = 0; k <= min; k++ ) {
1556: bik = k==i ? *dn : b[i][k];
1557: mulq(bik,b[k][j],&q0);
1558: addq(q,q0,&q1); q = q1;
1559: }
1560: mulq(a0[rinfo0[i]][j],dn2,&q1);
1561: if ( cmpq(q,q1) ) break;
1562: }
1563: if ( j < n ) break;
1564: }
1565: if ( i == n )
1566: return;
1567: }
1568: }
1569: }
1570: }
1571:
1.64 noro 1572: void nmat(N **m,int n)
1.53 noro 1573: {
1574: int i,j;
1575:
1576: for ( i = 0; i < n; i++ ) {
1577: for ( j = 0; j < n; j++ ) {
1578: printn(m[i][j]); printf(" ");
1579: }
1580: printf("\n");
1581: }
1582: }
1583:
1.24 noro 1584: int generic_gauss_elim_hensel(MAT mat,MAT *nmmat,Q *dn,int **rindp,int **cindp)
1.3 noro 1585: {
1586: MAT bmat,xmat;
1587: Q **a0,**a,**b,**x,**nm;
1588: Q *ai,*bi,*xi;
1589: int row,col;
1590: int **w;
1591: int *wi;
1592: int **wc;
1593: Q mdq,q,s,u;
1594: N tn;
1595: int ind,md,i,j,k,l,li,ri,rank;
1596: unsigned int t;
1597: int *cinfo,*rinfo;
1598: int *rind,*cind;
1599: int count;
1.41 noro 1600: int ret;
1601: struct oEGT eg_mul,eg_inv,eg_intrat,eg_check,tmp0,tmp1;
1.39 noro 1602: int period;
1.44 noro 1603: int *wx,*ptr;
1604: int wxsize,nsize;
1605: N wn;
1606: Q wq;
1.3 noro 1607:
1608: a0 = (Q **)mat->body;
1609: row = mat->row; col = mat->col;
1610: w = (int **)almat(row,col);
1611: for ( ind = 0; ; ind++ ) {
1.12 noro 1612: md = get_lprime(ind);
1.3 noro 1613: STOQ(md,mdq);
1614: for ( i = 0; i < row; i++ )
1615: for ( j = 0, ai = a0[i], wi = w[i]; j < col; j++ )
1616: if ( q = (Q)ai[j] ) {
1617: t = rem(NM(q),md);
1618: if ( t && SGN(q) < 0 )
1619: t = (md - t) % md;
1620: wi[j] = t;
1621: } else
1622: wi[j] = 0;
1623:
1.52 noro 1624: if ( DP_Print > 3 ) {
1.48 noro 1625: fprintf(asir_out,"LU decomposition.."); fflush(asir_out);
1626: }
1.27 noro 1627: rank = find_lhs_and_lu_mod((unsigned int **)w,row,col,md,&rinfo,&cinfo);
1.52 noro 1628: if ( DP_Print > 3 ) {
1.48 noro 1629: fprintf(asir_out,"done.\n"); fflush(asir_out);
1630: }
1.3 noro 1631: a = (Q **)almat_pointer(rank,rank); /* lhs mat */
1632: MKMAT(bmat,rank,col-rank); b = (Q **)bmat->body; /* lhs mat */
1633: for ( j = li = ri = 0; j < col; j++ )
1634: if ( cinfo[j] ) {
1635: /* the column is in lhs */
1636: for ( i = 0; i < rank; i++ ) {
1637: w[i][li] = w[i][j];
1638: a[i][li] = a0[rinfo[i]][j];
1639: }
1640: li++;
1641: } else {
1642: /* the column is in rhs */
1643: for ( i = 0; i < rank; i++ )
1644: b[i][ri] = a0[rinfo[i]][j];
1645: ri++;
1646: }
1647:
1648: /* solve Ax+B=0; A: rank x rank, B: rank x ri */
1649: MKMAT(xmat,rank,ri); x = (Q **)(xmat)->body;
1650: MKMAT(*nmmat,rank,ri); nm = (Q **)(*nmmat)->body;
1651: /* use the right part of w as work area */
1652: /* ri = col - rank */
1653: wc = (int **)almat(rank,ri);
1654: for ( i = 0; i < rank; i++ )
1655: wc[i] = w[i]+rank;
1656: *rindp = rind = (int *)MALLOC_ATOMIC(rank*sizeof(int));
1657: *cindp = cind = (int *)MALLOC_ATOMIC((ri)*sizeof(int));
1658:
1659: init_eg(&eg_mul); init_eg(&eg_inv);
1.41 noro 1660: init_eg(&eg_check); init_eg(&eg_intrat);
1.39 noro 1661: period = F4_INTRAT_PERIOD;
1.44 noro 1662: nsize = period;
1663: wxsize = rank*ri*nsize;
1664: wx = (int *)MALLOC_ATOMIC(wxsize*sizeof(int));
1665: for ( i = 0; i < wxsize; i++ ) wx[i] = 0;
1666: for ( q = ONE, count = 0; ; ) {
1.52 noro 1667: if ( DP_Print > 3 )
1.41 noro 1668: fprintf(stderr,"o");
1.3 noro 1669: /* wc = -b mod md */
1.44 noro 1670: get_eg(&tmp0);
1.3 noro 1671: for ( i = 0; i < rank; i++ )
1672: for ( j = 0, bi = b[i], wi = wc[i]; j < ri; j++ )
1673: if ( u = (Q)bi[j] ) {
1674: t = rem(NM(u),md);
1675: if ( t && SGN(u) > 0 )
1676: t = (md - t) % md;
1677: wi[j] = t;
1678: } else
1679: wi[j] = 0;
1.44 noro 1680: /* wc = A^(-1)wc; wc is not normalized */
1681: solve_by_lu_mod(w,rank,md,wc,ri,0);
1682: /* wx += q*wc */
1683: ptr = wx;
1684: for ( i = 0; i < rank; i++ )
1685: for ( j = 0, wi = wc[i]; j < ri; j++ ) {
1686: if ( wi[j] )
1687: muln_1(BD(NM(q)),PL(NM(q)),wi[j],ptr);
1688: ptr += nsize;
1689: }
1690: count++;
1.1 noro 1691: get_eg(&tmp1);
1.3 noro 1692: add_eg(&eg_inv,&tmp0,&tmp1);
1693: get_eg(&tmp0);
1694: for ( i = 0; i < rank; i++ )
1695: for ( j = 0; j < ri; j++ ) {
1696: inner_product_mat_int_mod(a,wc,rank,i,j,&u);
1697: addq(b[i][j],u,&s);
1698: if ( s ) {
1699: t = divin(NM(s),md,&tn);
1700: if ( t )
1701: error("generic_gauss_elim_hensel:incosistent");
1702: NTOQ(tn,SGN(s),b[i][j]);
1703: } else
1704: b[i][j] = 0;
1705: }
1706: get_eg(&tmp1);
1707: add_eg(&eg_mul,&tmp0,&tmp1);
1708: /* q = q*md */
1709: mulq(q,mdq,&u); q = u;
1.44 noro 1710: if ( count == period ) {
1.41 noro 1711: get_eg(&tmp0);
1.44 noro 1712: ptr = wx;
1713: for ( i = 0; i < rank; i++ )
1714: for ( j = 0, xi = x[i]; j < ri;
1715: j++, ptr += nsize ) {
1716: for ( k = nsize-1; k >= 0 && !ptr[k]; k-- );
1717: if ( k >= 0 ) {
1718: wn = NALLOC(k+1);
1719: PL(wn) = k+1;
1720: for ( l = 0; l <= k; l++ ) BD(wn)[l] = (unsigned int)ptr[l];
1721: NTOQ(wn,1,wq);
1722: subq(xi[j],wq,&u); xi[j] = u;
1723: }
1724: }
1.41 noro 1725: ret = intmtoratm_q(xmat,NM(q),*nmmat,dn);
1726: get_eg(&tmp1); add_eg(&eg_intrat,&tmp0,&tmp1);
1727: if ( ret ) {
1.50 noro 1728: rind = (int *)MALLOC_ATOMIC(rank*sizeof(int));
1729: cind = (int *)MALLOC_ATOMIC((col-rank)*sizeof(int));
1.39 noro 1730: for ( j = k = l = 0; j < col; j++ )
1731: if ( cinfo[j] )
1732: rind[k++] = j;
1733: else
1.50 noro 1734: cind[l++] = j;
1735: get_eg(&tmp0);
1736: ret = gensolve_check(mat,*nmmat,*dn,rind,cind);
1737: get_eg(&tmp1); add_eg(&eg_check,&tmp0,&tmp1);
1738: if ( ret ) {
1739: if ( DP_Print > 3 ) {
1740: fprintf(stderr,"\n");
1741: print_eg("INV",&eg_inv);
1742: print_eg("MUL",&eg_mul);
1743: print_eg("INTRAT",&eg_intrat);
1744: print_eg("CHECK",&eg_check);
1745: fflush(asir_out);
1746: }
1747: *rindp = rind;
1748: *cindp = cind;
1749: for ( j = k = 0; j < col; j++ )
1750: if ( !cinfo[j] )
1751: cind[k++] = j;
1752: return rank;
1753: }
1754: } else {
1755: period = period*3/2;
1756: count = 0;
1757: nsize += period;
1758: wxsize += rank*ri*nsize;
1759: wx = (int *)REALLOC(wx,wxsize*sizeof(int));
1760: for ( i = 0; i < wxsize; i++ ) wx[i] = 0;
1761: }
1762: }
1763: }
1764: }
1765: }
1766:
1.55 noro 1767: int generic_gauss_elim_hensel_dalg(MAT mat,DP *mb,MAT *nmmat,Q *dn,int **rindp,int **cindp)
1.50 noro 1768: {
1769: MAT bmat,xmat;
1770: Q **a0,**a,**b,**x,**nm;
1771: Q *ai,*bi,*xi;
1772: int row,col;
1773: int **w;
1774: int *wi;
1775: int **wc;
1776: Q mdq,q,s,u;
1777: N tn;
1778: int ind,md,i,j,k,l,li,ri,rank;
1779: unsigned int t;
1780: int *cinfo,*rinfo;
1781: int *rind,*cind;
1782: int count;
1783: int ret;
1784: struct oEGT eg_mul,eg_inv,eg_intrat,eg_check,tmp0,tmp1;
1785: int period;
1786: int *wx,*ptr;
1787: int wxsize,nsize;
1788: N wn;
1789: Q wq;
1790: NumberField nf;
1791: DP m;
1792: int col1;
1793:
1794: a0 = (Q **)mat->body;
1795: row = mat->row; col = mat->col;
1796: w = (int **)almat(row,col);
1797: for ( ind = 0; ; ind++ ) {
1798: md = get_lprime(ind);
1799: STOQ(md,mdq);
1800: for ( i = 0; i < row; i++ )
1801: for ( j = 0, ai = a0[i], wi = w[i]; j < col; j++ )
1802: if ( q = (Q)ai[j] ) {
1803: t = rem(NM(q),md);
1804: if ( t && SGN(q) < 0 )
1805: t = (md - t) % md;
1806: wi[j] = t;
1807: } else
1808: wi[j] = 0;
1809:
1810: if ( DP_Print ) {
1811: fprintf(asir_out,"LU decomposition.."); fflush(asir_out);
1812: }
1813: rank = find_lhs_and_lu_mod((unsigned int **)w,row,col,md,&rinfo,&cinfo);
1814: if ( DP_Print ) {
1815: fprintf(asir_out,"done.\n"); fflush(asir_out);
1816: }
1817: for ( i = 0; i < col-1; i++ ) {
1818: if ( !cinfo[i] ) {
1819: m = mb[i];
1820: for ( j = i+1; j < col-1; j++ )
1821: if ( dp_redble(mb[j],m) )
1822: cinfo[j] = -1;
1823: }
1824: }
1825: a = (Q **)almat_pointer(rank,rank); /* lhs mat */
1826: MKMAT(bmat,rank,col-rank); b = (Q **)bmat->body; /* lhs mat */
1827: for ( j = li = ri = 0; j < col; j++ )
1828: if ( cinfo[j] > 0 ) {
1829: /* the column is in lhs */
1830: for ( i = 0; i < rank; i++ ) {
1831: w[i][li] = w[i][j];
1832: a[i][li] = a0[rinfo[i]][j];
1833: }
1834: li++;
1835: } else if ( !cinfo[j] ) {
1836: /* the column is in rhs */
1837: for ( i = 0; i < rank; i++ )
1838: b[i][ri] = a0[rinfo[i]][j];
1839: ri++;
1840: }
1841:
1842: /* solve Ax+B=0; A: rank x rank, B: rank x ri */
1843: MKMAT(xmat,rank,ri); x = (Q **)(xmat)->body;
1844: MKMAT(*nmmat,rank,ri); nm = (Q **)(*nmmat)->body;
1845: /* use the right part of w as work area */
1846: wc = (int **)almat(rank,ri);
1847: for ( i = 0; i < rank; i++ )
1848: wc[i] = w[i]+rank;
1849: *rindp = rind = (int *)MALLOC_ATOMIC(rank*sizeof(int));
1850: *cindp = cind = (int *)MALLOC_ATOMIC((col-rank)*sizeof(int));
1851: init_eg(&eg_mul); init_eg(&eg_inv);
1852: init_eg(&eg_check); init_eg(&eg_intrat);
1853: period = F4_INTRAT_PERIOD;
1854: nsize = period;
1855: wxsize = rank*ri*nsize;
1856: wx = (int *)MALLOC_ATOMIC(wxsize*sizeof(int));
1857: for ( i = 0; i < wxsize; i++ ) wx[i] = 0;
1858: for ( q = ONE, count = 0; ; ) {
1859: if ( DP_Print )
1860: fprintf(stderr,"o");
1861: /* wc = -b mod md */
1862: get_eg(&tmp0);
1863: for ( i = 0; i < rank; i++ )
1864: for ( j = 0, bi = b[i], wi = wc[i]; j < ri; j++ )
1865: if ( u = (Q)bi[j] ) {
1866: t = rem(NM(u),md);
1867: if ( t && SGN(u) > 0 )
1868: t = (md - t) % md;
1869: wi[j] = t;
1870: } else
1871: wi[j] = 0;
1872: /* wc = A^(-1)wc; wc is not normalized */
1873: solve_by_lu_mod(w,rank,md,wc,ri,0);
1874: /* wx += q*wc */
1875: ptr = wx;
1876: for ( i = 0; i < rank; i++ )
1877: for ( j = 0, wi = wc[i]; j < ri; j++ ) {
1878: if ( wi[j] )
1879: muln_1(BD(NM(q)),PL(NM(q)),wi[j],ptr);
1880: ptr += nsize;
1881: }
1882: count++;
1883: get_eg(&tmp1);
1884: add_eg(&eg_inv,&tmp0,&tmp1);
1885: get_eg(&tmp0);
1886: for ( i = 0; i < rank; i++ )
1887: for ( j = 0; j < ri; j++ ) {
1888: inner_product_mat_int_mod(a,wc,rank,i,j,&u);
1889: addq(b[i][j],u,&s);
1890: if ( s ) {
1891: t = divin(NM(s),md,&tn);
1892: if ( t )
1893: error("generic_gauss_elim_hensel:incosistent");
1894: NTOQ(tn,SGN(s),b[i][j]);
1895: } else
1896: b[i][j] = 0;
1897: }
1898: get_eg(&tmp1);
1899: add_eg(&eg_mul,&tmp0,&tmp1);
1900: /* q = q*md */
1901: mulq(q,mdq,&u); q = u;
1902: if ( count == period ) {
1903: get_eg(&tmp0);
1904: ptr = wx;
1905: for ( i = 0; i < rank; i++ )
1906: for ( j = 0, xi = x[i]; j < ri;
1907: j++, ptr += nsize ) {
1908: for ( k = nsize-1; k >= 0 && !ptr[k]; k-- );
1909: if ( k >= 0 ) {
1910: wn = NALLOC(k+1);
1911: PL(wn) = k+1;
1912: for ( l = 0; l <= k; l++ ) BD(wn)[l] = (unsigned int)ptr[l];
1913: NTOQ(wn,1,wq);
1914: subq(xi[j],wq,&u); xi[j] = u;
1915: }
1916: }
1917: ret = intmtoratm_q(xmat,NM(q),*nmmat,dn);
1918: get_eg(&tmp1); add_eg(&eg_intrat,&tmp0,&tmp1);
1919: if ( ret ) {
1920: for ( j = k = l = 0; j < col; j++ )
1921: if ( cinfo[j] > 0 )
1922: rind[k++] = j;
1923: else if ( !cinfo[j] )
1.39 noro 1924: cind[l++] = j;
1.41 noro 1925: get_eg(&tmp0);
1926: ret = gensolve_check(mat,*nmmat,*dn,rind,cind);
1927: get_eg(&tmp1); add_eg(&eg_check,&tmp0,&tmp1);
1928: if ( ret ) {
1.42 noro 1929: if ( DP_Print > 3 ) {
1.40 noro 1930: fprintf(stderr,"\n");
1931: print_eg("INV",&eg_inv);
1932: print_eg("MUL",&eg_mul);
1.41 noro 1933: print_eg("INTRAT",&eg_intrat);
1934: print_eg("CHECK",&eg_check);
1.40 noro 1935: fflush(asir_out);
1936: }
1.39 noro 1937: return rank;
1938: }
1.44 noro 1939: } else {
1940: period = period*3/2;
1941: count = 0;
1942: nsize += period;
1943: wxsize += rank*ri*nsize;
1944: wx = (int *)REALLOC(wx,wxsize*sizeof(int));
1945: for ( i = 0; i < wxsize; i++ ) wx[i] = 0;
1946: }
1.41 noro 1947: }
1.1 noro 1948: }
1949: }
1950: }
1951:
1952: int f4_nocheck;
1953:
1.24 noro 1954: int gensolve_check(MAT mat,MAT nm,Q dn,int *rind,int *cind)
1.1 noro 1955: {
1956: int row,col,rank,clen,i,j,k,l;
1.24 noro 1957: Q s,t;
1.1 noro 1958: Q *w;
1959: Q *mati,*nmk;
1960:
1961: if ( f4_nocheck )
1962: return 1;
1963: row = mat->row; col = mat->col;
1964: rank = nm->row; clen = nm->col;
1965: w = (Q *)MALLOC(clen*sizeof(Q));
1966: for ( i = 0; i < row; i++ ) {
1967: mati = (Q *)mat->body[i];
1968: #if 1
1969: bzero(w,clen*sizeof(Q));
1970: for ( k = 0; k < rank; k++ )
1971: for ( l = 0, nmk = (Q *)nm->body[k]; l < clen; l++ ) {
1972: mulq(mati[rind[k]],nmk[l],&t);
1973: addq(w[l],t,&s); w[l] = s;
1974: }
1975: for ( j = 0; j < clen; j++ ) {
1976: mulq(dn,mati[cind[j]],&t);
1977: if ( cmpq(w[j],t) )
1978: break;
1979: }
1980: #else
1981: for ( j = 0; j < clen; j++ ) {
1982: for ( k = 0, s = 0; k < rank; k++ ) {
1983: mulq(mati[rind[k]],nm->body[k][j],&t);
1984: addq(s,t,&u); s = u;
1985: }
1986: mulq(dn,mati[cind[j]],&t);
1987: if ( cmpq(s,t) )
1988: break;
1989: }
1990: #endif
1991: if ( j != clen )
1992: break;
1993: }
1994: if ( i != row )
1995: return 0;
1996: else
1997: return 1;
1998: }
1999:
2000: /* assuming 0 < c < m */
2001:
1.24 noro 2002: int inttorat(N c,N m,N b,int *sgnp,N *nmp,N *dnp)
1.1 noro 2003: {
1.24 noro 2004: Q qq,t,u1,v1,r1;
2005: N q,u2,v2,r2;
1.1 noro 2006:
2007: u1 = 0; v1 = ONE; u2 = m; v2 = c;
2008: while ( cmpn(v2,b) >= 0 ) {
2009: divn(u2,v2,&q,&r2); u2 = v2; v2 = r2;
2010: NTOQ(q,1,qq); mulq(qq,v1,&t); subq(u1,t,&r1); u1 = v1; v1 = r1;
2011: }
2012: if ( cmpn(NM(v1),b) >= 0 )
2013: return 0;
2014: else {
2015: *nmp = v2;
2016: *dnp = NM(v1);
2017: *sgnp = SGN(v1);
2018: return 1;
2019: }
2020: }
2021:
2022: /* mat->body = N ** */
2023:
1.24 noro 2024: int intmtoratm(MAT mat,N md,MAT nm,Q *dn)
1.1 noro 2025: {
2026: N t,s,b;
1.24 noro 2027: Q dn0,dn1,nm1,q;
1.1 noro 2028: int i,j,k,l,row,col;
2029: Q **rmat;
2030: N **tmat;
2031: N *tmi;
2032: Q *nmk;
2033: N u,unm,udn;
2034: int sgn,ret;
2035:
1.3 noro 2036: if ( UNIN(md) )
2037: return 0;
1.1 noro 2038: row = mat->row; col = mat->col;
2039: bshiftn(md,1,&t);
2040: isqrt(t,&s);
2041: bshiftn(s,64,&b);
2042: if ( !b )
2043: b = ONEN;
2044: dn0 = ONE;
2045: tmat = (N **)mat->body;
2046: rmat = (Q **)nm->body;
2047: for ( i = 0; i < row; i++ )
2048: for ( j = 0, tmi = tmat[i]; j < col; j++ )
2049: if ( tmi[j] ) {
2050: muln(tmi[j],NM(dn0),&s);
2051: remn(s,md,&u);
2052: ret = inttorat(u,md,b,&sgn,&unm,&udn);
2053: if ( !ret )
2054: return 0;
2055: else {
2056: NTOQ(unm,sgn,nm1);
2057: NTOQ(udn,1,dn1);
2058: if ( !UNIQ(dn1) ) {
2059: for ( k = 0; k < i; k++ )
2060: for ( l = 0, nmk = rmat[k]; l < col; l++ ) {
2061: mulq(nmk[l],dn1,&q); nmk[l] = q;
2062: }
2063: for ( l = 0, nmk = rmat[i]; l < j; l++ ) {
2064: mulq(nmk[l],dn1,&q); nmk[l] = q;
2065: }
2066: }
2067: rmat[i][j] = nm1;
2068: mulq(dn0,dn1,&q); dn0 = q;
2069: }
2070: }
2071: *dn = dn0;
2072: return 1;
2073: }
2074:
1.3 noro 2075: /* mat->body = Q ** */
2076:
1.24 noro 2077: int intmtoratm_q(MAT mat,N md,MAT nm,Q *dn)
1.3 noro 2078: {
2079: N t,s,b;
1.24 noro 2080: Q dn0,dn1,nm1,q;
1.3 noro 2081: int i,j,k,l,row,col;
2082: Q **rmat;
2083: Q **tmat;
2084: Q *tmi;
2085: Q *nmk;
2086: N u,unm,udn;
2087: int sgn,ret;
2088:
2089: if ( UNIN(md) )
2090: return 0;
2091: row = mat->row; col = mat->col;
2092: bshiftn(md,1,&t);
2093: isqrt(t,&s);
2094: bshiftn(s,64,&b);
2095: if ( !b )
2096: b = ONEN;
2097: dn0 = ONE;
2098: tmat = (Q **)mat->body;
2099: rmat = (Q **)nm->body;
2100: for ( i = 0; i < row; i++ )
2101: for ( j = 0, tmi = tmat[i]; j < col; j++ )
2102: if ( tmi[j] ) {
2103: muln(NM(tmi[j]),NM(dn0),&s);
2104: remn(s,md,&u);
2105: ret = inttorat(u,md,b,&sgn,&unm,&udn);
2106: if ( !ret )
2107: return 0;
2108: else {
2109: if ( SGN(tmi[j])<0 )
2110: sgn = -sgn;
2111: NTOQ(unm,sgn,nm1);
2112: NTOQ(udn,1,dn1);
2113: if ( !UNIQ(dn1) ) {
2114: for ( k = 0; k < i; k++ )
2115: for ( l = 0, nmk = rmat[k]; l < col; l++ ) {
2116: mulq(nmk[l],dn1,&q); nmk[l] = q;
2117: }
2118: for ( l = 0, nmk = rmat[i]; l < j; l++ ) {
2119: mulq(nmk[l],dn1,&q); nmk[l] = q;
2120: }
2121: }
2122: rmat[i][j] = nm1;
2123: mulq(dn0,dn1,&q); dn0 = q;
2124: }
2125: }
2126: *dn = dn0;
2127: return 1;
2128: }
2129:
1.4 noro 2130: #define ONE_STEP1 if ( zzz = *s ) { DMAR(zzz,hc,*tj,md,*tj) } tj++; s++;
2131:
1.24 noro 2132: void reduce_reducers_mod(int **mat,int row,int col,int md)
1.4 noro 2133: {
2134: int i,j,k,l,hc,zzz;
2135: int *t,*s,*tj,*ind;
2136:
2137: /* reduce the reducers */
2138: ind = (int *)ALLOCA(row*sizeof(int));
2139: for ( i = 0; i < row; i++ ) {
2140: t = mat[i];
2141: for ( j = 0; j < col && !t[j]; j++ );
2142: /* register the position of the head term */
2143: ind[i] = j;
2144: for ( l = i-1; l >= 0; l-- ) {
2145: /* reduce mat[i] by mat[l] */
2146: if ( hc = t[ind[l]] ) {
2147: /* mat[i] = mat[i]-hc*mat[l] */
2148: j = ind[l];
2149: s = mat[l]+j;
2150: tj = t+j;
2151: hc = md-hc;
2152: k = col-j;
2153: for ( ; k >= 64; k -= 64 ) {
2154: ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
2155: ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
2156: ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
2157: ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
2158: ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
2159: ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
2160: ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
2161: ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
2162: ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
2163: ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
2164: ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
2165: ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
2166: ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
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: }
1.16 noro 2171: for ( ; k > 0; k-- ) {
1.4 noro 2172: if ( zzz = *s ) { DMAR(zzz,hc,*tj,md,*tj) } tj++; s++;
2173: }
2174: }
2175: }
2176: }
2177: }
2178:
2179: /*
2180: mat[i] : reducers (i=0,...,nred-1)
2181: spolys (i=nred,...,row-1)
2182: mat[0] < mat[1] < ... < mat[nred-1] w.r.t the term order
2183: 1. reduce the reducers
2184: 2. reduce spolys by the reduced reducers
2185: */
2186:
1.24 noro 2187: void pre_reduce_mod(int **mat,int row,int col,int nred,int md)
1.4 noro 2188: {
2189: int i,j,k,l,hc,inv;
2190: int *t,*s,*tk,*ind;
2191:
2192: #if 1
2193: /* reduce the reducers */
2194: ind = (int *)ALLOCA(row*sizeof(int));
2195: for ( i = 0; i < nred; i++ ) {
2196: /* make mat[i] monic and mat[i] by mat[0],...,mat[i-1] */
2197: t = mat[i];
2198: for ( j = 0; j < col && !t[j]; j++ );
2199: /* register the position of the head term */
2200: ind[i] = j;
2201: inv = invm(t[j],md);
2202: for ( k = j; k < col; k++ )
2203: if ( t[k] )
2204: DMAR(t[k],inv,0,md,t[k])
2205: for ( l = i-1; l >= 0; l-- ) {
2206: /* reduce mat[i] by mat[l] */
2207: if ( hc = t[ind[l]] ) {
2208: /* mat[i] = mat[i]-hc*mat[l] */
2209: for ( k = ind[l], hc = md-hc, s = mat[l]+k, tk = t+k;
2210: k < col; k++, tk++, s++ )
2211: if ( *s )
2212: DMAR(*s,hc,*tk,md,*tk)
2213: }
2214: }
2215: }
2216: /* reduce the spolys */
2217: for ( i = nred; i < row; i++ ) {
2218: t = mat[i];
2219: for ( l = nred-1; l >= 0; l-- ) {
2220: /* reduce mat[i] by mat[l] */
2221: if ( hc = t[ind[l]] ) {
2222: /* mat[i] = mat[i]-hc*mat[l] */
2223: for ( k = ind[l], hc = md-hc, s = mat[l]+k, tk = t+k;
2224: k < col; k++, tk++, s++ )
2225: if ( *s )
2226: DMAR(*s,hc,*tk,md,*tk)
2227: }
2228: }
2229: }
2230: #endif
2231: }
2232: /*
2233: mat[i] : reducers (i=0,...,nred-1)
2234: mat[0] < mat[1] < ... < mat[nred-1] w.r.t the term order
2235: */
2236:
1.24 noro 2237: void reduce_sp_by_red_mod(int *sp,int **redmat,int *ind,int nred,int col,int md)
1.4 noro 2238: {
2239: int i,j,k,hc,zzz;
1.24 noro 2240: int *s,*tj;
1.4 noro 2241:
2242: /* reduce the spolys by redmat */
2243: for ( i = nred-1; i >= 0; i-- ) {
2244: /* reduce sp by redmat[i] */
2245: if ( hc = sp[ind[i]] ) {
2246: /* sp = sp-hc*redmat[i] */
2247: j = ind[i];
2248: hc = md-hc;
2249: s = redmat[i]+j;
2250: tj = sp+j;
1.16 noro 2251: for ( k = col-j; k > 0; k-- ) {
1.4 noro 2252: if ( zzz = *s ) { DMAR(zzz,hc,*tj,md,*tj) } tj++; s++;
1.15 noro 2253: }
2254: }
1.17 noro 2255: }
2256: }
2257:
2258: /*
1.15 noro 2259: mat[i] : compressed reducers (i=0,...,nred-1)
2260: mat[0] < mat[1] < ... < mat[nred-1] w.r.t the term order
2261: */
2262:
1.24 noro 2263: void red_by_compress(int m,unsigned int *p,unsigned int *r,
2264: unsigned int *ri,unsigned int hc,int len)
1.18 noro 2265: {
1.19 noro 2266: unsigned int up,lo;
1.18 noro 2267: unsigned int dmy;
2268: unsigned int *pj;
2269:
1.21 noro 2270: p[*ri] = 0; r++; ri++;
2271: for ( len--; len; len--, r++, ri++ ) {
2272: pj = p+ *ri;
2273: DMA(*r,hc,*pj,up,lo);
1.18 noro 2274: if ( up ) {
2275: DSAB(m,up,lo,dmy,*pj);
2276: } else
2277: *pj = lo;
2278: }
2279: }
2280:
2281: /* p -= hc*r */
2282:
1.24 noro 2283: void red_by_vect(int m,unsigned int *p,unsigned int *r,unsigned int hc,int len)
1.18 noro 2284: {
1.58 ohara 2285: unsigned int up,lo,dmy;
1.18 noro 2286:
2287: *p++ = 0; r++; len--;
2288: for ( ; len; len--, r++, p++ )
2289: if ( *r ) {
1.20 noro 2290: DMA(*r,hc,*p,up,lo);
1.18 noro 2291: if ( up ) {
2292: DSAB(m,up,lo,dmy,*p);
2293: } else
2294: *p = lo;
2295: }
2296: }
2297:
1.32 noro 2298: void red_by_vect_sf(int m,unsigned int *p,unsigned int *r,unsigned int hc,int len)
2299: {
2300: *p++ = 0; r++; len--;
2301: for ( ; len; len--, r++, p++ )
2302: if ( *r )
2303: *p = _addsf(_mulsf(*r,hc),*p);
2304: }
2305:
1.21 noro 2306: extern unsigned int **psca;
2307:
1.24 noro 2308: void reduce_sp_by_red_mod_compress (int *sp,CDP *redmat,int *ind,
2309: int nred,int col,int md)
1.15 noro 2310: {
1.24 noro 2311: int i,len;
1.15 noro 2312: CDP ri;
1.24 noro 2313: unsigned int hc;
1.18 noro 2314: unsigned int *usp;
1.15 noro 2315:
1.18 noro 2316: usp = (unsigned int *)sp;
1.15 noro 2317: /* reduce the spolys by redmat */
2318: for ( i = nred-1; i >= 0; i-- ) {
2319: /* reduce sp by redmat[i] */
1.18 noro 2320: usp[ind[i]] %= md;
2321: if ( hc = usp[ind[i]] ) {
1.15 noro 2322: /* sp = sp-hc*redmat[i] */
2323: hc = md-hc;
2324: ri = redmat[i];
2325: len = ri->len;
1.21 noro 2326: red_by_compress(md,usp,psca[ri->psindex],ri->body,hc,len);
1.4 noro 2327: }
2328: }
1.18 noro 2329: for ( i = 0; i < col; i++ )
1.24 noro 2330: if ( usp[i] >= (unsigned int)md )
1.18 noro 2331: usp[i] %= md;
1.4 noro 2332: }
2333:
2334: #define ONE_STEP2 if ( zzz = *pk ) { DMAR(zzz,a,*tk,md,*tk) } pk++; tk++;
2335:
1.24 noro 2336: int generic_gauss_elim_mod(int **mat0,int row,int col,int md,int *colstat)
1.1 noro 2337: {
1.24 noro 2338: int i,j,k,l,inv,a,rank;
2339: unsigned int *t,*pivot,*pk;
1.18 noro 2340: unsigned int **mat;
1.1 noro 2341:
1.18 noro 2342: mat = (unsigned int **)mat0;
1.1 noro 2343: for ( rank = 0, j = 0; j < col; j++ ) {
1.18 noro 2344: for ( i = rank; i < row; i++ )
2345: mat[i][j] %= md;
2346: for ( i = rank; i < row; i++ )
2347: if ( mat[i][j] )
2348: break;
1.1 noro 2349: if ( i == row ) {
2350: colstat[j] = 0;
2351: continue;
2352: } else
2353: colstat[j] = 1;
2354: if ( i != rank ) {
2355: t = mat[i]; mat[i] = mat[rank]; mat[rank] = t;
2356: }
2357: pivot = mat[rank];
2358: inv = invm(pivot[j],md);
1.4 noro 2359: for ( k = j, pk = pivot+k; k < col; k++, pk++ )
2360: if ( *pk ) {
1.24 noro 2361: if ( *pk >= (unsigned int)md )
1.18 noro 2362: *pk %= md;
1.4 noro 2363: DMAR(*pk,inv,0,md,*pk)
1.1 noro 2364: }
2365: for ( i = rank+1; i < row; i++ ) {
2366: t = mat[i];
1.18 noro 2367: if ( a = t[j] )
2368: red_by_vect(md,t+j,pivot+j,md-a,col-j);
1.1 noro 2369: }
2370: rank++;
2371: }
2372: for ( j = col-1, l = rank-1; j >= 0; j-- )
2373: if ( colstat[j] ) {
2374: pivot = mat[l];
2375: for ( i = 0; i < l; i++ ) {
2376: t = mat[i];
1.18 noro 2377: t[j] %= md;
2378: if ( a = t[j] )
2379: red_by_vect(md,t+j,pivot+j,md-a,col-j);
1.1 noro 2380: }
2381: l--;
1.18 noro 2382: }
2383: for ( j = 0, l = 0; l < rank; j++ )
2384: if ( colstat[j] ) {
2385: t = mat[l];
2386: for ( k = j; k < col; k++ )
1.24 noro 2387: if ( t[k] >= (unsigned int)md )
1.18 noro 2388: t[k] %= md;
2389: l++;
1.32 noro 2390: }
2391: return rank;
2392: }
2393:
1.65 ! noro 2394: int generic_gauss_elim_mod2(int **mat0,int row,int col,int md,int *colstat,int *rowstat)
! 2395: {
! 2396: int i,j,k,l,inv,a,rank;
! 2397: unsigned int *t,*pivot,*pk;
! 2398: unsigned int **mat;
! 2399:
! 2400: for ( i = 0; i < row; i++ ) rowstat[i] = i;
! 2401: mat = (unsigned int **)mat0;
! 2402: for ( rank = 0, j = 0; j < col; j++ ) {
! 2403: for ( i = rank; i < row; i++ )
! 2404: mat[i][j] %= md;
! 2405: for ( i = rank; i < row; i++ )
! 2406: if ( mat[i][j] )
! 2407: break;
! 2408: if ( i == row ) {
! 2409: colstat[j] = 0;
! 2410: continue;
! 2411: } else
! 2412: colstat[j] = 1;
! 2413: if ( i != rank ) {
! 2414: t = mat[i]; mat[i] = mat[rank]; mat[rank] = t;
! 2415: k = rowstat[i]; rowstat[i] = rowstat[rank]; rowstat[rank] = k;
! 2416: }
! 2417: pivot = mat[rank];
! 2418: inv = invm(pivot[j],md);
! 2419: for ( k = j, pk = pivot+k; k < col; k++, pk++ )
! 2420: if ( *pk ) {
! 2421: if ( *pk >= (unsigned int)md )
! 2422: *pk %= md;
! 2423: DMAR(*pk,inv,0,md,*pk)
! 2424: }
! 2425: for ( i = rank+1; i < row; i++ ) {
! 2426: t = mat[i];
! 2427: if ( a = t[j] )
! 2428: red_by_vect(md,t+j,pivot+j,md-a,col-j);
! 2429: }
! 2430: rank++;
! 2431: }
! 2432: for ( j = col-1, l = rank-1; j >= 0; j-- )
! 2433: if ( colstat[j] ) {
! 2434: pivot = mat[l];
! 2435: for ( i = 0; i < l; i++ ) {
! 2436: t = mat[i];
! 2437: t[j] %= md;
! 2438: if ( a = t[j] )
! 2439: red_by_vect(md,t+j,pivot+j,md-a,col-j);
! 2440: }
! 2441: l--;
! 2442: }
! 2443: for ( j = 0, l = 0; l < rank; j++ )
! 2444: if ( colstat[j] ) {
! 2445: t = mat[l];
! 2446: for ( k = j; k < col; k++ )
! 2447: if ( t[k] >= (unsigned int)md )
! 2448: t[k] %= md;
! 2449: l++;
! 2450: }
! 2451: return rank;
! 2452: }
! 2453:
1.32 noro 2454: int generic_gauss_elim_sf(int **mat0,int row,int col,int md,int *colstat)
2455: {
2456: int i,j,k,l,inv,a,rank;
2457: unsigned int *t,*pivot,*pk;
2458: unsigned int **mat;
2459:
2460: mat = (unsigned int **)mat0;
2461: for ( rank = 0, j = 0; j < col; j++ ) {
2462: for ( i = rank; i < row; i++ )
2463: if ( mat[i][j] )
2464: break;
2465: if ( i == row ) {
2466: colstat[j] = 0;
2467: continue;
2468: } else
2469: colstat[j] = 1;
2470: if ( i != rank ) {
2471: t = mat[i]; mat[i] = mat[rank]; mat[rank] = t;
2472: }
2473: pivot = mat[rank];
2474: inv = _invsf(pivot[j]);
2475: for ( k = j, pk = pivot+k; k < col; k++, pk++ )
2476: if ( *pk )
2477: *pk = _mulsf(*pk,inv);
2478: for ( i = rank+1; i < row; i++ ) {
2479: t = mat[i];
2480: if ( a = t[j] )
2481: red_by_vect_sf(md,t+j,pivot+j,_chsgnsf(a),col-j);
2482: }
2483: rank++;
2484: }
2485: for ( j = col-1, l = rank-1; j >= 0; j-- )
2486: if ( colstat[j] ) {
2487: pivot = mat[l];
2488: for ( i = 0; i < l; i++ ) {
2489: t = mat[i];
2490: if ( a = t[j] )
2491: red_by_vect_sf(md,t+j,pivot+j,_chsgnsf(a),col-j);
2492: }
2493: l--;
1.1 noro 2494: }
2495: return rank;
2496: }
2497:
2498: /* LU decomposition; a[i][i] = 1/U[i][i] */
2499:
1.24 noro 2500: int lu_gfmmat(GFMMAT mat,unsigned int md,int *perm)
1.1 noro 2501: {
2502: int row,col;
1.24 noro 2503: int i,j,k;
1.1 noro 2504: unsigned int *t,*pivot;
2505: unsigned int **a;
2506: unsigned int inv,m;
2507:
2508: row = mat->row; col = mat->col;
2509: a = mat->body;
2510: bzero(perm,row*sizeof(int));
2511:
2512: for ( i = 0; i < row; i++ )
2513: perm[i] = i;
2514: for ( k = 0; k < col; k++ ) {
2515: for ( i = k; i < row && !a[i][k]; i++ );
2516: if ( i == row )
2517: return 0;
2518: if ( i != k ) {
2519: j = perm[i]; perm[i] = perm[k]; perm[k] = j;
2520: t = a[i]; a[i] = a[k]; a[k] = t;
2521: }
2522: pivot = a[k];
2523: pivot[k] = inv = invm(pivot[k],md);
2524: for ( i = k+1; i < row; i++ ) {
2525: t = a[i];
2526: if ( m = t[k] ) {
2527: DMAR(inv,m,0,md,t[k])
2528: for ( j = k+1, m = md - t[k]; j < col; j++ )
2529: if ( pivot[j] ) {
1.8 noro 2530: unsigned int tj;
2531:
2532: DMAR(m,pivot[j],t[j],md,tj)
2533: t[j] = tj;
1.1 noro 2534: }
2535: }
2536: }
2537: }
2538: return 1;
2539: }
2540:
1.3 noro 2541: /*
2542: Input
2543: a: a row x col matrix
2544: md : a modulus
2545:
2546: Output:
2547: return : d = the rank of mat
2548: a[0..(d-1)][0..(d-1)] : LU decomposition (a[i][i] = 1/U[i][i])
2549: rinfo: array of length row
2550: cinfo: array of length col
2551: i-th row in new a <-> rinfo[i]-th row in old a
2552: cinfo[j]=1 <=> j-th column is contained in the LU decomp.
2553: */
2554:
1.24 noro 2555: int find_lhs_and_lu_mod(unsigned int **a,int row,int col,
2556: unsigned int md,int **rinfo,int **cinfo)
1.3 noro 2557: {
1.24 noro 2558: int i,j,k,d;
1.3 noro 2559: int *rp,*cp;
2560: unsigned int *t,*pivot;
2561: unsigned int inv,m;
2562:
2563: *rinfo = rp = (int *)MALLOC_ATOMIC(row*sizeof(int));
2564: *cinfo = cp = (int *)MALLOC_ATOMIC(col*sizeof(int));
2565: for ( i = 0; i < row; i++ )
2566: rp[i] = i;
2567: for ( k = 0, d = 0; k < col; k++ ) {
2568: for ( i = d; i < row && !a[i][k]; i++ );
2569: if ( i == row ) {
2570: cp[k] = 0;
2571: continue;
2572: } else
2573: cp[k] = 1;
2574: if ( i != d ) {
2575: j = rp[i]; rp[i] = rp[d]; rp[d] = j;
2576: t = a[i]; a[i] = a[d]; a[d] = t;
2577: }
2578: pivot = a[d];
2579: pivot[k] = inv = invm(pivot[k],md);
2580: for ( i = d+1; i < row; i++ ) {
2581: t = a[i];
2582: if ( m = t[k] ) {
2583: DMAR(inv,m,0,md,t[k])
2584: for ( j = k+1, m = md - t[k]; j < col; j++ )
2585: if ( pivot[j] ) {
1.8 noro 2586: unsigned int tj;
2587: DMAR(m,pivot[j],t[j],md,tj)
2588: t[j] = tj;
1.3 noro 2589: }
2590: }
2591: }
2592: d++;
2593: }
2594: return d;
2595: }
2596:
1.53 noro 2597: int lu_mod(unsigned int **a,int n,unsigned int md,int **rinfo)
2598: {
2599: int i,j,k;
2600: int *rp;
2601: unsigned int *t,*pivot;
2602: unsigned int inv,m;
2603:
2604: *rinfo = rp = (int *)MALLOC_ATOMIC(n*sizeof(int));
2605: for ( i = 0; i < n; i++ ) rp[i] = i;
2606: for ( k = 0; k < n; k++ ) {
2607: for ( i = k; i < n && !a[i][k]; i++ );
2608: if ( i == n ) return 0;
2609: if ( i != k ) {
2610: j = rp[i]; rp[i] = rp[k]; rp[k] = j;
2611: t = a[i]; a[i] = a[k]; a[k] = t;
2612: }
2613: pivot = a[k];
2614: inv = invm(pivot[k],md);
2615: for ( i = k+1; i < n; i++ ) {
2616: t = a[i];
2617: if ( m = t[k] ) {
2618: DMAR(inv,m,0,md,t[k])
2619: for ( j = k+1, m = md - t[k]; j < n; j++ )
2620: if ( pivot[j] ) {
2621: unsigned int tj;
2622: DMAR(m,pivot[j],t[j],md,tj)
2623: t[j] = tj;
2624: }
2625: }
2626: }
2627: }
2628: return 1;
2629: }
2630:
1.3 noro 2631: /*
2632: Input
2633: a : n x n matrix; a result of LU-decomposition
2634: md : modulus
2635: b : n x l matrix
2636: Output
2637: b = a^(-1)b
2638: */
2639:
1.44 noro 2640: void solve_by_lu_mod(int **a,int n,int md,int **b,int l,int normalize)
1.3 noro 2641: {
2642: unsigned int *y,*c;
2643: int i,j,k;
2644: unsigned int t,m,m2;
2645:
2646: y = (int *)MALLOC_ATOMIC(n*sizeof(int));
2647: c = (int *)MALLOC_ATOMIC(n*sizeof(int));
2648: m2 = md>>1;
2649: for ( k = 0; k < l; k++ ) {
2650: /* copy b[.][k] to c */
2651: for ( i = 0; i < n; i++ )
2652: c[i] = (unsigned int)b[i][k];
2653: /* solve Ly=c */
2654: for ( i = 0; i < n; i++ ) {
2655: for ( t = c[i], j = 0; j < i; j++ )
2656: if ( a[i][j] ) {
2657: m = md - a[i][j];
2658: DMAR(m,y[j],t,md,t)
2659: }
2660: y[i] = t;
2661: }
2662: /* solve Uc=y */
2663: for ( i = n-1; i >= 0; i-- ) {
2664: for ( t = y[i], j =i+1; j < n; j++ )
2665: if ( a[i][j] ) {
2666: m = md - a[i][j];
2667: DMAR(m,c[j],t,md,t)
2668: }
2669: /* a[i][i] = 1/U[i][i] */
2670: DMAR(t,a[i][i],0,md,c[i])
2671: }
2672: /* copy c to b[.][k] with normalization */
1.44 noro 2673: if ( normalize )
2674: for ( i = 0; i < n; i++ )
2675: b[i][k] = (int)(c[i]>m2 ? c[i]-md : c[i]);
2676: else
2677: for ( i = 0; i < n; i++ )
2678: b[i][k] = c[i];
1.3 noro 2679: }
2680: }
2681:
1.24 noro 2682: void Pleqm1(NODE arg,VECT *rp)
1.1 noro 2683: {
2684: MAT m;
2685: VECT vect;
2686: pointer **mat;
2687: Q *v;
2688: Q q;
2689: int **wmat;
2690: int md,i,j,row,col,t,n,status;
2691:
2692: asir_assert(ARG0(arg),O_MAT,"leqm1");
2693: asir_assert(ARG1(arg),O_N,"leqm1");
2694: m = (MAT)ARG0(arg); md = QTOS((Q)ARG1(arg));
2695: row = m->row; col = m->col; mat = m->body;
2696: wmat = (int **)almat(row,col);
2697: for ( i = 0; i < row; i++ )
2698: for ( j = 0; j < col; j++ )
2699: if ( q = (Q)mat[i][j] ) {
2700: t = rem(NM(q),md);
2701: if ( SGN(q) < 0 )
2702: t = (md - t) % md;
2703: wmat[i][j] = t;
2704: } else
2705: wmat[i][j] = 0;
2706: status = gauss_elim_mod1(wmat,row,col,md);
2707: if ( status < 0 )
2708: *rp = 0;
2709: else if ( status > 0 )
2710: *rp = (VECT)ONE;
2711: else {
2712: n = col - 1;
2713: MKVECT(vect,n);
2714: for ( i = 0, v = (Q *)vect->body; i < n; i++ ) {
2715: t = (md-wmat[i][n])%md; STOQ(t,v[i]);
2716: }
2717: *rp = vect;
2718: }
2719: }
2720:
1.24 noro 2721: int gauss_elim_mod1(int **mat,int row,int col,int md)
1.1 noro 2722: {
2723: int i,j,k,inv,a,n;
2724: int *t,*pivot;
2725:
2726: n = col - 1;
2727: for ( j = 0; j < n; j++ ) {
2728: for ( i = j; i < row && !mat[i][j]; i++ );
2729: if ( i == row )
2730: return 1;
2731: if ( i != j ) {
2732: t = mat[i]; mat[i] = mat[j]; mat[j] = t;
2733: }
2734: pivot = mat[j];
2735: inv = invm(pivot[j],md);
2736: for ( k = j; k <= n; k++ )
2737: pivot[k] = dmar(pivot[k],inv,0,md);
2738: for ( i = j+1; i < row; i++ ) {
2739: t = mat[i];
2740: if ( i != j && (a = t[j]) )
2741: for ( k = j, a = md - a; k <= n; k++ )
2742: t[k] = dmar(pivot[k],a,t[k],md);
2743: }
2744: }
2745: for ( i = n; i < row && !mat[i][n]; i++ );
2746: if ( i == row ) {
2747: for ( j = n-1; j >= 0; j-- ) {
2748: for ( i = j-1, a = (md-mat[j][n])%md; i >= 0; i-- ) {
2749: mat[i][n] = dmar(mat[i][j],a,mat[i][n],md);
2750: mat[i][j] = 0;
2751: }
2752: }
2753: return 0;
2754: } else
2755: return -1;
2756: }
2757:
1.24 noro 2758: void Pgeninvm(NODE arg,LIST *rp)
1.1 noro 2759: {
2760: MAT m;
2761: pointer **mat;
2762: Q **tmat;
2763: Q q;
2764: unsigned int **wmat;
2765: int md,i,j,row,col,t,status;
2766: MAT mat1,mat2;
2767: NODE node1,node2;
2768:
2769: asir_assert(ARG0(arg),O_MAT,"leqm1");
2770: asir_assert(ARG1(arg),O_N,"leqm1");
2771: m = (MAT)ARG0(arg); md = QTOS((Q)ARG1(arg));
2772: row = m->row; col = m->col; mat = m->body;
2773: wmat = (unsigned int **)almat(row,col+row);
2774: for ( i = 0; i < row; i++ ) {
2775: bzero((char *)wmat[i],(col+row)*sizeof(int));
2776: for ( j = 0; j < col; j++ )
2777: if ( q = (Q)mat[i][j] ) {
2778: t = rem(NM(q),md);
2779: if ( SGN(q) < 0 )
2780: t = (md - t) % md;
2781: wmat[i][j] = t;
2782: }
2783: wmat[i][col+i] = 1;
2784: }
2785: status = gauss_elim_geninv_mod(wmat,row,col,md);
2786: if ( status > 0 )
2787: *rp = 0;
2788: else {
2789: MKMAT(mat1,col,row); MKMAT(mat2,row-col,row);
2790: for ( i = 0, tmat = (Q **)mat1->body; i < col; i++ )
2791: for ( j = 0; j < row; j++ )
1.24 noro 2792: UTOQ(wmat[i][j+col],tmat[i][j]);
1.1 noro 2793: for ( tmat = (Q **)mat2->body; i < row; i++ )
2794: for ( j = 0; j < row; j++ )
1.24 noro 2795: UTOQ(wmat[i][j+col],tmat[i-col][j]);
1.1 noro 2796: MKNODE(node2,mat2,0); MKNODE(node1,mat1,node2); MKLIST(*rp,node1);
2797: }
2798: }
2799:
1.24 noro 2800: int gauss_elim_geninv_mod(unsigned int **mat,int row,int col,int md)
1.1 noro 2801: {
2802: int i,j,k,inv,a,n,m;
2803: unsigned int *t,*pivot;
2804:
2805: n = col; m = row+col;
2806: for ( j = 0; j < n; j++ ) {
2807: for ( i = j; i < row && !mat[i][j]; i++ );
2808: if ( i == row )
2809: return 1;
2810: if ( i != j ) {
2811: t = mat[i]; mat[i] = mat[j]; mat[j] = t;
2812: }
2813: pivot = mat[j];
2814: inv = invm(pivot[j],md);
2815: for ( k = j; k < m; k++ )
2816: pivot[k] = dmar(pivot[k],inv,0,md);
2817: for ( i = j+1; i < row; i++ ) {
2818: t = mat[i];
2819: if ( a = t[j] )
2820: for ( k = j, a = md - a; k < m; k++ )
2821: t[k] = dmar(pivot[k],a,t[k],md);
2822: }
2823: }
2824: for ( j = n-1; j >= 0; j-- ) {
2825: pivot = mat[j];
2826: for ( i = j-1; i >= 0; i-- ) {
2827: t = mat[i];
2828: if ( a = t[j] )
2829: for ( k = j, a = md - a; k < m; k++ )
2830: t[k] = dmar(pivot[k],a,t[k],md);
2831: }
2832: }
2833: return 0;
2834: }
2835:
1.24 noro 2836: void Psolve_by_lu_gfmmat(NODE arg,VECT *rp)
1.1 noro 2837: {
2838: GFMMAT lu;
2839: Q *perm,*rhs,*v;
2840: int n,i;
2841: unsigned int md;
2842: unsigned int *b,*sol;
2843: VECT r;
2844:
2845: lu = (GFMMAT)ARG0(arg);
2846: perm = (Q *)BDY((VECT)ARG1(arg));
2847: rhs = (Q *)BDY((VECT)ARG2(arg));
2848: md = (unsigned int)QTOS((Q)ARG3(arg));
2849: n = lu->col;
2850: b = (unsigned int *)MALLOC_ATOMIC(n*sizeof(int));
2851: sol = (unsigned int *)MALLOC_ATOMIC(n*sizeof(int));
2852: for ( i = 0; i < n; i++ )
2853: b[i] = QTOS(rhs[QTOS(perm[i])]);
2854: solve_by_lu_gfmmat(lu,md,b,sol);
2855: MKVECT(r,n);
2856: for ( i = 0, v = (Q *)r->body; i < n; i++ )
1.24 noro 2857: UTOQ(sol[i],v[i]);
1.1 noro 2858: *rp = r;
2859: }
2860:
1.24 noro 2861: void solve_by_lu_gfmmat(GFMMAT lu,unsigned int md,
2862: unsigned int *b,unsigned int *x)
1.1 noro 2863: {
2864: int n;
2865: unsigned int **a;
2866: unsigned int *y;
2867: int i,j;
2868: unsigned int t,m;
2869:
2870: n = lu->col;
2871: a = lu->body;
2872: y = (unsigned int *)MALLOC_ATOMIC(n*sizeof(int));
2873: /* solve Ly=b */
2874: for ( i = 0; i < n; i++ ) {
2875: for ( t = b[i], j = 0; j < i; j++ )
2876: if ( a[i][j] ) {
2877: m = md - a[i][j];
2878: DMAR(m,y[j],t,md,t)
2879: }
2880: y[i] = t;
2881: }
2882: /* solve Ux=y */
2883: for ( i = n-1; i >= 0; i-- ) {
2884: for ( t = y[i], j =i+1; j < n; j++ )
2885: if ( a[i][j] ) {
2886: m = md - a[i][j];
2887: DMAR(m,x[j],t,md,t)
2888: }
2889: /* a[i][i] = 1/U[i][i] */
2890: DMAR(t,a[i][i],0,md,x[i])
2891: }
2892: }
2893:
1.53 noro 2894: void Plu_mat(NODE arg,LIST *rp)
2895: {
2896: MAT m,lu;
2897: Q dn;
2898: Q *v;
2899: int n,i;
2900: int *iperm;
2901: VECT perm;
2902: NODE n0;
2903:
2904: asir_assert(ARG0(arg),O_MAT,"lu_mat");
2905: m = (MAT)ARG0(arg);
2906: n = m->row;
2907: MKMAT(lu,n,n);
2908: lu_dec_cr(m,lu,&dn,&iperm);
2909: MKVECT(perm,n);
2910: for ( i = 0, v = (Q *)perm->body; i < n; i++ )
2911: STOQ(iperm[i],v[i]);
2912: n0 = mknode(3,lu,dn,perm);
2913: MKLIST(*rp,n0);
2914: }
2915:
1.24 noro 2916: void Plu_gfmmat(NODE arg,LIST *rp)
1.1 noro 2917: {
2918: MAT m;
2919: GFMMAT mm;
2920: unsigned int md;
2921: int i,row,col,status;
2922: int *iperm;
2923: Q *v;
2924: VECT perm;
2925: NODE n0;
2926:
1.53 noro 2927: asir_assert(ARG0(arg),O_MAT,"lu_gfmmat");
2928: asir_assert(ARG1(arg),O_N,"lu_gfmmat");
1.1 noro 2929: m = (MAT)ARG0(arg); md = (unsigned int)QTOS((Q)ARG1(arg));
2930: mat_to_gfmmat(m,md,&mm);
2931: row = m->row;
2932: col = m->col;
2933: iperm = (int *)MALLOC_ATOMIC(row*sizeof(int));
2934: status = lu_gfmmat(mm,md,iperm);
2935: if ( !status )
2936: n0 = 0;
2937: else {
2938: MKVECT(perm,row);
2939: for ( i = 0, v = (Q *)perm->body; i < row; i++ )
2940: STOQ(iperm[i],v[i]);
2941: n0 = mknode(2,mm,perm);
2942: }
2943: MKLIST(*rp,n0);
2944: }
2945:
1.24 noro 2946: void Pmat_to_gfmmat(NODE arg,GFMMAT *rp)
1.1 noro 2947: {
2948: MAT m;
2949: unsigned int md;
2950:
2951: asir_assert(ARG0(arg),O_MAT,"mat_to_gfmmat");
2952: asir_assert(ARG1(arg),O_N,"mat_to_gfmmat");
2953: m = (MAT)ARG0(arg); md = (unsigned int)QTOS((Q)ARG1(arg));
2954: mat_to_gfmmat(m,md,rp);
2955: }
2956:
1.24 noro 2957: void mat_to_gfmmat(MAT m,unsigned int md,GFMMAT *rp)
1.1 noro 2958: {
2959: unsigned int **wmat;
2960: unsigned int t;
2961: Q **mat;
2962: Q q;
2963: int i,j,row,col;
2964:
2965: row = m->row; col = m->col; mat = (Q **)m->body;
2966: wmat = (unsigned int **)almat(row,col);
2967: for ( i = 0; i < row; i++ ) {
2968: bzero((char *)wmat[i],col*sizeof(unsigned int));
2969: for ( j = 0; j < col; j++ )
2970: if ( q = mat[i][j] ) {
2971: t = (unsigned int)rem(NM(q),md);
2972: if ( SGN(q) < 0 )
2973: t = (md - t) % md;
2974: wmat[i][j] = t;
2975: }
2976: }
2977: TOGFMMAT(row,col,wmat,*rp);
2978: }
2979:
1.27 noro 2980: void Pgeninvm_swap(arg,rp)
2981: NODE arg;
2982: LIST *rp;
1.1 noro 2983: {
2984: MAT m;
2985: pointer **mat;
2986: Q **tmat;
2987: Q *tvect;
2988: Q q;
2989: unsigned int **wmat,**invmat;
2990: int *index;
2991: unsigned int t,md;
2992: int i,j,row,col,status;
2993: MAT mat1;
2994: VECT vect1;
2995: NODE node1,node2;
2996:
2997: asir_assert(ARG0(arg),O_MAT,"geninvm_swap");
2998: asir_assert(ARG1(arg),O_N,"geninvm_swap");
2999: m = (MAT)ARG0(arg); md = QTOS((Q)ARG1(arg));
3000: row = m->row; col = m->col; mat = m->body;
3001: wmat = (unsigned int **)almat(row,col+row);
3002: for ( i = 0; i < row; i++ ) {
3003: bzero((char *)wmat[i],(col+row)*sizeof(int));
3004: for ( j = 0; j < col; j++ )
3005: if ( q = (Q)mat[i][j] ) {
3006: t = (unsigned int)rem(NM(q),md);
3007: if ( SGN(q) < 0 )
3008: t = (md - t) % md;
3009: wmat[i][j] = t;
3010: }
3011: wmat[i][col+i] = 1;
3012: }
3013: status = gauss_elim_geninv_mod_swap(wmat,row,col,md,&invmat,&index);
3014: if ( status > 0 )
3015: *rp = 0;
3016: else {
3017: MKMAT(mat1,col,col);
3018: for ( i = 0, tmat = (Q **)mat1->body; i < col; i++ )
3019: for ( j = 0; j < col; j++ )
3020: UTOQ(invmat[i][j],tmat[i][j]);
3021: MKVECT(vect1,row);
3022: for ( i = 0, tvect = (Q *)vect1->body; i < row; i++ )
3023: STOQ(index[i],tvect[i]);
3024: MKNODE(node2,vect1,0); MKNODE(node1,mat1,node2); MKLIST(*rp,node1);
3025: }
3026: }
3027:
1.27 noro 3028: gauss_elim_geninv_mod_swap(mat,row,col,md,invmatp,indexp)
3029: unsigned int **mat;
3030: int row,col;
3031: unsigned int md;
3032: unsigned int ***invmatp;
3033: int **indexp;
1.1 noro 3034: {
3035: int i,j,k,inv,a,n,m;
3036: unsigned int *t,*pivot,*s;
3037: int *index;
3038: unsigned int **invmat;
3039:
3040: n = col; m = row+col;
3041: *indexp = index = (int *)MALLOC_ATOMIC(row*sizeof(int));
3042: for ( i = 0; i < row; i++ )
3043: index[i] = i;
3044: for ( j = 0; j < n; j++ ) {
3045: for ( i = j; i < row && !mat[i][j]; i++ );
3046: if ( i == row ) {
3047: *indexp = 0; *invmatp = 0; return 1;
3048: }
3049: if ( i != j ) {
3050: t = mat[i]; mat[i] = mat[j]; mat[j] = t;
3051: k = index[i]; index[i] = index[j]; index[j] = k;
3052: }
3053: pivot = mat[j];
3054: inv = (unsigned int)invm(pivot[j],md);
3055: for ( k = j; k < m; k++ )
3056: if ( pivot[k] )
3057: pivot[k] = (unsigned int)dmar(pivot[k],inv,0,md);
3058: for ( i = j+1; i < row; i++ ) {
3059: t = mat[i];
3060: if ( a = t[j] )
3061: for ( k = j, a = md - a; k < m; k++ )
3062: if ( pivot[k] )
3063: t[k] = dmar(pivot[k],a,t[k],md);
3064: }
3065: }
3066: for ( j = n-1; j >= 0; j-- ) {
3067: pivot = mat[j];
3068: for ( i = j-1; i >= 0; i-- ) {
3069: t = mat[i];
3070: if ( a = t[j] )
3071: for ( k = j, a = md - a; k < m; k++ )
3072: if ( pivot[k] )
3073: t[k] = dmar(pivot[k],a,t[k],md);
3074: }
3075: }
3076: *invmatp = invmat = (unsigned int **)almat(col,col);
1.27 noro 3077: for ( i = 0; i < col; i++ )
3078: for ( j = 0, s = invmat[i], t = mat[i]; j < col; j++ )
3079: s[j] = t[col+index[j]];
3080: return 0;
3081: }
3082:
3083: void Pgeninv_sf_swap(NODE arg,LIST *rp)
3084: {
3085: MAT m;
3086: GFS **mat,**tmat;
3087: Q *tvect;
3088: GFS q;
3089: int **wmat,**invmat;
3090: int *index;
3091: unsigned int t;
3092: int i,j,row,col,status;
3093: MAT mat1;
3094: VECT vect1;
3095: NODE node1,node2;
3096:
3097: asir_assert(ARG0(arg),O_MAT,"geninv_sf_swap");
3098: m = (MAT)ARG0(arg);
3099: row = m->row; col = m->col; mat = (GFS **)m->body;
3100: wmat = (int **)almat(row,col+row);
3101: for ( i = 0; i < row; i++ ) {
3102: bzero((char *)wmat[i],(col+row)*sizeof(int));
3103: for ( j = 0; j < col; j++ )
3104: if ( q = (GFS)mat[i][j] )
3105: wmat[i][j] = FTOIF(CONT(q));
3106: wmat[i][col+i] = _onesf();
3107: }
3108: status = gauss_elim_geninv_sf_swap(wmat,row,col,&invmat,&index);
3109: if ( status > 0 )
3110: *rp = 0;
3111: else {
3112: MKMAT(mat1,col,col);
3113: for ( i = 0, tmat = (GFS **)mat1->body; i < col; i++ )
3114: for ( j = 0; j < col; j++ )
3115: if ( t = invmat[i][j] ) {
3116: MKGFS(IFTOF(t),tmat[i][j]);
3117: }
3118: MKVECT(vect1,row);
3119: for ( i = 0, tvect = (Q *)vect1->body; i < row; i++ )
3120: STOQ(index[i],tvect[i]);
3121: MKNODE(node2,vect1,0); MKNODE(node1,mat1,node2); MKLIST(*rp,node1);
3122: }
3123: }
3124:
3125: int gauss_elim_geninv_sf_swap(int **mat,int row,int col,
3126: int ***invmatp,int **indexp)
3127: {
3128: int i,j,k,inv,a,n,m,u;
3129: int *t,*pivot,*s;
3130: int *index;
3131: int **invmat;
3132:
3133: n = col; m = row+col;
3134: *indexp = index = (int *)MALLOC_ATOMIC(row*sizeof(int));
3135: for ( i = 0; i < row; i++ )
3136: index[i] = i;
3137: for ( j = 0; j < n; j++ ) {
3138: for ( i = j; i < row && !mat[i][j]; i++ );
3139: if ( i == row ) {
3140: *indexp = 0; *invmatp = 0; return 1;
3141: }
3142: if ( i != j ) {
3143: t = mat[i]; mat[i] = mat[j]; mat[j] = t;
3144: k = index[i]; index[i] = index[j]; index[j] = k;
3145: }
3146: pivot = mat[j];
3147: inv = _invsf(pivot[j]);
3148: for ( k = j; k < m; k++ )
3149: if ( pivot[k] )
3150: pivot[k] = _mulsf(pivot[k],inv);
3151: for ( i = j+1; i < row; i++ ) {
3152: t = mat[i];
3153: if ( a = t[j] )
3154: for ( k = j, a = _chsgnsf(a); k < m; k++ )
3155: if ( pivot[k] ) {
3156: u = _mulsf(pivot[k],a);
3157: t[k] = _addsf(u,t[k]);
3158: }
3159: }
3160: }
3161: for ( j = n-1; j >= 0; j-- ) {
3162: pivot = mat[j];
3163: for ( i = j-1; i >= 0; i-- ) {
3164: t = mat[i];
3165: if ( a = t[j] )
3166: for ( k = j, a = _chsgnsf(a); k < m; k++ )
3167: if ( pivot[k] ) {
3168: u = _mulsf(pivot[k],a);
3169: t[k] = _addsf(u,t[k]);
3170: }
3171: }
3172: }
3173: *invmatp = invmat = (int **)almat(col,col);
1.1 noro 3174: for ( i = 0; i < col; i++ )
3175: for ( j = 0, s = invmat[i], t = mat[i]; j < col; j++ )
3176: s[j] = t[col+index[j]];
3177: return 0;
3178: }
3179:
3180: void _addn(N,N,N);
3181: int _subn(N,N,N);
3182: void _muln(N,N,N);
3183:
1.24 noro 3184: void inner_product_int(Q *a,Q *b,int n,Q *r)
1.1 noro 3185: {
3186: int la,lb,i;
3187: int sgn,sgn1;
3188: N wm,wma,sum,t;
3189:
3190: for ( la = lb = 0, i = 0; i < n; i++ ) {
3191: if ( a[i] )
3192: if ( DN(a[i]) )
3193: error("inner_product_int : invalid argument");
3194: else
3195: la = MAX(PL(NM(a[i])),la);
3196: if ( b[i] )
3197: if ( DN(b[i]) )
3198: error("inner_product_int : invalid argument");
3199: else
3200: lb = MAX(PL(NM(b[i])),lb);
3201: }
3202: sgn = 0;
3203: sum= NALLOC(la+lb+2);
3204: bzero((char *)sum,(la+lb+3)*sizeof(unsigned int));
3205: wm = NALLOC(la+lb+2);
3206: wma = NALLOC(la+lb+2);
3207: for ( i = 0; i < n; i++ ) {
3208: if ( !a[i] || !b[i] )
3209: continue;
3210: _muln(NM(a[i]),NM(b[i]),wm);
3211: sgn1 = SGN(a[i])*SGN(b[i]);
3212: if ( !sgn ) {
3213: sgn = sgn1;
3214: t = wm; wm = sum; sum = t;
3215: } else if ( sgn == sgn1 ) {
3216: _addn(sum,wm,wma);
3217: if ( !PL(wma) )
3218: sgn = 0;
3219: t = wma; wma = sum; sum = t;
3220: } else {
3221: /* sgn*sum+sgn1*wm = sgn*(sum-wm) */
3222: sgn *= _subn(sum,wm,wma);
3223: t = wma; wma = sum; sum = t;
3224: }
3225: }
1.61 noro 3226: GCFREE(wm);
3227: GCFREE(wma);
1.1 noro 3228: if ( !sgn ) {
1.61 noro 3229: GCFREE(sum);
1.1 noro 3230: *r = 0;
3231: } else
3232: NTOQ(sum,sgn,*r);
3233: }
3234:
1.3 noro 3235: /* (k,l) element of a*b where a: .x n matrix, b: n x . integer matrix */
3236:
1.24 noro 3237: void inner_product_mat_int_mod(Q **a,int **b,int n,int k,int l,Q *r)
1.3 noro 3238: {
3239: int la,lb,i;
3240: int sgn,sgn1;
3241: N wm,wma,sum,t;
3242: Q aki;
3243: int bil,bilsgn;
3244: struct oN tn;
3245:
3246: for ( la = 0, i = 0; i < n; i++ ) {
3247: if ( aki = a[k][i] )
3248: if ( DN(aki) )
3249: error("inner_product_int : invalid argument");
3250: else
3251: la = MAX(PL(NM(aki)),la);
3252: }
3253: lb = 1;
3254: sgn = 0;
3255: sum= NALLOC(la+lb+2);
3256: bzero((char *)sum,(la+lb+3)*sizeof(unsigned int));
3257: wm = NALLOC(la+lb+2);
3258: wma = NALLOC(la+lb+2);
3259: for ( i = 0; i < n; i++ ) {
3260: if ( !(aki = a[k][i]) || !(bil = b[i][l]) )
3261: continue;
3262: tn.p = 1;
3263: if ( bil > 0 ) {
3264: tn.b[0] = bil; bilsgn = 1;
3265: } else {
3266: tn.b[0] = -bil; bilsgn = -1;
3267: }
3268: _muln(NM(aki),&tn,wm);
3269: sgn1 = SGN(aki)*bilsgn;
3270: if ( !sgn ) {
3271: sgn = sgn1;
3272: t = wm; wm = sum; sum = t;
3273: } else if ( sgn == sgn1 ) {
3274: _addn(sum,wm,wma);
3275: if ( !PL(wma) )
3276: sgn = 0;
3277: t = wma; wma = sum; sum = t;
3278: } else {
3279: /* sgn*sum+sgn1*wm = sgn*(sum-wm) */
3280: sgn *= _subn(sum,wm,wma);
3281: t = wma; wma = sum; sum = t;
3282: }
3283: }
1.61 noro 3284: GCFREE(wm);
3285: GCFREE(wma);
1.3 noro 3286: if ( !sgn ) {
1.61 noro 3287: GCFREE(sum);
1.3 noro 3288: *r = 0;
3289: } else
3290: NTOQ(sum,sgn,*r);
3291: }
3292:
1.24 noro 3293: void Pmul_mat_vect_int(NODE arg,VECT *rp)
1.1 noro 3294: {
3295: MAT mat;
3296: VECT vect,r;
3297: int row,col,i;
3298:
3299: mat = (MAT)ARG0(arg);
3300: vect = (VECT)ARG1(arg);
3301: row = mat->row;
3302: col = mat->col;
3303: MKVECT(r,row);
1.24 noro 3304: for ( i = 0; i < row; i++ ) {
3305: inner_product_int((Q *)mat->body[i],(Q *)vect->body,col,(Q *)&r->body[i]);
3306: }
1.1 noro 3307: *rp = r;
3308: }
3309:
1.24 noro 3310: void Pnbpoly_up2(NODE arg,GF2N *rp)
1.1 noro 3311: {
3312: int m,type,ret;
3313: UP2 r;
3314:
3315: m = QTOS((Q)ARG0(arg));
3316: type = QTOS((Q)ARG1(arg));
3317: ret = generate_ONB_polynomial(&r,m,type);
3318: if ( ret == 0 )
3319: MKGF2N(r,*rp);
3320: else
3321: *rp = 0;
3322: }
3323:
1.24 noro 3324: void Px962_irredpoly_up2(NODE arg,GF2N *rp)
1.1 noro 3325: {
1.24 noro 3326: int m,ret,w;
1.1 noro 3327: GF2N prev;
3328: UP2 r;
3329:
3330: m = QTOS((Q)ARG0(arg));
3331: prev = (GF2N)ARG1(arg);
3332: if ( !prev ) {
3333: w = (m>>5)+1; NEWUP2(r,w); r->w = 0;
3334: bzero((char *)r->b,w*sizeof(unsigned int));
3335: } else {
3336: r = prev->body;
3337: if ( degup2(r) != m ) {
3338: w = (m>>5)+1; NEWUP2(r,w); r->w = 0;
3339: bzero((char *)r->b,w*sizeof(unsigned int));
3340: }
3341: }
1.24 noro 3342: ret = _generate_irreducible_polynomial(r,m);
1.1 noro 3343: if ( ret == 0 )
3344: MKGF2N(r,*rp);
3345: else
3346: *rp = 0;
3347: }
3348:
1.24 noro 3349: void Pirredpoly_up2(NODE arg,GF2N *rp)
1.1 noro 3350: {
1.24 noro 3351: int m,ret,w;
1.1 noro 3352: GF2N prev;
3353: UP2 r;
3354:
3355: m = QTOS((Q)ARG0(arg));
3356: prev = (GF2N)ARG1(arg);
3357: if ( !prev ) {
3358: w = (m>>5)+1; NEWUP2(r,w); r->w = 0;
3359: bzero((char *)r->b,w*sizeof(unsigned int));
3360: } else {
3361: r = prev->body;
3362: if ( degup2(r) != m ) {
3363: w = (m>>5)+1; NEWUP2(r,w); r->w = 0;
3364: bzero((char *)r->b,w*sizeof(unsigned int));
3365: }
3366: }
1.24 noro 3367: ret = _generate_good_irreducible_polynomial(r,m);
1.1 noro 3368: if ( ret == 0 )
3369: MKGF2N(r,*rp);
3370: else
3371: *rp = 0;
3372: }
3373:
1.26 noro 3374: void Pmat_swap_row_destructive(NODE arg, MAT *m)
3375: {
3376: int i1,i2;
3377: pointer *t;
3378: MAT mat;
3379:
3380: asir_assert(ARG0(arg),O_MAT,"mat_swap_row_destructive");
3381: asir_assert(ARG1(arg),O_N,"mat_swap_row_destructive");
3382: asir_assert(ARG2(arg),O_N,"mat_swap_row_destructive");
3383: mat = (MAT)ARG0(arg);
3384: i1 = QTOS((Q)ARG1(arg));
3385: i2 = QTOS((Q)ARG2(arg));
3386: if ( i1 < 0 || i2 < 0 || i1 >= mat->row || i2 >= mat->row )
3387: error("mat_swap_row_destructive : Out of range");
3388: t = mat->body[i1];
3389: mat->body[i1] = mat->body[i2];
3390: mat->body[i2] = t;
3391: *m = mat;
3392: }
3393:
3394: void Pmat_swap_col_destructive(NODE arg, MAT *m)
3395: {
3396: int j1,j2,i,n;
3397: pointer *mi;
3398: pointer t;
3399: MAT mat;
3400:
3401: asir_assert(ARG0(arg),O_MAT,"mat_swap_col_destructive");
3402: asir_assert(ARG1(arg),O_N,"mat_swap_col_destructive");
3403: asir_assert(ARG2(arg),O_N,"mat_swap_col_destructive");
3404: mat = (MAT)ARG0(arg);
3405: j1 = QTOS((Q)ARG1(arg));
3406: j2 = QTOS((Q)ARG2(arg));
3407: if ( j1 < 0 || j2 < 0 || j1 >= mat->col || j2 >= mat->col )
3408: error("mat_swap_col_destructive : Out of range");
3409: n = mat->row;
3410: for ( i = 0; i < n; i++ ) {
3411: mi = mat->body[i];
3412: t = mi[j1]; mi[j1] = mi[j2]; mi[j2] = t;
3413: }
3414: *m = mat;
3415: }
1.1 noro 3416: /*
3417: * f = type 'type' normal polynomial of degree m if exists
3418: * IEEE P1363 A.7.2
3419: *
3420: * return value : 0 --- exists
3421: * 1 --- does not exist
3422: * -1 --- failure (memory allocation error)
3423: */
3424:
3425: int generate_ONB_polynomial(UP2 *rp,int m,int type)
3426: {
3427: int i,r;
3428: int w;
3429: UP2 f,f0,f1,f2,t;
3430:
3431: w = (m>>5)+1;
3432: switch ( type ) {
3433: case 1:
3434: if ( !TypeT_NB_check(m,1) ) return 1;
3435: NEWUP2(f,w); *rp = f; f->w = w;
3436: /* set all the bits */
3437: for ( i = 0; i < w; i++ )
3438: f->b[i] = 0xffffffff;
3439: /* mask the top word if necessary */
3440: if ( r = (m+1)&31 )
3441: f->b[w-1] &= (1<<r)-1;
3442: return 0;
3443: break;
3444: case 2:
3445: if ( !TypeT_NB_check(m,2) ) return 1;
3446: NEWUP2(f,w); *rp = f;
3447: W_NEWUP2(f0,w);
3448: W_NEWUP2(f1,w);
3449: W_NEWUP2(f2,w);
3450:
3451: /* recursion for genrating Type II normal polynomial */
3452:
3453: /* f0 = 1, f1 = t+1 */
3454: f0->w = 1; f0->b[0] = 1;
3455: f1->w = 1; f1->b[0] = 3;
3456: for ( i = 2; i <= m; i++ ) {
3457: /* f2 = t*f1+f0 */
3458: _bshiftup2(f1,-1,f2);
3459: _addup2_destructive(f2,f0);
3460: /* cyclic change of the variables */
3461: t = f0; f0 = f1; f1 = f2; f2 = t;
3462: }
3463: _copyup2(f1,f);
3464: return 0;
3465: break;
3466: default:
3467: return -1;
3468: break;
3469: }
3470: }
3471:
3472: /*
3473: * f = an irreducible trinomial or pentanomial of degree d 'after' f
3474: * return value : 0 --- exists
3475: * 1 --- does not exist (exhaustion)
3476: */
3477:
3478: int _generate_irreducible_polynomial(UP2 f,int d)
3479: {
3480: int ret,i,j,k,nz,i0,j0,k0;
3481: int w;
3482: unsigned int *fd;
3483:
3484: /*
3485: * if f = x^d+x^i+1 then i0 <- i, j0 <- 0, k0 <-0.
3486: * if f = x^d+x^k+x^j+x^i+1 (k>j>i) then i0 <- i, j0 <- j, k0 <-k.
3487: * otherwise i0,j0,k0 is set to 0.
3488: */
3489:
3490: fd = f->b;
3491: w = (d>>5)+1;
3492: if ( f->w && (d==degup2(f)) ) {
3493: for ( nz = 0, i = d; i >= 0; i-- )
3494: if ( fd[i>>5]&(1<<(i&31)) ) nz++;
3495: switch ( nz ) {
3496: case 3:
3497: for ( i0 = 1; !(fd[i0>>5]&(1<<(i0&31))) ; i0++ );
3498: /* reset i0-th bit */
3499: fd[i0>>5] &= ~(1<<(i0&31));
3500: j0 = k0 = 0;
3501: break;
3502: case 5:
3503: for ( i0 = 1; !(fd[i0>>5]&(1<<(i0&31))) ; i0++ );
3504: /* reset i0-th bit */
3505: fd[i0>>5] &= ~(1<<(i0&31));
3506: for ( j0 = i0+1; !(fd[j0>>5]&(1<<(j0&31))) ; j0++ );
3507: /* reset j0-th bit */
3508: fd[j0>>5] &= ~(1<<(j0&31));
3509: for ( k0 = j0+1; !(fd[k0>>5]&(1<<(k0&31))) ; k0++ );
3510: /* reset k0-th bit */
3511: fd[k0>>5] &= ~(1<<(k0&31));
3512: break;
3513: default:
3514: f->w = 0; break;
3515: }
3516: } else
3517: f->w = 0;
3518:
3519: if ( !f->w ) {
3520: fd = f->b;
3521: f->w = w; fd[0] |= 1; fd[d>>5] |= (1<<(d&31));
3522: i0 = j0 = k0 = 0;
3523: }
3524: /* if j0 > 0 then f is already a pentanomial */
3525: if ( j0 > 0 ) goto PENTA;
3526:
3527: /* searching for an irreducible trinomial */
3528:
3529: for ( i = 1; 2*i <= d; i++ ) {
3530: /* skip the polynomials 'before' f */
3531: if ( i < i0 ) continue;
3532: if ( i == i0 ) { i0 = 0; continue; }
3533: /* set i-th bit */
3534: fd[i>>5] |= (1<<(i&31));
3535: ret = irredcheck_dddup2(f);
3536: if ( ret == 1 ) return 0;
3537: /* reset i-th bit */
3538: fd[i>>5] &= ~(1<<(i&31));
3539: }
3540:
3541: /* searching for an irreducible pentanomial */
3542: PENTA:
3543: for ( i = 1; i < d; i++ ) {
3544: /* skip the polynomials 'before' f */
3545: if ( i < i0 ) continue;
3546: if ( i == i0 ) i0 = 0;
3547: /* set i-th bit */
3548: fd[i>>5] |= (1<<(i&31));
3549: for ( j = i+1; j < d; j++ ) {
3550: /* skip the polynomials 'before' f */
3551: if ( j < j0 ) continue;
3552: if ( j == j0 ) j0 = 0;
3553: /* set j-th bit */
3554: fd[j>>5] |= (1<<(j&31));
3555: for ( k = j+1; k < d; k++ ) {
3556: /* skip the polynomials 'before' f */
3557: if ( k < k0 ) continue;
3558: else if ( k == k0 ) { k0 = 0; continue; }
3559: /* set k-th bit */
3560: fd[k>>5] |= (1<<(k&31));
3561: ret = irredcheck_dddup2(f);
3562: if ( ret == 1 ) return 0;
3563: /* reset k-th bit */
3564: fd[k>>5] &= ~(1<<(k&31));
3565: }
3566: /* reset j-th bit */
3567: fd[j>>5] &= ~(1<<(j&31));
3568: }
3569: /* reset i-th bit */
3570: fd[i>>5] &= ~(1<<(i&31));
3571: }
3572: /* exhausted */
3573: return 1;
3574: }
3575:
3576: /*
3577: * f = an irreducible trinomial or pentanomial of degree d 'after' f
3578: *
3579: * searching strategy:
3580: * trinomial x^d+x^i+1:
3581: * i is as small as possible.
3582: * trinomial x^d+x^i+x^j+x^k+1:
3583: * i is as small as possible.
3584: * For such i, j is as small as possible.
3585: * For such i and j, 'k' is as small as possible.
3586: *
3587: * return value : 0 --- exists
3588: * 1 --- does not exist (exhaustion)
3589: */
3590:
3591: int _generate_good_irreducible_polynomial(UP2 f,int d)
3592: {
3593: int ret,i,j,k,nz,i0,j0,k0;
3594: int w;
3595: unsigned int *fd;
3596:
3597: /*
3598: * if f = x^d+x^i+1 then i0 <- i, j0 <- 0, k0 <-0.
3599: * if f = x^d+x^k+x^j+x^i+1 (k>j>i) then i0 <- i, j0 <- j, k0 <-k.
3600: * otherwise i0,j0,k0 is set to 0.
3601: */
3602:
3603: fd = f->b;
3604: w = (d>>5)+1;
3605: if ( f->w && (d==degup2(f)) ) {
3606: for ( nz = 0, i = d; i >= 0; i-- )
3607: if ( fd[i>>5]&(1<<(i&31)) ) nz++;
3608: switch ( nz ) {
3609: case 3:
3610: for ( i0 = 1; !(fd[i0>>5]&(1<<(i0&31))) ; i0++ );
3611: /* reset i0-th bit */
3612: fd[i0>>5] &= ~(1<<(i0&31));
3613: j0 = k0 = 0;
3614: break;
3615: case 5:
3616: for ( i0 = 1; !(fd[i0>>5]&(1<<(i0&31))) ; i0++ );
3617: /* reset i0-th bit */
3618: fd[i0>>5] &= ~(1<<(i0&31));
3619: for ( j0 = i0+1; !(fd[j0>>5]&(1<<(j0&31))) ; j0++ );
3620: /* reset j0-th bit */
3621: fd[j0>>5] &= ~(1<<(j0&31));
3622: for ( k0 = j0+1; !(fd[k0>>5]&(1<<(k0&31))) ; k0++ );
3623: /* reset k0-th bit */
3624: fd[k0>>5] &= ~(1<<(k0&31));
3625: break;
3626: default:
3627: f->w = 0; break;
3628: }
3629: } else
3630: f->w = 0;
3631:
3632: if ( !f->w ) {
3633: fd = f->b;
3634: f->w = w; fd[0] |= 1; fd[d>>5] |= (1<<(d&31));
3635: i0 = j0 = k0 = 0;
3636: }
3637: /* if j0 > 0 then f is already a pentanomial */
3638: if ( j0 > 0 ) goto PENTA;
3639:
3640: /* searching for an irreducible trinomial */
3641:
3642: for ( i = 1; 2*i <= d; i++ ) {
3643: /* skip the polynomials 'before' f */
3644: if ( i < i0 ) continue;
3645: if ( i == i0 ) { i0 = 0; continue; }
3646: /* set i-th bit */
3647: fd[i>>5] |= (1<<(i&31));
3648: ret = irredcheck_dddup2(f);
3649: if ( ret == 1 ) return 0;
3650: /* reset i-th bit */
3651: fd[i>>5] &= ~(1<<(i&31));
3652: }
3653:
3654: /* searching for an irreducible pentanomial */
3655: PENTA:
3656: for ( i = 3; i < d; i++ ) {
3657: /* skip the polynomials 'before' f */
3658: if ( i < i0 ) continue;
3659: if ( i == i0 ) i0 = 0;
3660: /* set i-th bit */
3661: fd[i>>5] |= (1<<(i&31));
3662: for ( j = 2; j < i; j++ ) {
3663: /* skip the polynomials 'before' f */
3664: if ( j < j0 ) continue;
3665: if ( j == j0 ) j0 = 0;
3666: /* set j-th bit */
3667: fd[j>>5] |= (1<<(j&31));
3668: for ( k = 1; k < j; k++ ) {
3669: /* skip the polynomials 'before' f */
3670: if ( k < k0 ) continue;
3671: else if ( k == k0 ) { k0 = 0; continue; }
3672: /* set k-th bit */
3673: fd[k>>5] |= (1<<(k&31));
3674: ret = irredcheck_dddup2(f);
3675: if ( ret == 1 ) return 0;
3676: /* reset k-th bit */
3677: fd[k>>5] &= ~(1<<(k&31));
3678: }
3679: /* reset j-th bit */
3680: fd[j>>5] &= ~(1<<(j&31));
3681: }
3682: /* reset i-th bit */
3683: fd[i>>5] &= ~(1<<(i&31));
3684: }
3685: /* exhausted */
3686: return 1;
1.3 noro 3687: }
3688:
1.24 noro 3689: void printqmat(Q **mat,int row,int col)
1.3 noro 3690: {
3691: int i,j;
3692:
3693: for ( i = 0; i < row; i++ ) {
3694: for ( j = 0; j < col; j++ ) {
1.8 noro 3695: printnum((Num)mat[i][j]); printf(" ");
1.3 noro 3696: }
3697: printf("\n");
3698: }
3699: }
3700:
1.24 noro 3701: void printimat(int **mat,int row,int col)
1.3 noro 3702: {
3703: int i,j;
3704:
3705: for ( i = 0; i < row; i++ ) {
3706: for ( j = 0; j < col; j++ ) {
3707: printf("%d ",mat[i][j]);
3708: }
3709: printf("\n");
3710: }
1.36 noro 3711: }
3712:
3713: void Pnd_det(NODE arg,P *rp)
3714: {
1.37 noro 3715: if ( argc(arg) == 1 )
3716: nd_det(0,ARG0(arg),rp);
3717: else
3718: nd_det(QTOS((Q)ARG1(arg)),ARG0(arg),rp);
1.1 noro 3719: }
1.59 ohara 3720:
1.62 ohara 3721: void Pmat_col(NODE arg,VECT *rp)
1.59 ohara 3722: {
3723: int i,j,n;
3724: MAT mat;
3725: VECT vect;
3726:
3727: asir_assert(ARG0(arg),O_MAT,"mat_col");
3728: asir_assert(ARG1(arg),O_N,"mat_col");
3729: mat = (MAT)ARG0(arg);
3730: j = QTOS((Q)ARG1(arg));
3731: if ( j < 0 || j >= mat->col) {
3732: error("mat_col : Out of range");
3733: }
3734: n = mat->row;
3735: MKVECT(vect,n);
3736: for(i=0; i<n; i++) {
3737: BDY(vect)[i] = BDY(mat)[i][j];
3738: }
3739: *rp = vect;
3740: }