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