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