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