![[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.70 ! noro 48: * $OpenXM: OpenXM_contrib2/asir2000/builtin/array.c,v 1.69 2015/09/03 23:05:35 noro Exp $
1.6 noro 49: */
1.1 noro 50: #include "ca.h"
51: #include "base.h"
52: #include "parse.h"
53: #include "inline.h"
1.4 noro 54:
1.51 noro 55: #include <sys/types.h>
56: #include <sys/stat.h>
1.58 ohara 57: #if !defined(_MSC_VER)
1.51 noro 58: #include <unistd.h>
1.58 ohara 59: #endif
1.51 noro 60:
1.38 noro 61: #define F4_INTRAT_PERIOD 8
62:
1.4 noro 63: #if 0
1.1 noro 64: #undef DMAR
65: #define DMAR(a1,a2,a3,d,r) (r)=dmar(a1,a2,a3,d);
1.4 noro 66: #endif
1.1 noro 67:
1.11 noro 68: extern int DP_Print; /* XXX */
1.1 noro 69:
1.24 noro 70:
1.1 noro 71: void Pnewvect(), Pnewmat(), Psepvect(), Psize(), Pdet(), Pleqm(), Pleqm1(), Pgeninvm();
1.23 noro 72: void Pinvmat();
1.49 noro 73: void Pnewbytearray(),Pmemoryplot_to_coord();
1.1 noro 74:
1.25 noro 75: void Pgeneric_gauss_elim();
1.1 noro 76: void Pgeneric_gauss_elim_mod();
77:
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.70 ! noro 2362: void red_by_vect_lf(mpz_t *p,mpz_t *r,mpz_t hc,int len)
! 2363: {
! 2364: mpz_set_ui(*p++,0); r++; len--;
! 2365: for ( ; len; len--, r++, p++ )
! 2366: mpz_addmul(*p,*r,hc);
! 2367: }
! 2368:
! 2369:
1.21 noro 2370: extern unsigned int **psca;
2371:
1.24 noro 2372: void reduce_sp_by_red_mod_compress (int *sp,CDP *redmat,int *ind,
2373: int nred,int col,int md)
1.15 noro 2374: {
1.24 noro 2375: int i,len;
1.15 noro 2376: CDP ri;
1.24 noro 2377: unsigned int hc;
1.18 noro 2378: unsigned int *usp;
1.15 noro 2379:
1.18 noro 2380: usp = (unsigned int *)sp;
1.15 noro 2381: /* reduce the spolys by redmat */
2382: for ( i = nred-1; i >= 0; i-- ) {
2383: /* reduce sp by redmat[i] */
1.18 noro 2384: usp[ind[i]] %= md;
2385: if ( hc = usp[ind[i]] ) {
1.15 noro 2386: /* sp = sp-hc*redmat[i] */
2387: hc = md-hc;
2388: ri = redmat[i];
2389: len = ri->len;
1.21 noro 2390: red_by_compress(md,usp,psca[ri->psindex],ri->body,hc,len);
1.4 noro 2391: }
2392: }
1.18 noro 2393: for ( i = 0; i < col; i++ )
1.24 noro 2394: if ( usp[i] >= (unsigned int)md )
1.18 noro 2395: usp[i] %= md;
1.4 noro 2396: }
2397:
2398: #define ONE_STEP2 if ( zzz = *pk ) { DMAR(zzz,a,*tk,md,*tk) } pk++; tk++;
2399:
1.24 noro 2400: int generic_gauss_elim_mod(int **mat0,int row,int col,int md,int *colstat)
1.1 noro 2401: {
1.24 noro 2402: int i,j,k,l,inv,a,rank;
2403: unsigned int *t,*pivot,*pk;
1.18 noro 2404: unsigned int **mat;
1.1 noro 2405:
1.18 noro 2406: mat = (unsigned int **)mat0;
1.1 noro 2407: for ( rank = 0, j = 0; j < col; j++ ) {
1.18 noro 2408: for ( i = rank; i < row; i++ )
2409: mat[i][j] %= md;
2410: for ( i = rank; i < row; i++ )
2411: if ( mat[i][j] )
2412: break;
1.1 noro 2413: if ( i == row ) {
2414: colstat[j] = 0;
2415: continue;
2416: } else
2417: colstat[j] = 1;
2418: if ( i != rank ) {
2419: t = mat[i]; mat[i] = mat[rank]; mat[rank] = t;
2420: }
2421: pivot = mat[rank];
2422: inv = invm(pivot[j],md);
1.4 noro 2423: for ( k = j, pk = pivot+k; k < col; k++, pk++ )
2424: if ( *pk ) {
1.24 noro 2425: if ( *pk >= (unsigned int)md )
1.18 noro 2426: *pk %= md;
1.4 noro 2427: DMAR(*pk,inv,0,md,*pk)
1.1 noro 2428: }
2429: for ( i = rank+1; i < row; i++ ) {
2430: t = mat[i];
1.18 noro 2431: if ( a = t[j] )
2432: red_by_vect(md,t+j,pivot+j,md-a,col-j);
1.1 noro 2433: }
2434: rank++;
2435: }
2436: for ( j = col-1, l = rank-1; j >= 0; j-- )
2437: if ( colstat[j] ) {
2438: pivot = mat[l];
2439: for ( i = 0; i < l; i++ ) {
2440: t = mat[i];
1.18 noro 2441: t[j] %= md;
2442: if ( a = t[j] )
2443: red_by_vect(md,t+j,pivot+j,md-a,col-j);
1.1 noro 2444: }
2445: l--;
1.18 noro 2446: }
2447: for ( j = 0, l = 0; l < rank; j++ )
2448: if ( colstat[j] ) {
2449: t = mat[l];
2450: for ( k = j; k < col; k++ )
1.24 noro 2451: if ( t[k] >= (unsigned int)md )
1.18 noro 2452: t[k] %= md;
2453: l++;
1.32 noro 2454: }
2455: return rank;
2456: }
2457:
1.65 noro 2458: int generic_gauss_elim_mod2(int **mat0,int row,int col,int md,int *colstat,int *rowstat)
2459: {
2460: int i,j,k,l,inv,a,rank;
2461: unsigned int *t,*pivot,*pk;
2462: unsigned int **mat;
2463:
2464: for ( i = 0; i < row; i++ ) rowstat[i] = i;
2465: mat = (unsigned int **)mat0;
2466: for ( rank = 0, j = 0; j < col; j++ ) {
2467: for ( i = rank; i < row; i++ )
2468: mat[i][j] %= md;
2469: for ( i = rank; i < row; i++ )
2470: if ( mat[i][j] )
2471: break;
2472: if ( i == row ) {
2473: colstat[j] = 0;
2474: continue;
2475: } else
2476: colstat[j] = 1;
2477: if ( i != rank ) {
2478: t = mat[i]; mat[i] = mat[rank]; mat[rank] = t;
2479: k = rowstat[i]; rowstat[i] = rowstat[rank]; rowstat[rank] = k;
2480: }
2481: pivot = mat[rank];
2482: inv = invm(pivot[j],md);
2483: for ( k = j, pk = pivot+k; k < col; k++, pk++ )
2484: if ( *pk ) {
2485: if ( *pk >= (unsigned int)md )
2486: *pk %= md;
2487: DMAR(*pk,inv,0,md,*pk)
2488: }
2489: for ( i = rank+1; i < row; i++ ) {
2490: t = mat[i];
2491: if ( a = t[j] )
2492: red_by_vect(md,t+j,pivot+j,md-a,col-j);
2493: }
2494: rank++;
2495: }
2496: for ( j = col-1, l = rank-1; j >= 0; j-- )
2497: if ( colstat[j] ) {
2498: pivot = mat[l];
2499: for ( i = 0; i < l; i++ ) {
2500: t = mat[i];
2501: t[j] %= md;
2502: if ( a = t[j] )
2503: red_by_vect(md,t+j,pivot+j,md-a,col-j);
2504: }
2505: l--;
2506: }
2507: for ( j = 0, l = 0; l < rank; j++ )
2508: if ( colstat[j] ) {
2509: t = mat[l];
2510: for ( k = j; k < col; k++ )
2511: if ( t[k] >= (unsigned int)md )
2512: t[k] %= md;
2513: l++;
2514: }
2515: return rank;
2516: }
2517:
1.69 noro 2518: int indep_rows_mod(int **mat0,int row,int col,int md,int *rowstat)
2519: {
2520: int i,j,k,l,inv,a,rank;
2521: unsigned int *t,*pivot,*pk;
2522: unsigned int **mat;
2523:
2524: for ( i = 0; i < row; i++ ) rowstat[i] = i;
2525: mat = (unsigned int **)mat0;
2526: for ( rank = 0, j = 0; j < col; j++ ) {
2527: for ( i = rank; i < row; i++ )
2528: mat[i][j] %= md;
2529: for ( i = rank; i < row; i++ )
2530: if ( mat[i][j] )
2531: break;
2532: if ( i == row ) continue;
2533: if ( i != rank ) {
2534: t = mat[i]; mat[i] = mat[rank]; mat[rank] = t;
2535: k = rowstat[i]; rowstat[i] = rowstat[rank]; rowstat[rank] = k;
2536: }
2537: pivot = mat[rank];
2538: inv = invm(pivot[j],md);
2539: for ( k = j, pk = pivot+k; k < col; k++, pk++ )
2540: if ( *pk ) {
2541: if ( *pk >= (unsigned int)md )
2542: *pk %= md;
2543: DMAR(*pk,inv,0,md,*pk)
2544: }
2545: for ( i = rank+1; i < row; i++ ) {
2546: t = mat[i];
2547: if ( a = t[j] )
2548: red_by_vect(md,t+j,pivot+j,md-a,col-j);
2549: }
2550: rank++;
2551: }
2552: return rank;
2553: }
2554:
1.32 noro 2555: int generic_gauss_elim_sf(int **mat0,int row,int col,int md,int *colstat)
2556: {
2557: int i,j,k,l,inv,a,rank;
2558: unsigned int *t,*pivot,*pk;
2559: unsigned int **mat;
2560:
2561: mat = (unsigned int **)mat0;
2562: for ( rank = 0, j = 0; j < col; j++ ) {
2563: for ( i = rank; i < row; i++ )
2564: if ( mat[i][j] )
2565: break;
2566: if ( i == row ) {
2567: colstat[j] = 0;
2568: continue;
2569: } else
2570: colstat[j] = 1;
2571: if ( i != rank ) {
2572: t = mat[i]; mat[i] = mat[rank]; mat[rank] = t;
2573: }
2574: pivot = mat[rank];
2575: inv = _invsf(pivot[j]);
2576: for ( k = j, pk = pivot+k; k < col; k++, pk++ )
2577: if ( *pk )
2578: *pk = _mulsf(*pk,inv);
2579: for ( i = rank+1; i < row; i++ ) {
2580: t = mat[i];
2581: if ( a = t[j] )
2582: red_by_vect_sf(md,t+j,pivot+j,_chsgnsf(a),col-j);
2583: }
2584: rank++;
2585: }
2586: for ( j = col-1, l = rank-1; j >= 0; j-- )
2587: if ( colstat[j] ) {
2588: pivot = mat[l];
2589: for ( i = 0; i < l; i++ ) {
2590: t = mat[i];
2591: if ( a = t[j] )
2592: red_by_vect_sf(md,t+j,pivot+j,_chsgnsf(a),col-j);
2593: }
2594: l--;
1.1 noro 2595: }
2596: return rank;
2597: }
2598:
2599: /* LU decomposition; a[i][i] = 1/U[i][i] */
2600:
1.24 noro 2601: int lu_gfmmat(GFMMAT mat,unsigned int md,int *perm)
1.1 noro 2602: {
2603: int row,col;
1.24 noro 2604: int i,j,k;
1.1 noro 2605: unsigned int *t,*pivot;
2606: unsigned int **a;
2607: unsigned int inv,m;
2608:
2609: row = mat->row; col = mat->col;
2610: a = mat->body;
2611: bzero(perm,row*sizeof(int));
2612:
2613: for ( i = 0; i < row; i++ )
2614: perm[i] = i;
2615: for ( k = 0; k < col; k++ ) {
2616: for ( i = k; i < row && !a[i][k]; i++ );
2617: if ( i == row )
2618: return 0;
2619: if ( i != k ) {
2620: j = perm[i]; perm[i] = perm[k]; perm[k] = j;
2621: t = a[i]; a[i] = a[k]; a[k] = t;
2622: }
2623: pivot = a[k];
2624: pivot[k] = inv = invm(pivot[k],md);
2625: for ( i = k+1; i < row; i++ ) {
2626: t = a[i];
2627: if ( m = t[k] ) {
2628: DMAR(inv,m,0,md,t[k])
2629: for ( j = k+1, m = md - t[k]; j < col; j++ )
2630: if ( pivot[j] ) {
1.8 noro 2631: unsigned int tj;
2632:
2633: DMAR(m,pivot[j],t[j],md,tj)
2634: t[j] = tj;
1.1 noro 2635: }
2636: }
2637: }
2638: }
2639: return 1;
2640: }
2641:
1.3 noro 2642: /*
2643: Input
2644: a: a row x col matrix
2645: md : a modulus
2646:
2647: Output:
2648: return : d = the rank of mat
2649: a[0..(d-1)][0..(d-1)] : LU decomposition (a[i][i] = 1/U[i][i])
2650: rinfo: array of length row
2651: cinfo: array of length col
2652: i-th row in new a <-> rinfo[i]-th row in old a
2653: cinfo[j]=1 <=> j-th column is contained in the LU decomp.
2654: */
2655:
1.24 noro 2656: int find_lhs_and_lu_mod(unsigned int **a,int row,int col,
2657: unsigned int md,int **rinfo,int **cinfo)
1.3 noro 2658: {
1.24 noro 2659: int i,j,k,d;
1.3 noro 2660: int *rp,*cp;
2661: unsigned int *t,*pivot;
2662: unsigned int inv,m;
2663:
2664: *rinfo = rp = (int *)MALLOC_ATOMIC(row*sizeof(int));
2665: *cinfo = cp = (int *)MALLOC_ATOMIC(col*sizeof(int));
2666: for ( i = 0; i < row; i++ )
2667: rp[i] = i;
2668: for ( k = 0, d = 0; k < col; k++ ) {
2669: for ( i = d; i < row && !a[i][k]; i++ );
2670: if ( i == row ) {
2671: cp[k] = 0;
2672: continue;
2673: } else
2674: cp[k] = 1;
2675: if ( i != d ) {
2676: j = rp[i]; rp[i] = rp[d]; rp[d] = j;
2677: t = a[i]; a[i] = a[d]; a[d] = t;
2678: }
2679: pivot = a[d];
2680: pivot[k] = inv = invm(pivot[k],md);
2681: for ( i = d+1; i < row; i++ ) {
2682: t = a[i];
2683: if ( m = t[k] ) {
2684: DMAR(inv,m,0,md,t[k])
2685: for ( j = k+1, m = md - t[k]; j < col; j++ )
2686: if ( pivot[j] ) {
1.8 noro 2687: unsigned int tj;
2688: DMAR(m,pivot[j],t[j],md,tj)
2689: t[j] = tj;
1.3 noro 2690: }
2691: }
2692: }
2693: d++;
2694: }
2695: return d;
2696: }
2697:
1.53 noro 2698: int lu_mod(unsigned int **a,int n,unsigned int md,int **rinfo)
2699: {
2700: int i,j,k;
2701: int *rp;
2702: unsigned int *t,*pivot;
2703: unsigned int inv,m;
2704:
2705: *rinfo = rp = (int *)MALLOC_ATOMIC(n*sizeof(int));
2706: for ( i = 0; i < n; i++ ) rp[i] = i;
2707: for ( k = 0; k < n; k++ ) {
2708: for ( i = k; i < n && !a[i][k]; i++ );
2709: if ( i == n ) return 0;
2710: if ( i != k ) {
2711: j = rp[i]; rp[i] = rp[k]; rp[k] = j;
2712: t = a[i]; a[i] = a[k]; a[k] = t;
2713: }
2714: pivot = a[k];
2715: inv = invm(pivot[k],md);
2716: for ( i = k+1; i < n; i++ ) {
2717: t = a[i];
2718: if ( m = t[k] ) {
2719: DMAR(inv,m,0,md,t[k])
2720: for ( j = k+1, m = md - t[k]; j < n; j++ )
2721: if ( pivot[j] ) {
2722: unsigned int tj;
2723: DMAR(m,pivot[j],t[j],md,tj)
2724: t[j] = tj;
2725: }
2726: }
2727: }
2728: }
2729: return 1;
2730: }
2731:
1.3 noro 2732: /*
2733: Input
2734: a : n x n matrix; a result of LU-decomposition
2735: md : modulus
2736: b : n x l matrix
2737: Output
2738: b = a^(-1)b
2739: */
2740:
1.44 noro 2741: void solve_by_lu_mod(int **a,int n,int md,int **b,int l,int normalize)
1.3 noro 2742: {
2743: unsigned int *y,*c;
2744: int i,j,k;
2745: unsigned int t,m,m2;
2746:
2747: y = (int *)MALLOC_ATOMIC(n*sizeof(int));
2748: c = (int *)MALLOC_ATOMIC(n*sizeof(int));
2749: m2 = md>>1;
2750: for ( k = 0; k < l; k++ ) {
2751: /* copy b[.][k] to c */
2752: for ( i = 0; i < n; i++ )
2753: c[i] = (unsigned int)b[i][k];
2754: /* solve Ly=c */
2755: for ( i = 0; i < n; i++ ) {
2756: for ( t = c[i], j = 0; j < i; j++ )
2757: if ( a[i][j] ) {
2758: m = md - a[i][j];
2759: DMAR(m,y[j],t,md,t)
2760: }
2761: y[i] = t;
2762: }
2763: /* solve Uc=y */
2764: for ( i = n-1; i >= 0; i-- ) {
2765: for ( t = y[i], j =i+1; j < n; j++ )
2766: if ( a[i][j] ) {
2767: m = md - a[i][j];
2768: DMAR(m,c[j],t,md,t)
2769: }
2770: /* a[i][i] = 1/U[i][i] */
2771: DMAR(t,a[i][i],0,md,c[i])
2772: }
2773: /* copy c to b[.][k] with normalization */
1.44 noro 2774: if ( normalize )
2775: for ( i = 0; i < n; i++ )
2776: b[i][k] = (int)(c[i]>m2 ? c[i]-md : c[i]);
2777: else
2778: for ( i = 0; i < n; i++ )
2779: b[i][k] = c[i];
1.3 noro 2780: }
2781: }
2782:
1.24 noro 2783: void Pleqm1(NODE arg,VECT *rp)
1.1 noro 2784: {
2785: MAT m;
2786: VECT vect;
2787: pointer **mat;
2788: Q *v;
2789: Q q;
2790: int **wmat;
2791: int md,i,j,row,col,t,n,status;
2792:
2793: asir_assert(ARG0(arg),O_MAT,"leqm1");
2794: asir_assert(ARG1(arg),O_N,"leqm1");
2795: m = (MAT)ARG0(arg); md = QTOS((Q)ARG1(arg));
2796: row = m->row; col = m->col; mat = m->body;
2797: wmat = (int **)almat(row,col);
2798: for ( i = 0; i < row; i++ )
2799: for ( j = 0; j < col; j++ )
2800: if ( q = (Q)mat[i][j] ) {
2801: t = rem(NM(q),md);
2802: if ( SGN(q) < 0 )
2803: t = (md - t) % md;
2804: wmat[i][j] = t;
2805: } else
2806: wmat[i][j] = 0;
2807: status = gauss_elim_mod1(wmat,row,col,md);
2808: if ( status < 0 )
2809: *rp = 0;
2810: else if ( status > 0 )
2811: *rp = (VECT)ONE;
2812: else {
2813: n = col - 1;
2814: MKVECT(vect,n);
2815: for ( i = 0, v = (Q *)vect->body; i < n; i++ ) {
2816: t = (md-wmat[i][n])%md; STOQ(t,v[i]);
2817: }
2818: *rp = vect;
2819: }
2820: }
2821:
1.24 noro 2822: int gauss_elim_mod1(int **mat,int row,int col,int md)
1.1 noro 2823: {
2824: int i,j,k,inv,a,n;
2825: int *t,*pivot;
2826:
2827: n = col - 1;
2828: for ( j = 0; j < n; j++ ) {
2829: for ( i = j; i < row && !mat[i][j]; i++ );
2830: if ( i == row )
2831: return 1;
2832: if ( i != j ) {
2833: t = mat[i]; mat[i] = mat[j]; mat[j] = t;
2834: }
2835: pivot = mat[j];
2836: inv = invm(pivot[j],md);
2837: for ( k = j; k <= n; k++ )
2838: pivot[k] = dmar(pivot[k],inv,0,md);
2839: for ( i = j+1; i < row; i++ ) {
2840: t = mat[i];
2841: if ( i != j && (a = t[j]) )
2842: for ( k = j, a = md - a; k <= n; k++ )
2843: t[k] = dmar(pivot[k],a,t[k],md);
2844: }
2845: }
2846: for ( i = n; i < row && !mat[i][n]; i++ );
2847: if ( i == row ) {
2848: for ( j = n-1; j >= 0; j-- ) {
2849: for ( i = j-1, a = (md-mat[j][n])%md; i >= 0; i-- ) {
2850: mat[i][n] = dmar(mat[i][j],a,mat[i][n],md);
2851: mat[i][j] = 0;
2852: }
2853: }
2854: return 0;
2855: } else
2856: return -1;
2857: }
2858:
1.24 noro 2859: void Pgeninvm(NODE arg,LIST *rp)
1.1 noro 2860: {
2861: MAT m;
2862: pointer **mat;
2863: Q **tmat;
2864: Q q;
2865: unsigned int **wmat;
2866: int md,i,j,row,col,t,status;
2867: MAT mat1,mat2;
2868: NODE node1,node2;
2869:
2870: asir_assert(ARG0(arg),O_MAT,"leqm1");
2871: asir_assert(ARG1(arg),O_N,"leqm1");
2872: m = (MAT)ARG0(arg); md = QTOS((Q)ARG1(arg));
2873: row = m->row; col = m->col; mat = m->body;
2874: wmat = (unsigned int **)almat(row,col+row);
2875: for ( i = 0; i < row; i++ ) {
2876: bzero((char *)wmat[i],(col+row)*sizeof(int));
2877: for ( j = 0; j < col; j++ )
2878: if ( q = (Q)mat[i][j] ) {
2879: t = rem(NM(q),md);
2880: if ( SGN(q) < 0 )
2881: t = (md - t) % md;
2882: wmat[i][j] = t;
2883: }
2884: wmat[i][col+i] = 1;
2885: }
2886: status = gauss_elim_geninv_mod(wmat,row,col,md);
2887: if ( status > 0 )
2888: *rp = 0;
2889: else {
2890: MKMAT(mat1,col,row); MKMAT(mat2,row-col,row);
2891: for ( i = 0, tmat = (Q **)mat1->body; i < col; i++ )
2892: for ( j = 0; j < row; j++ )
1.24 noro 2893: UTOQ(wmat[i][j+col],tmat[i][j]);
1.1 noro 2894: for ( tmat = (Q **)mat2->body; i < row; i++ )
2895: for ( j = 0; j < row; j++ )
1.24 noro 2896: UTOQ(wmat[i][j+col],tmat[i-col][j]);
1.1 noro 2897: MKNODE(node2,mat2,0); MKNODE(node1,mat1,node2); MKLIST(*rp,node1);
2898: }
2899: }
2900:
1.24 noro 2901: int gauss_elim_geninv_mod(unsigned int **mat,int row,int col,int md)
1.1 noro 2902: {
2903: int i,j,k,inv,a,n,m;
2904: unsigned int *t,*pivot;
2905:
2906: n = col; m = row+col;
2907: for ( j = 0; j < n; j++ ) {
2908: for ( i = j; i < row && !mat[i][j]; i++ );
2909: if ( i == row )
2910: return 1;
2911: if ( i != j ) {
2912: t = mat[i]; mat[i] = mat[j]; mat[j] = t;
2913: }
2914: pivot = mat[j];
2915: inv = invm(pivot[j],md);
2916: for ( k = j; k < m; k++ )
2917: pivot[k] = dmar(pivot[k],inv,0,md);
2918: for ( i = j+1; i < row; i++ ) {
2919: t = mat[i];
2920: if ( a = t[j] )
2921: for ( k = j, a = md - a; k < m; k++ )
2922: t[k] = dmar(pivot[k],a,t[k],md);
2923: }
2924: }
2925: for ( j = n-1; j >= 0; j-- ) {
2926: pivot = mat[j];
2927: for ( i = j-1; i >= 0; i-- ) {
2928: t = mat[i];
2929: if ( a = t[j] )
2930: for ( k = j, a = md - a; k < m; k++ )
2931: t[k] = dmar(pivot[k],a,t[k],md);
2932: }
2933: }
2934: return 0;
2935: }
2936:
1.24 noro 2937: void Psolve_by_lu_gfmmat(NODE arg,VECT *rp)
1.1 noro 2938: {
2939: GFMMAT lu;
2940: Q *perm,*rhs,*v;
2941: int n,i;
2942: unsigned int md;
2943: unsigned int *b,*sol;
2944: VECT r;
2945:
2946: lu = (GFMMAT)ARG0(arg);
2947: perm = (Q *)BDY((VECT)ARG1(arg));
2948: rhs = (Q *)BDY((VECT)ARG2(arg));
2949: md = (unsigned int)QTOS((Q)ARG3(arg));
2950: n = lu->col;
2951: b = (unsigned int *)MALLOC_ATOMIC(n*sizeof(int));
2952: sol = (unsigned int *)MALLOC_ATOMIC(n*sizeof(int));
2953: for ( i = 0; i < n; i++ )
2954: b[i] = QTOS(rhs[QTOS(perm[i])]);
2955: solve_by_lu_gfmmat(lu,md,b,sol);
2956: MKVECT(r,n);
2957: for ( i = 0, v = (Q *)r->body; i < n; i++ )
1.24 noro 2958: UTOQ(sol[i],v[i]);
1.1 noro 2959: *rp = r;
2960: }
2961:
1.24 noro 2962: void solve_by_lu_gfmmat(GFMMAT lu,unsigned int md,
2963: unsigned int *b,unsigned int *x)
1.1 noro 2964: {
2965: int n;
2966: unsigned int **a;
2967: unsigned int *y;
2968: int i,j;
2969: unsigned int t,m;
2970:
2971: n = lu->col;
2972: a = lu->body;
2973: y = (unsigned int *)MALLOC_ATOMIC(n*sizeof(int));
2974: /* solve Ly=b */
2975: for ( i = 0; i < n; i++ ) {
2976: for ( t = b[i], j = 0; j < i; j++ )
2977: if ( a[i][j] ) {
2978: m = md - a[i][j];
2979: DMAR(m,y[j],t,md,t)
2980: }
2981: y[i] = t;
2982: }
2983: /* solve Ux=y */
2984: for ( i = n-1; i >= 0; i-- ) {
2985: for ( t = y[i], j =i+1; j < n; j++ )
2986: if ( a[i][j] ) {
2987: m = md - a[i][j];
2988: DMAR(m,x[j],t,md,t)
2989: }
2990: /* a[i][i] = 1/U[i][i] */
2991: DMAR(t,a[i][i],0,md,x[i])
2992: }
2993: }
2994:
1.53 noro 2995: void Plu_mat(NODE arg,LIST *rp)
2996: {
2997: MAT m,lu;
2998: Q dn;
2999: Q *v;
3000: int n,i;
3001: int *iperm;
3002: VECT perm;
3003: NODE n0;
3004:
3005: asir_assert(ARG0(arg),O_MAT,"lu_mat");
3006: m = (MAT)ARG0(arg);
3007: n = m->row;
3008: MKMAT(lu,n,n);
3009: lu_dec_cr(m,lu,&dn,&iperm);
3010: MKVECT(perm,n);
3011: for ( i = 0, v = (Q *)perm->body; i < n; i++ )
3012: STOQ(iperm[i],v[i]);
3013: n0 = mknode(3,lu,dn,perm);
3014: MKLIST(*rp,n0);
3015: }
3016:
1.24 noro 3017: void Plu_gfmmat(NODE arg,LIST *rp)
1.1 noro 3018: {
3019: MAT m;
3020: GFMMAT mm;
3021: unsigned int md;
3022: int i,row,col,status;
3023: int *iperm;
3024: Q *v;
3025: VECT perm;
3026: NODE n0;
3027:
1.53 noro 3028: asir_assert(ARG0(arg),O_MAT,"lu_gfmmat");
3029: asir_assert(ARG1(arg),O_N,"lu_gfmmat");
1.1 noro 3030: m = (MAT)ARG0(arg); md = (unsigned int)QTOS((Q)ARG1(arg));
3031: mat_to_gfmmat(m,md,&mm);
3032: row = m->row;
3033: col = m->col;
3034: iperm = (int *)MALLOC_ATOMIC(row*sizeof(int));
3035: status = lu_gfmmat(mm,md,iperm);
3036: if ( !status )
3037: n0 = 0;
3038: else {
3039: MKVECT(perm,row);
3040: for ( i = 0, v = (Q *)perm->body; i < row; i++ )
3041: STOQ(iperm[i],v[i]);
3042: n0 = mknode(2,mm,perm);
3043: }
3044: MKLIST(*rp,n0);
3045: }
3046:
1.24 noro 3047: void Pmat_to_gfmmat(NODE arg,GFMMAT *rp)
1.1 noro 3048: {
3049: MAT m;
3050: unsigned int md;
3051:
3052: asir_assert(ARG0(arg),O_MAT,"mat_to_gfmmat");
3053: asir_assert(ARG1(arg),O_N,"mat_to_gfmmat");
3054: m = (MAT)ARG0(arg); md = (unsigned int)QTOS((Q)ARG1(arg));
3055: mat_to_gfmmat(m,md,rp);
3056: }
3057:
1.24 noro 3058: void mat_to_gfmmat(MAT m,unsigned int md,GFMMAT *rp)
1.1 noro 3059: {
3060: unsigned int **wmat;
3061: unsigned int t;
3062: Q **mat;
3063: Q q;
3064: int i,j,row,col;
3065:
3066: row = m->row; col = m->col; mat = (Q **)m->body;
3067: wmat = (unsigned int **)almat(row,col);
3068: for ( i = 0; i < row; i++ ) {
3069: bzero((char *)wmat[i],col*sizeof(unsigned int));
3070: for ( j = 0; j < col; j++ )
3071: if ( q = mat[i][j] ) {
3072: t = (unsigned int)rem(NM(q),md);
3073: if ( SGN(q) < 0 )
3074: t = (md - t) % md;
3075: wmat[i][j] = t;
3076: }
3077: }
3078: TOGFMMAT(row,col,wmat,*rp);
3079: }
3080:
1.27 noro 3081: void Pgeninvm_swap(arg,rp)
3082: NODE arg;
3083: LIST *rp;
1.1 noro 3084: {
3085: MAT m;
3086: pointer **mat;
3087: Q **tmat;
3088: Q *tvect;
3089: Q q;
3090: unsigned int **wmat,**invmat;
3091: int *index;
3092: unsigned int t,md;
3093: int i,j,row,col,status;
3094: MAT mat1;
3095: VECT vect1;
3096: NODE node1,node2;
3097:
3098: asir_assert(ARG0(arg),O_MAT,"geninvm_swap");
3099: asir_assert(ARG1(arg),O_N,"geninvm_swap");
3100: m = (MAT)ARG0(arg); md = QTOS((Q)ARG1(arg));
3101: row = m->row; col = m->col; mat = m->body;
3102: wmat = (unsigned int **)almat(row,col+row);
3103: for ( i = 0; i < row; i++ ) {
3104: bzero((char *)wmat[i],(col+row)*sizeof(int));
3105: for ( j = 0; j < col; j++ )
3106: if ( q = (Q)mat[i][j] ) {
3107: t = (unsigned int)rem(NM(q),md);
3108: if ( SGN(q) < 0 )
3109: t = (md - t) % md;
3110: wmat[i][j] = t;
3111: }
3112: wmat[i][col+i] = 1;
3113: }
3114: status = gauss_elim_geninv_mod_swap(wmat,row,col,md,&invmat,&index);
3115: if ( status > 0 )
3116: *rp = 0;
3117: else {
3118: MKMAT(mat1,col,col);
3119: for ( i = 0, tmat = (Q **)mat1->body; i < col; i++ )
3120: for ( j = 0; j < col; j++ )
3121: UTOQ(invmat[i][j],tmat[i][j]);
3122: MKVECT(vect1,row);
3123: for ( i = 0, tvect = (Q *)vect1->body; i < row; i++ )
3124: STOQ(index[i],tvect[i]);
3125: MKNODE(node2,vect1,0); MKNODE(node1,mat1,node2); MKLIST(*rp,node1);
3126: }
3127: }
3128:
1.27 noro 3129: gauss_elim_geninv_mod_swap(mat,row,col,md,invmatp,indexp)
3130: unsigned int **mat;
3131: int row,col;
3132: unsigned int md;
3133: unsigned int ***invmatp;
3134: int **indexp;
1.1 noro 3135: {
3136: int i,j,k,inv,a,n,m;
3137: unsigned int *t,*pivot,*s;
3138: int *index;
3139: unsigned int **invmat;
3140:
3141: n = col; m = row+col;
3142: *indexp = index = (int *)MALLOC_ATOMIC(row*sizeof(int));
3143: for ( i = 0; i < row; i++ )
3144: index[i] = i;
3145: for ( j = 0; j < n; j++ ) {
3146: for ( i = j; i < row && !mat[i][j]; i++ );
3147: if ( i == row ) {
3148: *indexp = 0; *invmatp = 0; return 1;
3149: }
3150: if ( i != j ) {
3151: t = mat[i]; mat[i] = mat[j]; mat[j] = t;
3152: k = index[i]; index[i] = index[j]; index[j] = k;
3153: }
3154: pivot = mat[j];
3155: inv = (unsigned int)invm(pivot[j],md);
3156: for ( k = j; k < m; k++ )
3157: if ( pivot[k] )
3158: pivot[k] = (unsigned int)dmar(pivot[k],inv,0,md);
3159: for ( i = j+1; i < row; i++ ) {
3160: t = mat[i];
3161: if ( a = t[j] )
3162: for ( k = j, a = md - a; k < m; k++ )
3163: if ( pivot[k] )
3164: t[k] = dmar(pivot[k],a,t[k],md);
3165: }
3166: }
3167: for ( j = n-1; j >= 0; j-- ) {
3168: pivot = mat[j];
3169: for ( i = j-1; i >= 0; i-- ) {
3170: t = mat[i];
3171: if ( a = t[j] )
3172: for ( k = j, a = md - a; k < m; k++ )
3173: if ( pivot[k] )
3174: t[k] = dmar(pivot[k],a,t[k],md);
3175: }
3176: }
3177: *invmatp = invmat = (unsigned int **)almat(col,col);
1.27 noro 3178: for ( i = 0; i < col; i++ )
3179: for ( j = 0, s = invmat[i], t = mat[i]; j < col; j++ )
3180: s[j] = t[col+index[j]];
3181: return 0;
3182: }
3183:
3184: void Pgeninv_sf_swap(NODE arg,LIST *rp)
3185: {
3186: MAT m;
3187: GFS **mat,**tmat;
3188: Q *tvect;
3189: GFS q;
3190: int **wmat,**invmat;
3191: int *index;
3192: unsigned int t;
3193: int i,j,row,col,status;
3194: MAT mat1;
3195: VECT vect1;
3196: NODE node1,node2;
3197:
3198: asir_assert(ARG0(arg),O_MAT,"geninv_sf_swap");
3199: m = (MAT)ARG0(arg);
3200: row = m->row; col = m->col; mat = (GFS **)m->body;
3201: wmat = (int **)almat(row,col+row);
3202: for ( i = 0; i < row; i++ ) {
3203: bzero((char *)wmat[i],(col+row)*sizeof(int));
3204: for ( j = 0; j < col; j++ )
3205: if ( q = (GFS)mat[i][j] )
3206: wmat[i][j] = FTOIF(CONT(q));
3207: wmat[i][col+i] = _onesf();
3208: }
3209: status = gauss_elim_geninv_sf_swap(wmat,row,col,&invmat,&index);
3210: if ( status > 0 )
3211: *rp = 0;
3212: else {
3213: MKMAT(mat1,col,col);
3214: for ( i = 0, tmat = (GFS **)mat1->body; i < col; i++ )
3215: for ( j = 0; j < col; j++ )
3216: if ( t = invmat[i][j] ) {
3217: MKGFS(IFTOF(t),tmat[i][j]);
3218: }
3219: MKVECT(vect1,row);
3220: for ( i = 0, tvect = (Q *)vect1->body; i < row; i++ )
3221: STOQ(index[i],tvect[i]);
3222: MKNODE(node2,vect1,0); MKNODE(node1,mat1,node2); MKLIST(*rp,node1);
3223: }
3224: }
3225:
3226: int gauss_elim_geninv_sf_swap(int **mat,int row,int col,
3227: int ***invmatp,int **indexp)
3228: {
3229: int i,j,k,inv,a,n,m,u;
3230: int *t,*pivot,*s;
3231: int *index;
3232: int **invmat;
3233:
3234: n = col; m = row+col;
3235: *indexp = index = (int *)MALLOC_ATOMIC(row*sizeof(int));
3236: for ( i = 0; i < row; i++ )
3237: index[i] = i;
3238: for ( j = 0; j < n; j++ ) {
3239: for ( i = j; i < row && !mat[i][j]; i++ );
3240: if ( i == row ) {
3241: *indexp = 0; *invmatp = 0; return 1;
3242: }
3243: if ( i != j ) {
3244: t = mat[i]; mat[i] = mat[j]; mat[j] = t;
3245: k = index[i]; index[i] = index[j]; index[j] = k;
3246: }
3247: pivot = mat[j];
3248: inv = _invsf(pivot[j]);
3249: for ( k = j; k < m; k++ )
3250: if ( pivot[k] )
3251: pivot[k] = _mulsf(pivot[k],inv);
3252: for ( i = j+1; i < row; i++ ) {
3253: t = mat[i];
3254: if ( a = t[j] )
3255: for ( k = j, a = _chsgnsf(a); k < m; k++ )
3256: if ( pivot[k] ) {
3257: u = _mulsf(pivot[k],a);
3258: t[k] = _addsf(u,t[k]);
3259: }
3260: }
3261: }
3262: for ( j = n-1; j >= 0; j-- ) {
3263: pivot = mat[j];
3264: for ( i = j-1; i >= 0; i-- ) {
3265: t = mat[i];
3266: if ( a = t[j] )
3267: for ( k = j, a = _chsgnsf(a); k < m; k++ )
3268: if ( pivot[k] ) {
3269: u = _mulsf(pivot[k],a);
3270: t[k] = _addsf(u,t[k]);
3271: }
3272: }
3273: }
3274: *invmatp = invmat = (int **)almat(col,col);
1.1 noro 3275: for ( i = 0; i < col; i++ )
3276: for ( j = 0, s = invmat[i], t = mat[i]; j < col; j++ )
3277: s[j] = t[col+index[j]];
3278: return 0;
3279: }
3280:
3281: void _addn(N,N,N);
3282: int _subn(N,N,N);
3283: void _muln(N,N,N);
3284:
1.24 noro 3285: void inner_product_int(Q *a,Q *b,int n,Q *r)
1.1 noro 3286: {
3287: int la,lb,i;
3288: int sgn,sgn1;
3289: N wm,wma,sum,t;
3290:
3291: for ( la = lb = 0, i = 0; i < n; i++ ) {
3292: if ( a[i] )
3293: if ( DN(a[i]) )
3294: error("inner_product_int : invalid argument");
3295: else
3296: la = MAX(PL(NM(a[i])),la);
3297: if ( b[i] )
3298: if ( DN(b[i]) )
3299: error("inner_product_int : invalid argument");
3300: else
3301: lb = MAX(PL(NM(b[i])),lb);
3302: }
3303: sgn = 0;
3304: sum= NALLOC(la+lb+2);
3305: bzero((char *)sum,(la+lb+3)*sizeof(unsigned int));
3306: wm = NALLOC(la+lb+2);
3307: wma = NALLOC(la+lb+2);
3308: for ( i = 0; i < n; i++ ) {
3309: if ( !a[i] || !b[i] )
3310: continue;
3311: _muln(NM(a[i]),NM(b[i]),wm);
3312: sgn1 = SGN(a[i])*SGN(b[i]);
3313: if ( !sgn ) {
3314: sgn = sgn1;
3315: t = wm; wm = sum; sum = t;
3316: } else if ( sgn == sgn1 ) {
3317: _addn(sum,wm,wma);
3318: if ( !PL(wma) )
3319: sgn = 0;
3320: t = wma; wma = sum; sum = t;
3321: } else {
3322: /* sgn*sum+sgn1*wm = sgn*(sum-wm) */
3323: sgn *= _subn(sum,wm,wma);
3324: t = wma; wma = sum; sum = t;
3325: }
3326: }
1.61 noro 3327: GCFREE(wm);
3328: GCFREE(wma);
1.1 noro 3329: if ( !sgn ) {
1.61 noro 3330: GCFREE(sum);
1.1 noro 3331: *r = 0;
3332: } else
3333: NTOQ(sum,sgn,*r);
3334: }
3335:
1.3 noro 3336: /* (k,l) element of a*b where a: .x n matrix, b: n x . integer matrix */
3337:
1.24 noro 3338: void inner_product_mat_int_mod(Q **a,int **b,int n,int k,int l,Q *r)
1.3 noro 3339: {
3340: int la,lb,i;
3341: int sgn,sgn1;
3342: N wm,wma,sum,t;
3343: Q aki;
3344: int bil,bilsgn;
3345: struct oN tn;
3346:
3347: for ( la = 0, i = 0; i < n; i++ ) {
3348: if ( aki = a[k][i] )
3349: if ( DN(aki) )
3350: error("inner_product_int : invalid argument");
3351: else
3352: la = MAX(PL(NM(aki)),la);
3353: }
3354: lb = 1;
3355: sgn = 0;
3356: sum= NALLOC(la+lb+2);
3357: bzero((char *)sum,(la+lb+3)*sizeof(unsigned int));
3358: wm = NALLOC(la+lb+2);
3359: wma = NALLOC(la+lb+2);
3360: for ( i = 0; i < n; i++ ) {
3361: if ( !(aki = a[k][i]) || !(bil = b[i][l]) )
3362: continue;
3363: tn.p = 1;
3364: if ( bil > 0 ) {
3365: tn.b[0] = bil; bilsgn = 1;
3366: } else {
3367: tn.b[0] = -bil; bilsgn = -1;
3368: }
3369: _muln(NM(aki),&tn,wm);
3370: sgn1 = SGN(aki)*bilsgn;
3371: if ( !sgn ) {
3372: sgn = sgn1;
3373: t = wm; wm = sum; sum = t;
3374: } else if ( sgn == sgn1 ) {
3375: _addn(sum,wm,wma);
3376: if ( !PL(wma) )
3377: sgn = 0;
3378: t = wma; wma = sum; sum = t;
3379: } else {
3380: /* sgn*sum+sgn1*wm = sgn*(sum-wm) */
3381: sgn *= _subn(sum,wm,wma);
3382: t = wma; wma = sum; sum = t;
3383: }
3384: }
1.61 noro 3385: GCFREE(wm);
3386: GCFREE(wma);
1.3 noro 3387: if ( !sgn ) {
1.61 noro 3388: GCFREE(sum);
1.3 noro 3389: *r = 0;
3390: } else
3391: NTOQ(sum,sgn,*r);
3392: }
3393:
1.24 noro 3394: void Pmul_mat_vect_int(NODE arg,VECT *rp)
1.1 noro 3395: {
3396: MAT mat;
3397: VECT vect,r;
3398: int row,col,i;
3399:
3400: mat = (MAT)ARG0(arg);
3401: vect = (VECT)ARG1(arg);
3402: row = mat->row;
3403: col = mat->col;
3404: MKVECT(r,row);
1.24 noro 3405: for ( i = 0; i < row; i++ ) {
3406: inner_product_int((Q *)mat->body[i],(Q *)vect->body,col,(Q *)&r->body[i]);
3407: }
1.1 noro 3408: *rp = r;
3409: }
3410:
1.24 noro 3411: void Pnbpoly_up2(NODE arg,GF2N *rp)
1.1 noro 3412: {
3413: int m,type,ret;
3414: UP2 r;
3415:
3416: m = QTOS((Q)ARG0(arg));
3417: type = QTOS((Q)ARG1(arg));
3418: ret = generate_ONB_polynomial(&r,m,type);
3419: if ( ret == 0 )
3420: MKGF2N(r,*rp);
3421: else
3422: *rp = 0;
3423: }
3424:
1.24 noro 3425: void Px962_irredpoly_up2(NODE arg,GF2N *rp)
1.1 noro 3426: {
1.24 noro 3427: int m,ret,w;
1.1 noro 3428: GF2N prev;
3429: UP2 r;
3430:
3431: m = QTOS((Q)ARG0(arg));
3432: prev = (GF2N)ARG1(arg);
3433: if ( !prev ) {
3434: w = (m>>5)+1; NEWUP2(r,w); r->w = 0;
3435: bzero((char *)r->b,w*sizeof(unsigned int));
3436: } else {
3437: r = prev->body;
3438: if ( degup2(r) != m ) {
3439: w = (m>>5)+1; NEWUP2(r,w); r->w = 0;
3440: bzero((char *)r->b,w*sizeof(unsigned int));
3441: }
3442: }
1.24 noro 3443: ret = _generate_irreducible_polynomial(r,m);
1.1 noro 3444: if ( ret == 0 )
3445: MKGF2N(r,*rp);
3446: else
3447: *rp = 0;
3448: }
3449:
1.24 noro 3450: void Pirredpoly_up2(NODE arg,GF2N *rp)
1.1 noro 3451: {
1.24 noro 3452: int m,ret,w;
1.1 noro 3453: GF2N prev;
3454: UP2 r;
3455:
3456: m = QTOS((Q)ARG0(arg));
3457: prev = (GF2N)ARG1(arg);
3458: if ( !prev ) {
3459: w = (m>>5)+1; NEWUP2(r,w); r->w = 0;
3460: bzero((char *)r->b,w*sizeof(unsigned int));
3461: } else {
3462: r = prev->body;
3463: if ( degup2(r) != m ) {
3464: w = (m>>5)+1; NEWUP2(r,w); r->w = 0;
3465: bzero((char *)r->b,w*sizeof(unsigned int));
3466: }
3467: }
1.24 noro 3468: ret = _generate_good_irreducible_polynomial(r,m);
1.1 noro 3469: if ( ret == 0 )
3470: MKGF2N(r,*rp);
3471: else
3472: *rp = 0;
3473: }
3474:
1.26 noro 3475: void Pmat_swap_row_destructive(NODE arg, MAT *m)
3476: {
3477: int i1,i2;
3478: pointer *t;
3479: MAT mat;
3480:
3481: asir_assert(ARG0(arg),O_MAT,"mat_swap_row_destructive");
3482: asir_assert(ARG1(arg),O_N,"mat_swap_row_destructive");
3483: asir_assert(ARG2(arg),O_N,"mat_swap_row_destructive");
3484: mat = (MAT)ARG0(arg);
3485: i1 = QTOS((Q)ARG1(arg));
3486: i2 = QTOS((Q)ARG2(arg));
3487: if ( i1 < 0 || i2 < 0 || i1 >= mat->row || i2 >= mat->row )
3488: error("mat_swap_row_destructive : Out of range");
3489: t = mat->body[i1];
3490: mat->body[i1] = mat->body[i2];
3491: mat->body[i2] = t;
3492: *m = mat;
3493: }
3494:
3495: void Pmat_swap_col_destructive(NODE arg, MAT *m)
3496: {
3497: int j1,j2,i,n;
3498: pointer *mi;
3499: pointer t;
3500: MAT mat;
3501:
3502: asir_assert(ARG0(arg),O_MAT,"mat_swap_col_destructive");
3503: asir_assert(ARG1(arg),O_N,"mat_swap_col_destructive");
3504: asir_assert(ARG2(arg),O_N,"mat_swap_col_destructive");
3505: mat = (MAT)ARG0(arg);
3506: j1 = QTOS((Q)ARG1(arg));
3507: j2 = QTOS((Q)ARG2(arg));
3508: if ( j1 < 0 || j2 < 0 || j1 >= mat->col || j2 >= mat->col )
3509: error("mat_swap_col_destructive : Out of range");
3510: n = mat->row;
3511: for ( i = 0; i < n; i++ ) {
3512: mi = mat->body[i];
3513: t = mi[j1]; mi[j1] = mi[j2]; mi[j2] = t;
3514: }
3515: *m = mat;
3516: }
1.1 noro 3517: /*
3518: * f = type 'type' normal polynomial of degree m if exists
3519: * IEEE P1363 A.7.2
3520: *
3521: * return value : 0 --- exists
3522: * 1 --- does not exist
3523: * -1 --- failure (memory allocation error)
3524: */
3525:
3526: int generate_ONB_polynomial(UP2 *rp,int m,int type)
3527: {
3528: int i,r;
3529: int w;
3530: UP2 f,f0,f1,f2,t;
3531:
3532: w = (m>>5)+1;
3533: switch ( type ) {
3534: case 1:
3535: if ( !TypeT_NB_check(m,1) ) return 1;
3536: NEWUP2(f,w); *rp = f; f->w = w;
3537: /* set all the bits */
3538: for ( i = 0; i < w; i++ )
3539: f->b[i] = 0xffffffff;
3540: /* mask the top word if necessary */
3541: if ( r = (m+1)&31 )
3542: f->b[w-1] &= (1<<r)-1;
3543: return 0;
3544: break;
3545: case 2:
3546: if ( !TypeT_NB_check(m,2) ) return 1;
3547: NEWUP2(f,w); *rp = f;
3548: W_NEWUP2(f0,w);
3549: W_NEWUP2(f1,w);
3550: W_NEWUP2(f2,w);
3551:
3552: /* recursion for genrating Type II normal polynomial */
3553:
3554: /* f0 = 1, f1 = t+1 */
3555: f0->w = 1; f0->b[0] = 1;
3556: f1->w = 1; f1->b[0] = 3;
3557: for ( i = 2; i <= m; i++ ) {
3558: /* f2 = t*f1+f0 */
3559: _bshiftup2(f1,-1,f2);
3560: _addup2_destructive(f2,f0);
3561: /* cyclic change of the variables */
3562: t = f0; f0 = f1; f1 = f2; f2 = t;
3563: }
3564: _copyup2(f1,f);
3565: return 0;
3566: break;
3567: default:
3568: return -1;
3569: break;
3570: }
3571: }
3572:
3573: /*
3574: * f = an irreducible trinomial or pentanomial of degree d 'after' f
3575: * return value : 0 --- exists
3576: * 1 --- does not exist (exhaustion)
3577: */
3578:
3579: int _generate_irreducible_polynomial(UP2 f,int d)
3580: {
3581: int ret,i,j,k,nz,i0,j0,k0;
3582: int w;
3583: unsigned int *fd;
3584:
3585: /*
3586: * if f = x^d+x^i+1 then i0 <- i, j0 <- 0, k0 <-0.
3587: * if f = x^d+x^k+x^j+x^i+1 (k>j>i) then i0 <- i, j0 <- j, k0 <-k.
3588: * otherwise i0,j0,k0 is set to 0.
3589: */
3590:
3591: fd = f->b;
3592: w = (d>>5)+1;
3593: if ( f->w && (d==degup2(f)) ) {
3594: for ( nz = 0, i = d; i >= 0; i-- )
3595: if ( fd[i>>5]&(1<<(i&31)) ) nz++;
3596: switch ( nz ) {
3597: case 3:
3598: for ( i0 = 1; !(fd[i0>>5]&(1<<(i0&31))) ; i0++ );
3599: /* reset i0-th bit */
3600: fd[i0>>5] &= ~(1<<(i0&31));
3601: j0 = k0 = 0;
3602: break;
3603: case 5:
3604: for ( i0 = 1; !(fd[i0>>5]&(1<<(i0&31))) ; i0++ );
3605: /* reset i0-th bit */
3606: fd[i0>>5] &= ~(1<<(i0&31));
3607: for ( j0 = i0+1; !(fd[j0>>5]&(1<<(j0&31))) ; j0++ );
3608: /* reset j0-th bit */
3609: fd[j0>>5] &= ~(1<<(j0&31));
3610: for ( k0 = j0+1; !(fd[k0>>5]&(1<<(k0&31))) ; k0++ );
3611: /* reset k0-th bit */
3612: fd[k0>>5] &= ~(1<<(k0&31));
3613: break;
3614: default:
3615: f->w = 0; break;
3616: }
3617: } else
3618: f->w = 0;
3619:
3620: if ( !f->w ) {
3621: fd = f->b;
3622: f->w = w; fd[0] |= 1; fd[d>>5] |= (1<<(d&31));
3623: i0 = j0 = k0 = 0;
3624: }
3625: /* if j0 > 0 then f is already a pentanomial */
3626: if ( j0 > 0 ) goto PENTA;
3627:
3628: /* searching for an irreducible trinomial */
3629:
3630: for ( i = 1; 2*i <= d; i++ ) {
3631: /* skip the polynomials 'before' f */
3632: if ( i < i0 ) continue;
3633: if ( i == i0 ) { i0 = 0; continue; }
3634: /* set i-th bit */
3635: fd[i>>5] |= (1<<(i&31));
3636: ret = irredcheck_dddup2(f);
3637: if ( ret == 1 ) return 0;
3638: /* reset i-th bit */
3639: fd[i>>5] &= ~(1<<(i&31));
3640: }
3641:
3642: /* searching for an irreducible pentanomial */
3643: PENTA:
3644: for ( i = 1; i < d; i++ ) {
3645: /* skip the polynomials 'before' f */
3646: if ( i < i0 ) continue;
3647: if ( i == i0 ) i0 = 0;
3648: /* set i-th bit */
3649: fd[i>>5] |= (1<<(i&31));
3650: for ( j = i+1; j < d; j++ ) {
3651: /* skip the polynomials 'before' f */
3652: if ( j < j0 ) continue;
3653: if ( j == j0 ) j0 = 0;
3654: /* set j-th bit */
3655: fd[j>>5] |= (1<<(j&31));
3656: for ( k = j+1; k < d; k++ ) {
3657: /* skip the polynomials 'before' f */
3658: if ( k < k0 ) continue;
3659: else if ( k == k0 ) { k0 = 0; continue; }
3660: /* set k-th bit */
3661: fd[k>>5] |= (1<<(k&31));
3662: ret = irredcheck_dddup2(f);
3663: if ( ret == 1 ) return 0;
3664: /* reset k-th bit */
3665: fd[k>>5] &= ~(1<<(k&31));
3666: }
3667: /* reset j-th bit */
3668: fd[j>>5] &= ~(1<<(j&31));
3669: }
3670: /* reset i-th bit */
3671: fd[i>>5] &= ~(1<<(i&31));
3672: }
3673: /* exhausted */
3674: return 1;
3675: }
3676:
3677: /*
3678: * f = an irreducible trinomial or pentanomial of degree d 'after' f
3679: *
3680: * searching strategy:
3681: * trinomial x^d+x^i+1:
3682: * i is as small as possible.
3683: * trinomial x^d+x^i+x^j+x^k+1:
3684: * i is as small as possible.
3685: * For such i, j is as small as possible.
3686: * For such i and j, 'k' is as small as possible.
3687: *
3688: * return value : 0 --- exists
3689: * 1 --- does not exist (exhaustion)
3690: */
3691:
3692: int _generate_good_irreducible_polynomial(UP2 f,int d)
3693: {
3694: int ret,i,j,k,nz,i0,j0,k0;
3695: int w;
3696: unsigned int *fd;
3697:
3698: /*
3699: * if f = x^d+x^i+1 then i0 <- i, j0 <- 0, k0 <-0.
3700: * if f = x^d+x^k+x^j+x^i+1 (k>j>i) then i0 <- i, j0 <- j, k0 <-k.
3701: * otherwise i0,j0,k0 is set to 0.
3702: */
3703:
3704: fd = f->b;
3705: w = (d>>5)+1;
3706: if ( f->w && (d==degup2(f)) ) {
3707: for ( nz = 0, i = d; i >= 0; i-- )
3708: if ( fd[i>>5]&(1<<(i&31)) ) nz++;
3709: switch ( nz ) {
3710: case 3:
3711: for ( i0 = 1; !(fd[i0>>5]&(1<<(i0&31))) ; i0++ );
3712: /* reset i0-th bit */
3713: fd[i0>>5] &= ~(1<<(i0&31));
3714: j0 = k0 = 0;
3715: break;
3716: case 5:
3717: for ( i0 = 1; !(fd[i0>>5]&(1<<(i0&31))) ; i0++ );
3718: /* reset i0-th bit */
3719: fd[i0>>5] &= ~(1<<(i0&31));
3720: for ( j0 = i0+1; !(fd[j0>>5]&(1<<(j0&31))) ; j0++ );
3721: /* reset j0-th bit */
3722: fd[j0>>5] &= ~(1<<(j0&31));
3723: for ( k0 = j0+1; !(fd[k0>>5]&(1<<(k0&31))) ; k0++ );
3724: /* reset k0-th bit */
3725: fd[k0>>5] &= ~(1<<(k0&31));
3726: break;
3727: default:
3728: f->w = 0; break;
3729: }
3730: } else
3731: f->w = 0;
3732:
3733: if ( !f->w ) {
3734: fd = f->b;
3735: f->w = w; fd[0] |= 1; fd[d>>5] |= (1<<(d&31));
3736: i0 = j0 = k0 = 0;
3737: }
3738: /* if j0 > 0 then f is already a pentanomial */
3739: if ( j0 > 0 ) goto PENTA;
3740:
3741: /* searching for an irreducible trinomial */
3742:
3743: for ( i = 1; 2*i <= d; i++ ) {
3744: /* skip the polynomials 'before' f */
3745: if ( i < i0 ) continue;
3746: if ( i == i0 ) { i0 = 0; continue; }
3747: /* set i-th bit */
3748: fd[i>>5] |= (1<<(i&31));
3749: ret = irredcheck_dddup2(f);
3750: if ( ret == 1 ) return 0;
3751: /* reset i-th bit */
3752: fd[i>>5] &= ~(1<<(i&31));
3753: }
3754:
3755: /* searching for an irreducible pentanomial */
3756: PENTA:
3757: for ( i = 3; i < d; i++ ) {
3758: /* skip the polynomials 'before' f */
3759: if ( i < i0 ) continue;
3760: if ( i == i0 ) i0 = 0;
3761: /* set i-th bit */
3762: fd[i>>5] |= (1<<(i&31));
3763: for ( j = 2; j < i; j++ ) {
3764: /* skip the polynomials 'before' f */
3765: if ( j < j0 ) continue;
3766: if ( j == j0 ) j0 = 0;
3767: /* set j-th bit */
3768: fd[j>>5] |= (1<<(j&31));
3769: for ( k = 1; k < j; k++ ) {
3770: /* skip the polynomials 'before' f */
3771: if ( k < k0 ) continue;
3772: else if ( k == k0 ) { k0 = 0; continue; }
3773: /* set k-th bit */
3774: fd[k>>5] |= (1<<(k&31));
3775: ret = irredcheck_dddup2(f);
3776: if ( ret == 1 ) return 0;
3777: /* reset k-th bit */
3778: fd[k>>5] &= ~(1<<(k&31));
3779: }
3780: /* reset j-th bit */
3781: fd[j>>5] &= ~(1<<(j&31));
3782: }
3783: /* reset i-th bit */
3784: fd[i>>5] &= ~(1<<(i&31));
3785: }
3786: /* exhausted */
3787: return 1;
1.3 noro 3788: }
3789:
1.24 noro 3790: void printqmat(Q **mat,int row,int col)
1.3 noro 3791: {
3792: int i,j;
3793:
3794: for ( i = 0; i < row; i++ ) {
3795: for ( j = 0; j < col; j++ ) {
1.8 noro 3796: printnum((Num)mat[i][j]); printf(" ");
1.3 noro 3797: }
3798: printf("\n");
3799: }
3800: }
3801:
1.24 noro 3802: void printimat(int **mat,int row,int col)
1.3 noro 3803: {
3804: int i,j;
3805:
3806: for ( i = 0; i < row; i++ ) {
3807: for ( j = 0; j < col; j++ ) {
3808: printf("%d ",mat[i][j]);
3809: }
3810: printf("\n");
3811: }
1.36 noro 3812: }
3813:
3814: void Pnd_det(NODE arg,P *rp)
3815: {
1.37 noro 3816: if ( argc(arg) == 1 )
3817: nd_det(0,ARG0(arg),rp);
3818: else
3819: nd_det(QTOS((Q)ARG1(arg)),ARG0(arg),rp);
1.1 noro 3820: }
1.59 ohara 3821:
1.62 ohara 3822: void Pmat_col(NODE arg,VECT *rp)
1.59 ohara 3823: {
3824: int i,j,n;
3825: MAT mat;
3826: VECT vect;
3827:
3828: asir_assert(ARG0(arg),O_MAT,"mat_col");
3829: asir_assert(ARG1(arg),O_N,"mat_col");
3830: mat = (MAT)ARG0(arg);
3831: j = QTOS((Q)ARG1(arg));
3832: if ( j < 0 || j >= mat->col) {
3833: error("mat_col : Out of range");
3834: }
3835: n = mat->row;
3836: MKVECT(vect,n);
3837: for(i=0; i<n; i++) {
3838: BDY(vect)[i] = BDY(mat)[i][j];
3839: }
3840: *rp = vect;
3841: }