Annotation of OpenXM_contrib2/asir2018/builtin/array.c, Revision 1.2
1.1 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
26: * e-mail at risa-admin@sec.flab.fujitsu.co.jp of the detailed specification
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.2 ! noro 48: * $OpenXM: OpenXM_contrib2/asir2018/builtin/array.c,v 1.1 2018/09/19 05:45:05 noro Exp $
1.1 noro 49: */
50: #include "ca.h"
51: #include "base.h"
52: #include "parse.h"
53: #include "inline.h"
54:
55: #include <sys/types.h>
56: #include <sys/stat.h>
57: #if !defined(_MSC_VER)
58: #include <unistd.h>
59: #endif
60:
61: #define F4_INTRAT_PERIOD 8
62:
63: #if 0
64: #undef DMAR
65: #define DMAR(a1,a2,a3,d,r) (r)=dmar(a1,a2,a3,d);
66: #endif
67:
68: extern int DP_Print; /* XXX */
69:
70:
71: void Pnewvect(), Pnewmat(), Psepvect(), Psize(), Pdet(), Pleqm(), Pleqm1(), Pgeninvm(), Ptriangleq();
72: void Pinvmat();
73: void Pnewbytearray(),Pmemoryplot_to_coord();
74:
75: void Pgeneric_gauss_elim();
76: void Pgeneric_gauss_elim_mod();
77:
78: void Pindep_rows_mod();
79:
80: void Pmat_to_gfmmat(),Plu_gfmmat(),Psolve_by_lu_gfmmat();
81: void Pgeninvm_swap(), Premainder(), Psremainder(), Pvtol(), Pltov();
82: void Pgeninv_sf_swap();
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();
94: void Pexponent_vector();
95: void Pmat_swap_row_destructive();
96: void Pmat_swap_col_destructive();
97: void Pvect();
98: void Pmat();
99: void Pmatc();
100: void Pnd_det();
101: void Plu_mat();
102: void Pmat_col();
103: void Plusolve_prep();
104: void Plusolve_main();
105:
106: struct ftab array_tab[] = {
107: #if 0
108: {"lu_mat",Plu_mat,1},
109: #endif
110: {"solve_by_lu_gfmmat",Psolve_by_lu_gfmmat,4},
111: {"lu_gfmmat",Plu_gfmmat,2},
112: {"mat_to_gfmmat",Pmat_to_gfmmat,2},
113: {"generic_gauss_elim",Pgeneric_gauss_elim,1},
114: {"generic_gauss_elim_mod",Pgeneric_gauss_elim_mod,2},
115: {"indep_rows_mod",Pindep_rows_mod,2},
116: {"newvect",Pnewvect,-2},
117: {"vect",Pvect,-99999999},
118: {"vector",Pnewvect,-2},
119: {"exponent_vector",Pexponent_vector,-99999999},
120: {"newmat",Pnewmat,-3},
121: {"matrix",Pnewmat,-3},
122: {"mat",Pmat,-99999999},
123: {"matr",Pmat,-99999999},
124: {"matc",Pmatc,-99999999},
125: {"newbytearray",Pnewbytearray,-2},
126: {"memoryplot_to_coord",Pmemoryplot_to_coord,1},
127: {"sepmat_destructive",Psepmat_destructive,2},
128: {"sepvect",Psepvect,2},
129: {"qsort",Pqsort,-2},
130: {"vtol",Pvtol,1},
131: {"ltov",Pltov,1},
132: {"size",Psize,1},
133: {"det",Pdet,-2},
134: {"nd_det",Pnd_det,-2},
135: {"invmat",Pinvmat,-2},
136: {"leqm",Pleqm,2},
137: {"leqm1",Pleqm1,2},
138: {"geninvm",Pgeninvm,2},
139: {"geninvm_swap",Pgeninvm_swap,2},
140: {"geninv_sf_swap",Pgeninv_sf_swap,1},
141: {"remainder",Premainder,2},
142: {"sremainder",Psremainder,2},
143: {"mulmat_gf2n",Pmulmat_gf2n,1},
144: {"bconvmat_gf2n",Pbconvmat_gf2n,-4},
145: {"mul_vect_mat_gf2n",Pmul_vect_mat_gf2n,2},
146: {"mul_mat_vect_int",Pmul_mat_vect_int,2},
147: {"nbmul_gf2n",PNBmul_gf2n,3},
148: {"x962_irredpoly_up2",Px962_irredpoly_up2,2},
149: {"irredpoly_up2",Pirredpoly_up2,2},
150: {"nbpoly_up2",Pnbpoly_up2,2},
151: {"mat_swap_row_destructive",Pmat_swap_row_destructive,3},
152: {"mat_swap_col_destructive",Pmat_swap_col_destructive,3},
153: {"mat_col",Pmat_col,2},
154: {"lusolve_prep",Plusolve_prep,1},
155: {"lusolve_main",Plusolve_main,1},
156: {"triangleq",Ptriangleq,1},
157: {0,0,0},
158: };
159:
160: typedef struct _ent { int j; unsigned int e; } ent;
161:
162: ent *get_row(FILE *,int *l);
163: void put_row(FILE *out,int l,ent *a);
164: void lu_elim(int *l,ent **a,int k,int i,int mul,int mod);
165: void lu_append(int *,ent **,int *,int,int,int);
166: void solve_l(int *,ent **,int,int *,int);
167: void solve_u(int *,ent **,int,int *,int);
168:
169:
170: static int *ul,*ll;
171: static ent **u,**l;
172: static int modulus;
173:
174: void Plusolve_prep(NODE arg,Q *rp)
175: {
176: char *fname;
177: FILE *in;
178: int len,i,rank;
179: int *rhs;
180:
181: fname = BDY((STRING)ARG0(arg));
182: in = fopen(fname,"r");
183: modulus = getw(in);
184: len = getw(in);
185: ul = (int *)MALLOC_ATOMIC(len*sizeof(int));
186: u = (ent **)MALLOC(len*sizeof(ent *));
187: ll = (int *)MALLOC_ATOMIC(len*sizeof(int));
188: l = (ent **)MALLOC(len*sizeof(ent *));
189: for ( i = 0; i < len; i++ ) {
190: u[i] = get_row(in,&ul[i]);
191: }
192: for ( i = 0; i < len; i++ ) {
193: l[i] = get_row(in,&ll[i]);
194: }
195: fclose(in);
196: *rp = (Q)ONE;
197: }
198:
199: void Plusolve_main(NODE arg,VECT *rp)
200: {
201: Z *d,*p;
202: VECT v,r;
203: int len,i;
204: int *rhs;
205:
206: v = (VECT)ARG0(arg); len = v->len;
207: d = (Z *)BDY(v);
208: rhs = (int *)MALLOC_ATOMIC(len*sizeof(int));
1.2 ! noro 209: for ( i = 0; i < len; i++ ) rhs[i] = ZTOS(d[i]);
1.1 noro 210: solve_l(ll,l,len,rhs,modulus);
211: solve_u(ul,u,len,rhs,modulus);
212: NEWVECT(r); r->len = len;
213: r->body = (pointer *)MALLOC(len*sizeof(pointer));
214: p = (Z *)r->body;
215: for ( i = 0; i < len; i++ )
1.2 ! noro 216: STOZ(rhs[i],p[i]);
1.1 noro 217: *rp = r;
218: }
219:
220: ent *get_row(FILE *in,int *l)
221: {
222: int len,i;
223: ent *a;
224:
225: *l = len = getw(in);
226: a = (ent *)MALLOC_ATOMIC(len*sizeof(ent));
227: for ( i = 0; i < len; i++ ) {
228: a[i].j = getw(in);
229: a[i].e = getw(in);
230: }
231: return a;
232: }
233:
234: void lu_gauss(int *ul,ent **u,int *ll,ent **l,int n,int mod)
235: {
236: int i,j,k,s,mul;
237: unsigned int inv;
238: int *ll2;
239:
240: ll2 = (int *)MALLOC_ATOMIC(n*sizeof(int));
241: for ( i = 0; i < n; i++ ) ll2[i] = 0;
242: for ( i = 0; i < n; i++ ) {
243: fprintf(stderr,"i=%d\n",i);
244: inv = invm(u[i][0].e,mod);
245: for ( k = i+1; k < n; k++ )
246: if ( u[k][0].j == n-i ) {
247: s = u[k][0].e;
248: DMAR(s,inv,0,mod,mul);
249: lu_elim(ul,u,k,i,mul,mod);
250: lu_append(ll,l,ll2,k,i,mul);
251: }
252: }
253: }
254:
255: #define INITLEN 10
256:
257: void lu_append(int *l,ent **a,int *l2,int k,int i,int mul)
258: {
259: int len;
260: ent *p;
261:
262: len = l[k];
263: if ( !len ) {
264: a[k] = p = (ent *)MALLOC_ATOMIC(INITLEN*sizeof(ent));
265: p[0].j = i; p[0].e = mul;
266: l[k] = 1; l2[k] = INITLEN;
267: } else {
268: if ( l2[k] == l[k] ) {
269: l2[k] *= 2;
270: a[k] = REALLOC(a[k],l2[k]*sizeof(ent));
271: }
272: p =a[k];
273: p[l[k]].j = i; p[l[k]].e = mul;
274: l[k]++;
275: }
276: }
277:
278: /* a[k] = a[k]-mul*a[i] */
279:
280: void lu_elim(int *l,ent **a,int k,int i,int mul,int mod)
281: {
282: ent *ak,*ai,*w;
283: int lk,li,j,m,p,q,r,s,t,j0;
284:
285: ak = a[k]; ai = a[i]; lk = l[k]; li = l[i];
286: w = (ent *)alloca((lk+li)*sizeof(ent));
287: p = 0; q = 0; j = 0;
288: mul = mod-mul;
289: while ( p < lk && q < li ) {
290: if ( ak[p].j > ai[q].j ) {
291: w[j] = ak[p]; j++; p++;
292: } else if ( ak[p].j < ai[q].j ) {
293: w[j].j = ai[q].j;
294: t = ai[q].e;
295: DMAR(t,mul,0,mod,r);
296: w[j].e = r;
297: j++; q++;
298: } else {
299: t = ai[q].e; s = ak[p].e;
300: DMAR(t,mul,s,mod,r);
301: if ( r ) {
302: w[j].j = ai[q].j; w[j].e = r; j++;
303: }
304: p++; q++;
305: }
306: }
307: if ( q == li )
308: while ( p < lk ) {
309: w[j] = ak[p]; j++; p++;
310: }
311: else if ( p == lk )
312: while ( q < li ) {
313: w[j].j = ai[q].j;
314: t = ai[q].e;
315: DMAR(t,mul,0,mod,r);
316: w[j].e = r;
317: j++; q++;
318: }
319: if ( j <= lk ) {
320: for ( m = 0; m < j; m++ ) ak[m] = w[m];
321: } else {
322: a[k] = ak = (ent *)MALLOC_ATOMIC(j*sizeof(ent));
323: for ( m = 0; m < j; m++ ) ak[m] = w[m];
324: }
325: l[k] = j;
326: }
327:
328: void solve_l(int *ll,ent **l,int n,int *rhs,int mod)
329: {
330: int j,k,s,len;
331: ent *p;
332:
333: for ( j = 0; j < n; j++ ) {
334: len = ll[j]; p = l[j];
335: for ( k = 0, s = 0; k < len; k++ )
336: s = dmar(p[k].e,rhs[p[k].j],s,mod);
337: rhs[j] -= s;
338: if ( rhs[j] < 0 ) rhs[j] += mod;
339: }
340: }
341:
342: void solve_u(int *ul,ent **u,int n,int *rhs,int mod)
343: {
344: int j,k,s,len,inv;
345: ent *p;
346:
347: for ( j = n-1; j >= 0; j-- ) {
348: len = ul[j]; p = u[j];
349: for ( k = 1, s = 0; k < len; k++ )
350: s = dmar(p[k].e,rhs[p[k].j],s,mod);
351: rhs[j] -= s;
352: if ( rhs[j] < 0 ) rhs[j] += mod;
353: inv = invm((unsigned int)p[0].e,mod);
354: rhs[j] = dmar(rhs[j],inv,0,mod);
355: }
356: }
357:
358: int comp_obj(Obj *a,Obj *b)
359: {
360: return arf_comp(CO,*a,*b);
361: }
362:
363: static FUNC generic_comp_obj_func;
364: static NODE generic_comp_obj_arg;
365: static NODE generic_comp_obj_option;
366:
367: int generic_comp_obj(Obj *a,Obj *b)
368: {
369: Q r;
370:
371: BDY(generic_comp_obj_arg)=(pointer)(*a);
372: BDY(NEXT(generic_comp_obj_arg))=(pointer)(*b);
373: r = (Q)bevalf_with_opts(generic_comp_obj_func,generic_comp_obj_arg,generic_comp_obj_option);
374: if ( !r )
375: return 0;
376: else
377: return sgnq(r)>0?1:-1;
378: }
379:
380:
381: void Pqsort(NODE arg,LIST *rp)
382: {
383: VECT vect;
384: NODE n,n1;
385: P p;
386: V v;
387: FUNC func;
388: int len,i;
389: pointer *a;
390: Obj t;
391:
392: t = ARG0(arg);
393: if (OID(t) == O_LIST) {
394: n = (NODE)BDY((LIST)t);
395: len = length(n);
396: MKVECT(vect,len);
397: for ( i = 0; i < len; i++, n = NEXT(n) ) {
398: BDY(vect)[i] = BDY(n);
399: }
400:
401: }else if (OID(t) != O_VECT) {
402: error("qsort : invalid argument");
403: }else {
404: vect = (VECT)t;
405: }
406: if ( argc(arg) == 1 )
407: qsort(BDY(vect),vect->len,sizeof(Obj),(int (*)(const void *,const void *))comp_obj);
408: else {
409: p = (P)ARG1(arg);
410: if ( !p || OID(p)!=2 )
411: error("qsort : invalid argument");
412: v = VR(p);
413: gen_searchf(NAME(v),&func);
414: if ( !func ) {
415: if ( (long)v->attr != V_SR )
416: error("qsort : no such function");
417: func = (FUNC)v->priv;
418: }
419: generic_comp_obj_func = func;
420: MKNODE(n,0,0); MKNODE(generic_comp_obj_arg,0,n);
421: generic_comp_obj_option = current_option;
422: qsort(BDY(vect),vect->len,sizeof(Obj),(int (*)(const void *,const void *))generic_comp_obj);
423: }
424: if (OID(t) == O_LIST) {
425: a = BDY(vect);
426: for ( i = len - 1, n = 0; i >= 0; i-- ) {
427: MKNODE(n1,a[i],n); n = n1;
428: }
429: MKLIST(*rp,n);
430: }else {
431: *rp = (LIST)vect;
432: }
433: }
434:
435: void PNBmul_gf2n(NODE arg,GF2N *rp)
436: {
437: GF2N a,b;
438: GF2MAT mat;
439: int n,w;
440: unsigned int *ab,*bb;
441: UP2 r;
442:
443: a = (GF2N)ARG0(arg);
444: b = (GF2N)ARG1(arg);
445: mat = (GF2MAT)ARG2(arg);
446: if ( !a || !b )
447: *rp = 0;
448: else {
449: n = mat->row;
450: w = (n+BSH-1)/BSH;
451:
452: ab = (unsigned int *)ALLOCA(w*sizeof(unsigned int));
453: bzero((char *)ab,w*sizeof(unsigned int));
454: bcopy(a->body->b,ab,(a->body->w)*sizeof(unsigned int));
455:
456: bb = (unsigned int *)ALLOCA(w*sizeof(unsigned int));
457: bzero((char *)bb,w*sizeof(unsigned int));
458: bcopy(b->body->b,bb,(b->body->w)*sizeof(unsigned int));
459:
460: NEWUP2(r,w);
461: bzero((char *)r->b,w*sizeof(unsigned int));
462: mul_nb(mat,ab,bb,r->b);
463: r->w = w;
464: _adjup2(r);
465: if ( !r->w )
466: *rp = 0;
467: else
468: MKGF2N(r,*rp);
469: }
470: }
471:
472: void Pmul_vect_mat_gf2n(NODE arg,GF2N *rp)
473: {
474: GF2N a;
475: GF2MAT mat;
476: int n,w;
477: unsigned int *b;
478: UP2 r;
479:
480: a = (GF2N)ARG0(arg);
481: mat = (GF2MAT)ARG1(arg);
482: if ( !a )
483: *rp = 0;
484: else {
485: n = mat->row;
486: w = (n+BSH-1)/BSH;
487: b = (unsigned int *)ALLOCA(w*sizeof(unsigned int));
488: bzero((char *)b,w*sizeof(unsigned int));
489: bcopy(a->body->b,b,(a->body->w)*sizeof(unsigned int));
490: NEWUP2(r,w);
491: bzero((char *)r->b,w*sizeof(unsigned int));
492: mulgf2vectmat(mat->row,b,mat->body,r->b);
493: r->w = w;
494: _adjup2(r);
495: if ( !r->w )
496: *rp = 0;
497: else {
498: MKGF2N(r,*rp);
499: }
500: }
501: }
502:
503: void Pbconvmat_gf2n(NODE arg,LIST *rp)
504: {
505: P p0,p1;
506: int to;
507: GF2MAT p01,p10;
508: GF2N root;
509: NODE n0,n1;
510:
511: p0 = (P)ARG0(arg);
512: p1 = (P)ARG1(arg);
513: to = ARG2(arg)?1:0;
514: if ( argc(arg) == 4 ) {
515: root = (GF2N)ARG3(arg);
516: compute_change_of_basis_matrix_with_root(p0,p1,to,root,&p01,&p10);
517: } else
518: compute_change_of_basis_matrix(p0,p1,to,&p01,&p10);
519: MKNODE(n1,p10,0); MKNODE(n0,p01,n1);
520: MKLIST(*rp,n0);
521: }
522:
523: void Pmulmat_gf2n(NODE arg,GF2MAT *rp)
524: {
525: GF2MAT m;
526:
527: if ( !compute_multiplication_matrix((P)ARG0(arg),&m) )
528: error("mulmat_gf2n : input is not a normal polynomial");
529: *rp = m;
530: }
531:
532: void Psepmat_destructive(NODE arg,LIST *rp)
533: {
534: MAT mat,mat1;
535: int i,j,row,col;
536: Z **a,**a1;
537: Z ent;
538: Z nm,mod,rem,quo;
539: int sgn;
540: NODE n0,n1;
541:
542: mat = (MAT)ARG0(arg); mod = (Z)ARG1(arg);
543: row = mat->row; col = mat->col;
544: MKMAT(mat1,row,col);
545: a = (Z **)mat->body; a1 = (Z **)mat1->body;
546: for ( i = 0; i < row; i++ )
547: for ( j = 0; j < col; j++ ) {
548: ent = a[i][j];
549: if ( !ent )
550: continue;
551: sgn = sgnz(ent);
552: absz(ent,&nm);
553: divqrz(nm,mod,&quo,&rem);
554: if ( sgn > 0 ) {
555: a[i][j] = rem;
556: a1[i][j] = quo;
557: } else {
558: chsgnz(rem,&a[i][j]);
559: chsgnz(quo,&a1[i][j]);
560: }
561: }
562: MKNODE(n1,mat1,0); MKNODE(n0,mat,n1);
563: MKLIST(*rp,n0);
564: }
565:
566: void Psepvect(NODE arg,VECT *rp)
567: {
1.2 ! noro 568: sepvect((VECT)ARG0(arg),ZTOS((Q)ARG1(arg)),rp);
1.1 noro 569: }
570:
571: void sepvect(VECT v,int d,VECT *rp)
572: {
573: int i,j,k,n,q,q1,r;
574: pointer *pv,*pw,*pu;
575: VECT w,u;
576:
577: n = v->len;
578: if ( d > n )
579: d = n;
580: q = n/d; r = n%d; q1 = q+1;
581: MKVECT(w,d); *rp = w;
582: pv = BDY(v); pw = BDY(w); k = 0;
583: for ( i = 0; i < r; i++ ) {
584: MKVECT(u,q1); pw[i] = (pointer)u;
585: for ( pu = BDY(u), j = 0; j < q1; j++, k++ )
586: pu[j] = pv[k];
587: }
588: for ( ; i < d; i++ ) {
589: MKVECT(u,q); pw[i] = (pointer)u;
590: for ( pu = BDY(u), j = 0; j < q; j++, k++ )
591: pu[j] = pv[k];
592: }
593: }
594:
595: void Pnewvect(NODE arg,VECT *rp)
596: {
597: int len,i,r;
598: VECT vect;
599: pointer *vb;
600: LIST list;
601: NODE tn;
602:
603: asir_assert(ARG0(arg),O_N,"newvect");
1.2 ! noro 604: len = ZTOS((Q)ARG0(arg));
1.1 noro 605: if ( len < 0 )
606: error("newvect : invalid size");
607: MKVECT(vect,len);
608: if ( argc(arg) == 2 ) {
609: list = (LIST)ARG1(arg);
610: asir_assert(list,O_LIST,"newvect");
611: #if 0
612: for ( r = 0, tn = BDY(list); tn; r++, tn = NEXT(tn) );
613: if ( r > len ) {
614: *rp = vect;
615: return;
616: }
617: #endif
618: for ( i = 0, tn = BDY(list), vb = BDY(vect); tn; i++, tn = NEXT(tn) )
619: vb[i] = (pointer)BDY(tn);
620: }
621: *rp = vect;
622: }
623:
624: void Pvect(NODE arg,VECT *rp) {
625: int len,i;
626: VECT vect;
627: pointer *vb;
628: NODE tn;
629:
630: if ( !arg ) {
631: *rp =0;
632: return;
633: }
634:
635: for (len = 0, tn = arg; tn; tn = NEXT(tn), len++);
636: if ( len == 1 ) {
637: if ( ARG0(arg) != 0 ) {
638: switch ( OID(ARG0(arg)) ) {
639: case O_VECT:
640: *rp = ARG0(arg);
641: return;
642: case O_LIST:
643: for ( len = 0, tn = ARG0(arg); tn; tn = NEXT(tn), len++ );
644: MKVECT(vect,len-1);
645: for ( i = 0, tn = BDY((LIST)ARG0(arg)), vb =BDY(vect);
646: tn; i++, tn = NEXT(tn) )
647: vb[i] = (pointer)BDY(tn);
648: *rp=vect;
649: return;
650: }
651: }
652: }
653: MKVECT(vect,len);
654: for ( i = 0, tn = arg, vb = BDY(vect); tn; i++, tn = NEXT(tn) )
655: vb[i] = (pointer)BDY(tn);
656: *rp = vect;
657: }
658:
659: void Pexponent_vector(NODE arg,DP *rp)
660: {
661: nodetod(arg,rp);
662: }
663:
664: void Pnewbytearray(NODE arg,BYTEARRAY *rp)
665: {
666: int len,i,r;
667: BYTEARRAY array;
668: unsigned char *vb;
669: char *str;
670: LIST list;
671: NODE tn;
672: int ac;
673: struct stat sbuf;
674: char *fname;
675: FILE *fp;
676:
677: ac = argc(arg);
678: if ( ac == 1 ) {
679: if ( !OID((Obj)ARG0(arg)) ) error("newbytearray : invalid argument");
680: switch ( OID((Obj)ARG0(arg)) ) {
681: case O_STR:
682: fname = BDY((STRING)ARG0(arg));
683: fp = fopen(fname,"rb");
684: if ( !fp ) error("newbytearray : fopen failed");
685: if ( stat(fname,&sbuf) < 0 )
686: error("newbytearray : stat failed");
687: len = sbuf.st_size;
688: MKBYTEARRAY(array,len);
689: fread(BDY(array),len,sizeof(char),fp);
690: break;
691: case O_N:
692: if ( !RATN(ARG0(arg)) )
693: error("newbytearray : invalid argument");
1.2 ! noro 694: len = ZTOS((Q)ARG0(arg));
1.1 noro 695: if ( len < 0 )
696: error("newbytearray : invalid size");
697: MKBYTEARRAY(array,len);
698: break;
699: default:
700: error("newbytearray : invalid argument");
701: }
702: } else if ( ac == 2 ) {
703: asir_assert(ARG0(arg),O_N,"newbytearray");
1.2 ! noro 704: len = ZTOS((Q)ARG0(arg));
1.1 noro 705: if ( len < 0 )
706: error("newbytearray : invalid size");
707: MKBYTEARRAY(array,len);
708: if ( !ARG1(arg) )
709: error("newbytearray : invalid initialization");
710: switch ( OID((Obj)ARG1(arg)) ) {
711: case O_LIST:
712: list = (LIST)ARG1(arg);
713: asir_assert(list,O_LIST,"newbytearray");
714: for ( r = 0, tn = BDY(list); tn; r++, tn = NEXT(tn) );
715: if ( r <= len ) {
716: for ( i = 0, tn = BDY(list), vb = BDY(array); tn;
717: i++, tn = NEXT(tn) )
1.2 ! noro 718: vb[i] = (unsigned char)ZTOS((Q)BDY(tn));
1.1 noro 719: }
720: break;
721: case O_STR:
722: str = BDY((STRING)ARG1(arg));
723: r = strlen(str);
724: if ( r <= len )
725: bcopy(str,BDY(array),r);
726: break;
727: default:
728: if ( !ARG1(arg) )
729: error("newbytearray : invalid initialization");
730: }
731: } else
732: error("newbytearray : invalid argument");
733: *rp = array;
734: }
735:
736: #define MEMORY_GETPOINT(a,len,x,y) (((a)[(len)*(y)+((x)>>3)])&(1<<((x)&7)))
737:
738: void Pmemoryplot_to_coord(NODE arg,LIST *rp)
739: {
740: int len,blen,y,i,j;
741: unsigned char *a;
742: NODE r0,r,n;
743: LIST l;
744: BYTEARRAY ba;
745: Z iq,jq;
746:
747: asir_assert(ARG0(arg),O_LIST,"memoryplot_to_coord");
748: arg = BDY((LIST)ARG0(arg));
1.2 ! noro 749: len = ZTOS((Q)ARG0(arg));
1.1 noro 750: blen = (len+7)/8;
1.2 ! noro 751: y = ZTOS((Q)ARG1(arg));
1.1 noro 752: ba = (BYTEARRAY)ARG2(arg); a = ba->body;
753: r0 = 0;
754: for ( j = 0; j < y; j++ )
755: for ( i = 0; i < len; i++ )
756: if ( MEMORY_GETPOINT(a,blen,i,j) ) {
757: NEXTNODE(r0,r);
1.2 ! noro 758: STOZ(i,iq); STOZ(j,jq);
1.1 noro 759: n = mknode(2,iq,jq);
760: MKLIST(l,n);
761: BDY(r) = l;
762: }
763: if ( r0 ) NEXT(r) = 0;
764: MKLIST(*rp,r0);
765: }
766:
767: void Pnewmat(NODE arg,MAT *rp)
768: {
769: int row,col;
770: int i,j,r,c;
771: NODE tn,sn;
772: MAT m;
773: pointer **mb;
774: LIST list;
775:
776: asir_assert(ARG0(arg),O_N,"newmat");
777: asir_assert(ARG1(arg),O_N,"newmat");
1.2 ! noro 778: row = ZTOS((Q)ARG0(arg)); col = ZTOS((Q)ARG1(arg));
1.1 noro 779: if ( row < 0 || col < 0 )
780: error("newmat : invalid size");
781: MKMAT(m,row,col);
782: if ( argc(arg) == 3 ) {
783: list = (LIST)ARG2(arg);
784: asir_assert(list,O_LIST,"newmat");
785: for ( r = 0, c = 0, tn = BDY(list); tn; r++, tn = NEXT(tn) ) {
786: for ( j = 0, sn = BDY((LIST)BDY(tn)); sn; j++, sn = NEXT(sn) );
787: c = MAX(c,j);
788: }
789: if ( (r > row) || (c > col) ) {
790: *rp = m;
791: return;
792: }
793: for ( i = 0, tn = BDY(list), mb = BDY(m); tn; i++, tn = NEXT(tn) ) {
794: asir_assert(BDY(tn),O_LIST,"newmat");
795: for ( j = 0, sn = BDY((LIST)BDY(tn)); sn; j++, sn = NEXT(sn) )
796: mb[i][j] = (pointer)BDY(sn);
797: }
798: }
799: *rp = m;
800: }
801:
802: void Pmat(NODE arg, MAT *rp)
803: {
804: int row,col;
805: int i;
806: MAT m;
807: pointer **mb;
808: pointer *ent;
809: NODE tn, sn;
810: VECT v;
811:
812: if ( !arg ) {
813: *rp =0;
814: return;
815: }
816:
817: for (row = 0, tn = arg; tn; tn = NEXT(tn), row++);
818: if ( row == 1 ) {
819: if ( OID(ARG0(arg)) == O_MAT ) {
820: *rp=ARG0(arg);
821: return;
822: } else if ( !(OID(ARG0(arg)) == O_LIST || OID(ARG0(arg)) == O_VECT)) {
823: error("mat : invalid argument");
824: }
825: }
826: if ( OID(ARG0(arg)) == O_VECT ) {
827: v = ARG0(arg);
828: col = v->len;
829: } else if ( OID(ARG0(arg)) == O_LIST ) {
830: for (col = 0, tn = BDY((LIST)ARG0(arg)); tn ; tn = NEXT(tn), col++);
831: } else {
832: error("mat : invalid argument");
833: }
834:
835: MKMAT(m,row,col);
836: for (row = 0, tn = arg, mb = BDY(m); tn; tn = NEXT(tn), row++) {
837: if ( BDY(tn) == 0 ) {
838: error("mat : invalid argument");
839: } else if ( OID(BDY(tn)) == O_VECT ) {
840: v = tn->body;
841: ent = BDY(v);
842: for (i = 0; i < v->len; i++ ) mb[row][i] = (Obj)ent[i];
843: } else if ( OID(BDY(tn)) == O_LIST ) {
844: for (col = 0, sn = BDY((LIST)BDY(tn)); sn; col++, sn = NEXT(sn) )
845: mb[row][col] = (pointer)BDY(sn);
846: } else {
847: error("mat : invalid argument");
848: }
849: }
850: *rp = m;
851: }
852:
853: void Pmatc(NODE arg, MAT *rp)
854: {
855: int row,col;
856: int i;
857: MAT m;
858: pointer **mb;
859: pointer *ent;
860: NODE tn, sn;
861: VECT v;
862:
863: if ( !arg ) {
864: *rp =0;
865: return;
866: }
867:
868: for (col = 0, tn = arg; tn; tn = NEXT(tn), col++);
869: if ( col == 1 ) {
870: if ( OID(ARG0(arg)) == O_MAT ) {
871: *rp=ARG0(arg);
872: return;
873: } else if ( !(OID(ARG0(arg)) == O_LIST || OID(ARG0(arg)) == O_VECT)) {
874: error("matc : invalid argument");
875: }
876: }
877: if ( OID(ARG0(arg)) == O_VECT ) {
878: v = ARG0(arg);
879: row = v->len;
880: } else if ( OID(ARG0(arg)) == O_LIST ) {
881: for (row = 0, tn = BDY((LIST)ARG0(arg)); tn ; tn = NEXT(tn), row++);
882: } else {
883: error("matc : invalid argument");
884: }
885:
886: MKMAT(m,row,col);
887: for (col = 0, tn = arg, mb = BDY(m); tn; tn = NEXT(tn), col++) {
888: if ( BDY(tn) == 0 ) {
889: error("matc : invalid argument");
890: } else if ( OID(BDY(tn)) == O_VECT ) {
891: v = tn->body;
892: ent = BDY(v);
893: for (i = 0; i < v->len; i++ ) mb[i][col] = (Obj)ent[i];
894: } else if ( OID(BDY(tn)) == O_LIST ) {
895: for (row = 0, sn = BDY((LIST)BDY(tn)); sn; row++, sn = NEXT(sn) )
896: mb[row][col] = (pointer)BDY(sn);
897: } else {
898: error("matc : invalid argument");
899: }
900: }
901: *rp = m;
902: }
903:
904: void Pvtol(NODE arg,LIST *rp)
905: {
906: NODE n,n1;
907: VECT v;
908: pointer *a;
909: int len,i;
910:
911: if ( OID(ARG0(arg)) == O_LIST ) {
912: *rp = ARG0(arg);
913: return;
914: }
915: asir_assert(ARG0(arg),O_VECT,"vtol");
916: v = (VECT)ARG0(arg); len = v->len; a = BDY(v);
917: for ( i = len - 1, n = 0; i >= 0; i-- ) {
918: MKNODE(n1,a[i],n); n = n1;
919: }
920: MKLIST(*rp,n);
921: }
922:
923: void Pltov(NODE arg,VECT *rp)
924: {
925: NODE n;
926: VECT v,v0;
927: int len,i;
928:
929: if ( OID(ARG0(arg)) == O_VECT ) {
930: v0 = (VECT)ARG0(arg); len = v0->len;
931: MKVECT(v,len);
932: for ( i = 0; i < len; i++ ) {
933: BDY(v)[i] = BDY(v0)[i];
934: }
935: *rp = v;
936: return;
937: }
938: asir_assert(ARG0(arg),O_LIST,"ltov");
939: n = (NODE)BDY((LIST)ARG0(arg));
940: len = length(n);
941: MKVECT(v,len);
942: for ( i = 0; i < len; i++, n = NEXT(n) )
943: BDY(v)[i] = BDY(n);
944: *rp = v;
945: }
946:
947: /* XXX it seems that remainder is not used */
948:
949: void Premainder(NODE arg,Obj *rp)
950: {
951: Obj a;
952: VECT v,w;
953: MAT m,l;
954: pointer *vb,*wb;
955: pointer **mb,**lb;
956: int id,i,j,n,row,col,t,smd,sgn;
957: Z md,q,r;
958:
959: a = (Obj)ARG0(arg); md = (Z)ARG1(arg);
960: if ( !a )
961: *rp = 0;
962: else {
963: id = OID(a);
964: switch ( id ) {
965: case O_N:
966: case O_P:
967: cmp(md,(P)a,(P *)rp); break;
968: case O_VECT:
969: /* strage spec */
1.2 ! noro 970: smd = ZTOS(md);
1.1 noro 971: v = (VECT)a; n = v->len; vb = v->body;
972: MKVECT(w,n); wb = w->body;
973: for ( i = 0; i < n; i++ ) {
974: if ( sgnz(vb[i]) > 0 )
975: remz(vb[i],md,(Z *)&wb[i]);
976: else {
977: absz(vb[i],&q);
978: remz(q,md,&r);
979: chsgnz(r,(Z *)&wb[i]);
980: }
981: }
982: *rp = (Obj)w;
983: break;
984: case O_MAT:
985: m = (MAT)a; row = m->row; col = m->col; mb = m->body;
986: MKMAT(l,row,col); lb = l->body;
987: for ( i = 0; i < row; i++ )
988: for ( j = 0, vb = mb[i], wb = lb[i]; j < col; j++ )
989: cmp(md,(P)vb[j],(P *)&wb[j]);
990: *rp = (Obj)l;
991: break;
992: default:
993: error("remainder : invalid argument");
994: }
995: }
996: }
997:
998: void Psremainder(NODE arg,Obj *rp)
999: {
1000: Obj a;
1001: VECT v,w;
1002: MAT m,l;
1003: pointer *vb,*wb;
1004: pointer **mb,**lb;
1005: int id,i,j,n,row,col,t,smd,sgn;
1006: Z md,q,r;
1007:
1008: a = (Obj)ARG0(arg); md = (Z)ARG1(arg);
1009: if ( !a )
1010: *rp = 0;
1011: else {
1012: id = OID(a);
1013: switch ( id ) {
1014: case O_N:
1015: case O_P:
1016: cmp(md,(P)a,(P *)rp); break;
1017: case O_VECT:
1.2 ! noro 1018: smd = ZTOS(md);
1.1 noro 1019: v = (VECT)a; n = v->len; vb = v->body;
1020: MKVECT(w,n); wb = w->body;
1021: for ( i = 0; i < n; i++ )
1022: remz(vb[i],md,(Z *)&wb[i]);
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: }
1037: }
1038:
1039: void Psize(NODE arg,LIST *rp)
1040: {
1041:
1042: int n,m;
1043: Z 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;
1.2 ! noro 1052: STOZ(n,q); MKNODE(t,q,0);
1.1 noro 1053: break;
1054: case O_MAT:
1055: n = ((MAT)ARG0(arg))->row; m = ((MAT)ARG0(arg))->col;
1.2 ! noro 1056: STOZ(m,q); MKNODE(s,q,0); STOZ(n,q); MKNODE(t,q,s);
1.1 noro 1057: break;
1058: case O_IMAT:
1059: n = ((IMAT)ARG0(arg))->row; m = ((IMAT)ARG0(arg))->col;
1.2 ! noro 1060: STOZ(m,q); MKNODE(s,q,0); STOZ(n,q); MKNODE(t,q,s);
1.1 noro 1061: break;
1062: default:
1063: error("size : invalid argument"); break;
1064: }
1065: }
1066: MKLIST(*rp,t);
1067: }
1068:
1069: void Pdet(NODE arg,P *rp)
1070: {
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 {
1.2 ! noro 1083: n = m->row; mod = ZTOS((Q)ARG1(arg)); mat = (P **)BDY(m);
1.1 noro 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: }
1091: }
1092:
1093: void Pinvmat(NODE arg,LIST *rp)
1094: {
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 {
1.2 ! noro 1112: n = m->row; mod = ZTOS((Q)ARG1(arg)); mat = (P **)BDY(m);
1.1 noro 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]);
1117: #if 0
1118: detmp(CO,mod,w,n,&d);
1119: mptop(d,rp);
1120: #else
1121: error("not implemented yet");
1122: #endif
1123: }
1124: }
1125:
1126: /*
1127: input : a row x col matrix A
1128: A[I] <-> A[I][0]*x_0+A[I][1]*x_1+...
1129:
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]}+...
1136: */
1137:
1138: void Pgeneric_gauss_elim(NODE arg,LIST *rp)
1139: {
1140: NODE n0,opt,p;
1141: MAT m,nm;
1142: int *ri,*ci;
1143: VECT rind,cind;
1144: Z 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++ ) {
1.2 ! noro 1172: STOZ(ri[i],q);
1.1 noro 1173: BDY(rind)[i] = (pointer)q;
1174: }
1175: for ( i = 0; i < t; i++ ) {
1.2 ! noro 1176: STOZ(ci[i],q);
1.1 noro 1177: BDY(cind)[i] = (pointer)q;
1178: }
1179: n0 = mknode(4,nm,dn,rind,cind);
1180: MKLIST(*rp,n0);
1181: }
1182:
1183: int indep_rows_mod(int **mat0,int row,int col,int md,int *rowstat);
1184:
1185: void Pindep_rows_mod(NODE arg,VECT *rp)
1186: {
1187: MAT m,mat;
1188: VECT rind;
1189: Z **tmat;
1190: int **wmat,**row0;
1191: Z *rib;
1192: int *rowstat,*p;
1193: Z 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");
1.2 ! noro 1198: m = (MAT)ARG0(arg); md = ZTOS((Z)ARG1(arg));
1.1 noro 1199: row = m->row; col = m->col; tmat = (Z **)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: wmat[i][j] = remqi((Q)tmat[i][j],md);
1209: rank = indep_rows_mod(wmat,row,col,md,rowstat);
1210:
1211: MKVECT(rind,rank);
1212: rib = (Z *)rind->body;
1213: for ( j = 0; j < rank; j++ ) {
1.2 ! noro 1214: STOZ(rowstat[j],rib[j]);
1.1 noro 1215: }
1216: *rp = rind;
1217: }
1218:
1219: /*
1220: input : a row x col matrix A
1221: A[I] <-> A[I][0]*x_0+A[I][1]*x_1+...
1222:
1223: output : [B,R,C]
1224: B : a rank(A) x col-rank(A) matrix
1225: R : a vector of length rank(A)
1226: C : a vector of length col-rank(A)
1227: RN : a vector of length rank(A) indicating useful rows
1228:
1229: B[I] <-> x_{R[I]}+B[I][0]x_{C[0]}+B[I][1]x_{C[1]}+...
1230: */
1231:
1232: void Pgeneric_gauss_elim_mod(NODE arg,LIST *rp)
1233: {
1234: NODE n0;
1235: MAT m,mat;
1236: VECT rind,cind,rnum;
1237: Z **tmat;
1238: int **wmat,**row0;
1239: Z *rib,*cib,*rnb;
1240: int *colstat,*p;
1241: Z q;
1242: int md,i,j,k,l,row,col,t,rank;
1243:
1244: asir_assert(ARG0(arg),O_MAT,"generic_gauss_elim_mod");
1245: asir_assert(ARG1(arg),O_N,"generic_gauss_elim_mod");
1.2 ! noro 1246: m = (MAT)ARG0(arg); md = ZTOS((Z)ARG1(arg));
1.1 noro 1247: row = m->row; col = m->col; tmat = (Z **)m->body;
1248: wmat = (int **)almat(row,col);
1249:
1250: row0 = (int **)ALLOCA(row*sizeof(int *));
1251: for ( i = 0; i < row; i++ ) row0[i] = wmat[i];
1252:
1253: colstat = (int *)MALLOC_ATOMIC(col*sizeof(int));
1254: for ( i = 0; i < row; i++ )
1255: for ( j = 0; j < col; j++ )
1256: wmat[i][j] = remqi((Q)tmat[i][j],md);
1257: rank = generic_gauss_elim_mod(wmat,row,col,md,colstat);
1258:
1259: MKVECT(rnum,rank);
1260: rnb = (Z *)rnum->body;
1261: for ( i = 0; i < rank; i++ )
1262: for ( j = 0, p = wmat[i]; j < row; j++ )
1263: if ( p == row0[j] )
1.2 ! noro 1264: STOZ(j,rnb[i]);
1.1 noro 1265:
1266: MKMAT(mat,rank,col-rank);
1267: tmat = (Z **)mat->body;
1268: for ( i = 0; i < rank; i++ )
1269: for ( j = k = 0; j < col; j++ )
1270: if ( !colstat[j] ) {
1.2 ! noro 1271: UTOZ(wmat[i][j],tmat[i][k]); k++;
1.1 noro 1272: }
1273:
1274: MKVECT(rind,rank);
1275: MKVECT(cind,col-rank);
1276: rib = (Z *)rind->body; cib = (Z *)cind->body;
1277: for ( j = k = l = 0; j < col; j++ )
1278: if ( colstat[j] ) {
1.2 ! noro 1279: STOZ(j,rib[k]); k++;
1.1 noro 1280: } else {
1.2 ! noro 1281: STOZ(j,cib[l]); l++;
1.1 noro 1282: }
1283: n0 = mknode(4,mat,rind,cind,rnum);
1284: MKLIST(*rp,n0);
1285: }
1286:
1287: void Pleqm(NODE arg,VECT *rp)
1288: {
1289: MAT m;
1290: VECT vect;
1291: pointer **mat;
1292: Z *v;
1293: Z q;
1294: int **wmat;
1295: int md,i,j,row,col,t,n,status;
1296:
1297: asir_assert(ARG0(arg),O_MAT,"leqm");
1298: asir_assert(ARG1(arg),O_N,"leqm");
1.2 ! noro 1299: m = (MAT)ARG0(arg); md = ZTOS((Z)ARG1(arg));
1.1 noro 1300: row = m->row; col = m->col; mat = m->body;
1301: wmat = (int **)almat(row,col);
1302: for ( i = 0; i < row; i++ )
1303: for ( j = 0; j < col; j++ )
1304: wmat[i][j] = remqi((Q)mat[i][j],md);
1305: status = gauss_elim_mod(wmat,row,col,md);
1306: if ( status < 0 )
1307: *rp = 0;
1308: else if ( status > 0 )
1309: *rp = (VECT)ONE;
1310: else {
1311: n = col - 1;
1312: MKVECT(vect,n);
1313: for ( i = 0, v = (Z *)vect->body; i < n; i++ ) {
1.2 ! noro 1314: t = (md-wmat[i][n])%md; STOZ(t,v[i]);
1.1 noro 1315: }
1316: *rp = vect;
1317: }
1318: }
1319:
1320: int gauss_elim_mod(int **mat,int row,int col,int md)
1321: {
1322: int i,j,k,inv,a,n;
1323: int *t,*pivot;
1324:
1325: n = col - 1;
1326: for ( j = 0; j < n; j++ ) {
1327: for ( i = j; i < row && !mat[i][j]; i++ );
1328: if ( i == row )
1329: return 1;
1330: if ( i != j ) {
1331: t = mat[i]; mat[i] = mat[j]; mat[j] = t;
1332: }
1333: pivot = mat[j];
1334: inv = invm(pivot[j],md);
1335: for ( k = j; k <= n; k++ ) {
1336: /* pivot[k] = dmar(pivot[k],inv,0,md); */
1337: DMAR(pivot[k],inv,0,md,pivot[k])
1338: }
1339: for ( i = 0; i < row; i++ ) {
1340: t = mat[i];
1341: if ( i != j && (a = t[j]) )
1342: for ( k = j, a = md - a; k <= n; k++ ) {
1343: unsigned int tk;
1344: /* t[k] = dmar(pivot[k],a,t[k],md); */
1345: DMAR(pivot[k],a,t[k],md,tk)
1346: t[k] = tk;
1347: }
1348: }
1349: }
1350: for ( i = n; i < row && !mat[i][n]; i++ );
1351: if ( i == row )
1352: return 0;
1353: else
1354: return -1;
1355: }
1356:
1357: struct oEGT eg_mod,eg_elim,eg_elim1,eg_elim2,eg_chrem,eg_gschk,eg_intrat,eg_symb;
1358: struct oEGT eg_conv;
1359:
1360: #if 0
1361: void lu_dec_cr(MAT mat,MAT lu,Q *dn,int **perm);
1362:
1363: /* XXX broken */
1364: void lu_dec_cr(MAT mat,MAT lu,Q *dn,int **perm)
1365: {
1366: Q **a0,**b;
1367: Q *aiq;
1368: N **a;
1369: N *ai;
1370: Q q,q1,dn2,a1,q0,bik;
1371: MAT m;
1372: unsigned int md;
1373: int n,ind,i,j,rank,t,inv,t1,ret,min,k;
1374: int **w;
1375: int *wi,*rinfo0,*rinfo;
1376: N m1,m2,m3,u,s;
1377:
1378: a0 = (Q **)mat->body;
1379: n = mat->row;
1380: if ( n != mat->col )
1381: error("lu_dec_cr : non-square matrix");
1382: w = (int **)almat(n,n);
1383: MKMAT(m,n,n);
1384: a = (N **)m->body;
1385: UTON(1,m1);
1386: rinfo0 = 0;
1387: ind = 0;
1388: while ( 1 ) {
1389: md = get_lprime(ind);
1390: /* mat mod md */
1391: for ( i = 0; i < n; i++ )
1392: for ( j = 0, aiq = a0[i], wi = w[i]; j < n; j++ )
1393: if ( q = aiq[j] ) {
1394: t = rem(NM(q),md);
1395: if ( t && SGN(q) < 0 )
1396: t = (md - t) % md;
1397: wi[j] = t;
1398: } else
1399: wi[j] = 0;
1400:
1401: if ( !lu_mod((unsigned int **)w,n,md,&rinfo) ) continue;
1402: printf("."); fflush(stdout);
1403: if ( !rinfo0 )
1404: *perm = rinfo0 = rinfo;
1405: else {
1406: for ( i = 0; i < n; i++ )
1407: if ( rinfo[i] != rinfo0[i] ) break;
1408: if ( i < n ) continue;
1409: }
1410: if ( UNIN(m1) ) {
1411: for ( i = 0; i < n; i++ )
1412: for ( j = 0, ai = a[i], wi = w[i]; j < n; j++ ) {
1413: UTON(wi[j],u); ai[j] = u;
1414: }
1415: UTON(md,m1);
1416: } else {
1417: inv = invm(rem(m1,md),md);
1418: UTON(md,m2); muln(m1,m2,&m3);
1419: for ( i = 0; i < n; i++ )
1420: for ( j = 0, ai = a[i], wi = w[i]; j < n; j++ )
1421: if ( ai[i] ) {
1422: /* f3 = f1+m1*(m1 mod m2)^(-1)*(f2 - f1 mod m2) */
1423: t = rem(ai[j],md);
1424: if ( wi[j] >= t )
1425: t = wi[j]-t;
1426: else
1427: t = md-(t-wi[j]);
1428: DMAR(t,inv,0,md,t1)
1429: UTON(t1,u);
1430: muln(m1,u,&s);
1431: addn(ai[j],s,&u); ai[j] = u;
1432: } else if ( wi[j] ) {
1433: /* f3 = m1*(m1 mod m2)^(-1)*f2 */
1434: DMAR(wi[j],inv,0,md,t)
1435: UTON(t,u);
1436: muln(m1,u,&s); ai[j] = s;
1437: }
1438: m1 = m3;
1439: }
1440: if ( (++ind%8) == 0 ) {
1441: ret = intmtoratm(m,m1,lu,dn);
1442: if ( ret ) {
1443: b = (Q **)lu->body;
1444: mulq(*dn,*dn,&dn2);
1445: for ( i = 0; i < n; i++ ) {
1446: for ( j = 0; j < n; j++ ) {
1447: q = 0;
1448: min = MIN(i,j);
1449: for ( k = 0; k <= min; k++ ) {
1450: bik = k==i ? *dn : b[i][k];
1451: mulq(bik,b[k][j],&q0);
1452: addq(q,q0,&q1); q = q1;
1453: }
1454: mulq(a0[rinfo0[i]][j],dn2,&q1);
1455: if ( cmpq(q,q1) ) break;
1456: }
1457: if ( j < n ) break;
1458: }
1459: if ( i == n )
1460: return;
1461: }
1462: }
1463: }
1464: }
1465: #endif
1466:
1467:
1468: int f4_nocheck;
1469:
1470: #define ONE_STEP1 if ( zzz = *s ) { DMAR(zzz,hc,*tj,md,*tj) } tj++; s++;
1471:
1472: void reduce_reducers_mod(int **mat,int row,int col,int md)
1473: {
1474: int i,j,k,l,hc,zzz;
1475: int *t,*s,*tj,*ind;
1476:
1477: /* reduce the reducers */
1478: ind = (int *)ALLOCA(row*sizeof(int));
1479: for ( i = 0; i < row; i++ ) {
1480: t = mat[i];
1481: for ( j = 0; j < col && !t[j]; j++ );
1482: /* register the position of the head term */
1483: ind[i] = j;
1484: for ( l = i-1; l >= 0; l-- ) {
1485: /* reduce mat[i] by mat[l] */
1486: if ( hc = t[ind[l]] ) {
1487: /* mat[i] = mat[i]-hc*mat[l] */
1488: j = ind[l];
1489: s = mat[l]+j;
1490: tj = t+j;
1491: hc = md-hc;
1492: k = col-j;
1493: for ( ; k >= 64; k -= 64 ) {
1494: ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
1495: ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
1496: ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
1497: ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
1498: ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
1499: ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
1500: ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
1501: ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
1502: ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
1503: ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
1504: ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
1505: ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
1506: ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
1507: ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
1508: ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
1509: ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
1510: }
1511: for ( ; k > 0; k-- ) {
1512: if ( zzz = *s ) { DMAR(zzz,hc,*tj,md,*tj) } tj++; s++;
1513: }
1514: }
1515: }
1516: }
1517: }
1518:
1519: /*
1520: mat[i] : reducers (i=0,...,nred-1)
1521: spolys (i=nred,...,row-1)
1522: mat[0] < mat[1] < ... < mat[nred-1] w.r.t the term order
1523: 1. reduce the reducers
1524: 2. reduce spolys by the reduced reducers
1525: */
1526:
1527: void pre_reduce_mod(int **mat,int row,int col,int nred,int md)
1528: {
1529: int i,j,k,l,hc,inv;
1530: int *t,*s,*tk,*ind;
1531:
1532: #if 1
1533: /* reduce the reducers */
1534: ind = (int *)ALLOCA(row*sizeof(int));
1535: for ( i = 0; i < nred; i++ ) {
1536: /* make mat[i] monic and mat[i] by mat[0],...,mat[i-1] */
1537: t = mat[i];
1538: for ( j = 0; j < col && !t[j]; j++ );
1539: /* register the position of the head term */
1540: ind[i] = j;
1541: inv = invm(t[j],md);
1542: for ( k = j; k < col; k++ )
1543: if ( t[k] )
1544: DMAR(t[k],inv,0,md,t[k])
1545: for ( l = i-1; l >= 0; l-- ) {
1546: /* reduce mat[i] by mat[l] */
1547: if ( hc = t[ind[l]] ) {
1548: /* mat[i] = mat[i]-hc*mat[l] */
1549: for ( k = ind[l], hc = md-hc, s = mat[l]+k, tk = t+k;
1550: k < col; k++, tk++, s++ )
1551: if ( *s )
1552: DMAR(*s,hc,*tk,md,*tk)
1553: }
1554: }
1555: }
1556: /* reduce the spolys */
1557: for ( i = nred; i < row; i++ ) {
1558: t = mat[i];
1559: for ( l = nred-1; l >= 0; l-- ) {
1560: /* reduce mat[i] by mat[l] */
1561: if ( hc = t[ind[l]] ) {
1562: /* mat[i] = mat[i]-hc*mat[l] */
1563: for ( k = ind[l], hc = md-hc, s = mat[l]+k, tk = t+k;
1564: k < col; k++, tk++, s++ )
1565: if ( *s )
1566: DMAR(*s,hc,*tk,md,*tk)
1567: }
1568: }
1569: }
1570: #endif
1571: }
1572: /*
1573: mat[i] : reducers (i=0,...,nred-1)
1574: mat[0] < mat[1] < ... < mat[nred-1] w.r.t the term order
1575: */
1576:
1577: void reduce_sp_by_red_mod(int *sp,int **redmat,int *ind,int nred,int col,int md)
1578: {
1579: int i,j,k,hc,zzz;
1580: int *s,*tj;
1581:
1582: /* reduce the spolys by redmat */
1583: for ( i = nred-1; i >= 0; i-- ) {
1584: /* reduce sp by redmat[i] */
1585: if ( hc = sp[ind[i]] ) {
1586: /* sp = sp-hc*redmat[i] */
1587: j = ind[i];
1588: hc = md-hc;
1589: s = redmat[i]+j;
1590: tj = sp+j;
1591: for ( k = col-j; k > 0; k-- ) {
1592: if ( zzz = *s ) { DMAR(zzz,hc,*tj,md,*tj) } tj++; s++;
1593: }
1594: }
1595: }
1596: }
1597:
1598: /*
1599: mat[i] : compressed reducers (i=0,...,nred-1)
1600: mat[0] < mat[1] < ... < mat[nred-1] w.r.t the term order
1601: */
1602:
1603: void red_by_compress(int m,unsigned int *p,unsigned int *r,
1604: unsigned int *ri,unsigned int hc,int len)
1605: {
1606: unsigned int up,lo;
1607: unsigned int dmy;
1608: unsigned int *pj;
1609:
1610: p[*ri] = 0; r++; ri++;
1611: for ( len--; len; len--, r++, ri++ ) {
1612: pj = p+ *ri;
1613: DMA(*r,hc,*pj,up,lo);
1614: if ( up ) {
1615: DSAB(m,up,lo,dmy,*pj);
1616: } else
1617: *pj = lo;
1618: }
1619: }
1620:
1621: /* p -= hc*r */
1622:
1623: void red_by_vect(int m,unsigned int *p,unsigned int *r,unsigned int hc,int len)
1624: {
1625: unsigned int up,lo,dmy;
1626:
1627: *p++ = 0; r++; len--;
1628: for ( ; len; len--, r++, p++ )
1629: if ( *r ) {
1630: DMA(*r,hc,*p,up,lo);
1631: if ( up ) {
1632: DSAB(m,up,lo,dmy,*p);
1633: } else
1634: *p = lo;
1635: }
1636: }
1637:
1638: #if defined(__GNUC__) && SIZEOF_LONG==8
1639: /* 64bit vector += UNIT vector(normalized) */
1640:
1641: void red_by_vect64(int m, U64 *p,unsigned int *c,U64 *r,unsigned int hc,int len)
1642: {
1643: U64 t;
1644:
1645: /* (p[0],c[0]) is normalized */
1646: *p++ = 0; *c++ = 0; r++; len--;
1647: for ( ; len; len--, r++, p++, c++ )
1648: if ( *r ) {
1649: t = (*p)+(*r)*hc;
1650: if ( t < *p ) (*c)++;
1651: *p = t;
1652: }
1653: }
1654: #endif
1655:
1656: void red_by_vect_sf(int m,unsigned int *p,unsigned int *r,unsigned int hc,int len)
1657: {
1658: *p++ = 0; r++; len--;
1659: for ( ; len; len--, r++, p++ )
1660: if ( *r )
1661: *p = _addsf(_mulsf(*r,hc),*p);
1662: }
1663:
1664: extern Z current_mod_lf;
1665: extern int current_mod_lf_size;
1666:
1667: void red_by_vect_lf(mpz_t *p,mpz_t *r,mpz_t hc,int len)
1668: {
1669: mpz_set_ui(*p++,0); r++; len--;
1670: for ( ; len; len--, r++, p++ ) {
1671: mpz_addmul(*p,*r,hc);
1672: #if 0
1673: if ( mpz_size(*p) > current_mod_lf_size )
1674: mpz_mod(*p,*p,BDY(current_mod_lf));
1675: #endif
1676: }
1677: }
1678:
1679:
1680: extern unsigned int **psca;
1681:
1682: void reduce_sp_by_red_mod_compress (int *sp,CDP *redmat,int *ind,
1683: int nred,int col,int md)
1684: {
1685: int i,len;
1686: CDP ri;
1687: unsigned int hc;
1688: unsigned int *usp;
1689:
1690: usp = (unsigned int *)sp;
1691: /* reduce the spolys by redmat */
1692: for ( i = nred-1; i >= 0; i-- ) {
1693: /* reduce sp by redmat[i] */
1694: usp[ind[i]] %= md;
1695: if ( hc = usp[ind[i]] ) {
1696: /* sp = sp-hc*redmat[i] */
1697: hc = md-hc;
1698: ri = redmat[i];
1699: len = ri->len;
1700: red_by_compress(md,usp,psca[ri->psindex],ri->body,hc,len);
1701: }
1702: }
1703: for ( i = 0; i < col; i++ )
1704: if ( usp[i] >= (unsigned int)md )
1705: usp[i] %= md;
1706: }
1707:
1708: #define ONE_STEP2 if ( zzz = *pk ) { DMAR(zzz,a,*tk,md,*tk) } pk++; tk++;
1709:
1710: int generic_gauss_elim_mod(int **mat0,int row,int col,int md,int *colstat)
1711: {
1712: int i,j,k,l,inv,a,rank;
1713: unsigned int *t,*pivot,*pk;
1714: unsigned int **mat;
1715:
1716: mat = (unsigned int **)mat0;
1717: for ( rank = 0, j = 0; j < col; j++ ) {
1718: for ( i = rank; i < row; i++ )
1719: mat[i][j] %= md;
1720: for ( i = rank; i < row; i++ )
1721: if ( mat[i][j] )
1722: break;
1723: if ( i == row ) {
1724: colstat[j] = 0;
1725: continue;
1726: } else
1727: colstat[j] = 1;
1728: if ( i != rank ) {
1729: t = mat[i]; mat[i] = mat[rank]; mat[rank] = t;
1730: }
1731: pivot = mat[rank];
1732: inv = invm(pivot[j],md);
1733: for ( k = j, pk = pivot+k; k < col; k++, pk++ )
1734: if ( *pk ) {
1735: if ( *pk >= (unsigned int)md )
1736: *pk %= md;
1737: DMAR(*pk,inv,0,md,*pk)
1738: }
1739: for ( i = rank+1; i < row; i++ ) {
1740: t = mat[i];
1741: if ( a = t[j] )
1742: red_by_vect(md,t+j,pivot+j,md-a,col-j);
1743: }
1744: rank++;
1745: }
1746: for ( j = col-1, l = rank-1; j >= 0; j-- )
1747: if ( colstat[j] ) {
1748: pivot = mat[l];
1749: for ( i = 0; i < l; i++ ) {
1750: t = mat[i];
1751: t[j] %= md;
1752: if ( a = t[j] )
1753: red_by_vect(md,t+j,pivot+j,md-a,col-j);
1754: }
1755: l--;
1756: }
1757: for ( j = 0, l = 0; l < rank; j++ )
1758: if ( colstat[j] ) {
1759: t = mat[l];
1760: for ( k = j; k < col; k++ )
1761: if ( t[k] >= (unsigned int)md )
1762: t[k] %= md;
1763: l++;
1764: }
1765: return rank;
1766: }
1767:
1768: int generic_gauss_elim_mod2(int **mat0,int row,int col,int md,int *colstat,int *rowstat)
1769: {
1770: int i,j,k,l,inv,a,rank;
1771: unsigned int *t,*pivot,*pk;
1772: unsigned int **mat;
1773:
1774: for ( i = 0; i < row; i++ ) rowstat[i] = i;
1775: mat = (unsigned int **)mat0;
1776: for ( rank = 0, j = 0; j < col; j++ ) {
1777: for ( i = rank; i < row; i++ )
1778: mat[i][j] %= md;
1779: for ( i = rank; i < row; i++ )
1780: if ( mat[i][j] )
1781: break;
1782: if ( i == row ) {
1783: colstat[j] = 0;
1784: continue;
1785: } else
1786: colstat[j] = 1;
1787: if ( i != rank ) {
1788: t = mat[i]; mat[i] = mat[rank]; mat[rank] = t;
1789: k = rowstat[i]; rowstat[i] = rowstat[rank]; rowstat[rank] = k;
1790: }
1791: pivot = mat[rank];
1792: inv = invm(pivot[j],md);
1793: for ( k = j, pk = pivot+k; k < col; k++, pk++ )
1794: if ( *pk ) {
1795: if ( *pk >= (unsigned int)md )
1796: *pk %= md;
1797: DMAR(*pk,inv,0,md,*pk)
1798: }
1799: for ( i = rank+1; i < row; i++ ) {
1800: t = mat[i];
1801: if ( a = t[j] )
1802: red_by_vect(md,t+j,pivot+j,md-a,col-j);
1803: }
1804: rank++;
1805: }
1806: for ( j = col-1, l = rank-1; j >= 0; j-- )
1807: if ( colstat[j] ) {
1808: pivot = mat[l];
1809: for ( i = 0; i < l; i++ ) {
1810: t = mat[i];
1811: t[j] %= md;
1812: if ( a = t[j] )
1813: red_by_vect(md,t+j,pivot+j,md-a,col-j);
1814: }
1815: l--;
1816: }
1817: for ( j = 0, l = 0; l < rank; j++ )
1818: if ( colstat[j] ) {
1819: t = mat[l];
1820: for ( k = j; k < col; k++ )
1821: if ( t[k] >= (unsigned int)md )
1822: t[k] %= md;
1823: l++;
1824: }
1825: return rank;
1826: }
1827:
1828: int indep_rows_mod(int **mat0,int row,int col,int md,int *rowstat)
1829: {
1830: int i,j,k,l,inv,a,rank;
1831: unsigned int *t,*pivot,*pk;
1832: unsigned int **mat;
1833:
1834: for ( i = 0; i < row; i++ ) rowstat[i] = i;
1835: mat = (unsigned int **)mat0;
1836: for ( rank = 0, j = 0; j < col; j++ ) {
1837: for ( i = rank; i < row; i++ )
1838: mat[i][j] %= md;
1839: for ( i = rank; i < row; i++ )
1840: if ( mat[i][j] )
1841: break;
1842: if ( i == row ) continue;
1843: if ( i != rank ) {
1844: t = mat[i]; mat[i] = mat[rank]; mat[rank] = t;
1845: k = rowstat[i]; rowstat[i] = rowstat[rank]; rowstat[rank] = k;
1846: }
1847: pivot = mat[rank];
1848: inv = invm(pivot[j],md);
1849: for ( k = j, pk = pivot+k; k < col; k++, pk++ )
1850: if ( *pk ) {
1851: if ( *pk >= (unsigned int)md )
1852: *pk %= md;
1853: DMAR(*pk,inv,0,md,*pk)
1854: }
1855: for ( i = rank+1; i < row; i++ ) {
1856: t = mat[i];
1857: if ( a = t[j] )
1858: red_by_vect(md,t+j,pivot+j,md-a,col-j);
1859: }
1860: rank++;
1861: }
1862: return rank;
1863: }
1864:
1865: int generic_gauss_elim_sf(int **mat0,int row,int col,int md,int *colstat)
1866: {
1867: int i,j,k,l,inv,a,rank;
1868: unsigned int *t,*pivot,*pk;
1869: unsigned int **mat;
1870:
1871: mat = (unsigned int **)mat0;
1872: for ( rank = 0, j = 0; j < col; j++ ) {
1873: for ( i = rank; i < row; i++ )
1874: if ( mat[i][j] )
1875: break;
1876: if ( i == row ) {
1877: colstat[j] = 0;
1878: continue;
1879: } else
1880: colstat[j] = 1;
1881: if ( i != rank ) {
1882: t = mat[i]; mat[i] = mat[rank]; mat[rank] = t;
1883: }
1884: pivot = mat[rank];
1885: inv = _invsf(pivot[j]);
1886: for ( k = j, pk = pivot+k; k < col; k++, pk++ )
1887: if ( *pk )
1888: *pk = _mulsf(*pk,inv);
1889: for ( i = rank+1; i < row; i++ ) {
1890: t = mat[i];
1891: if ( a = t[j] )
1892: red_by_vect_sf(md,t+j,pivot+j,_chsgnsf(a),col-j);
1893: }
1894: rank++;
1895: }
1896: for ( j = col-1, l = rank-1; j >= 0; j-- )
1897: if ( colstat[j] ) {
1898: pivot = mat[l];
1899: for ( i = 0; i < l; i++ ) {
1900: t = mat[i];
1901: if ( a = t[j] )
1902: red_by_vect_sf(md,t+j,pivot+j,_chsgnsf(a),col-j);
1903: }
1904: l--;
1905: }
1906: return rank;
1907: }
1908:
1909: /* LU decomposition; a[i][i] = 1/U[i][i] */
1910:
1911: int lu_gfmmat(GFMMAT mat,unsigned int md,int *perm)
1912: {
1913: int row,col;
1914: int i,j,k;
1915: unsigned int *t,*pivot;
1916: unsigned int **a;
1917: unsigned int inv,m;
1918:
1919: row = mat->row; col = mat->col;
1920: a = mat->body;
1921: bzero(perm,row*sizeof(int));
1922:
1923: for ( i = 0; i < row; i++ )
1924: perm[i] = i;
1925: for ( k = 0; k < col; k++ ) {
1926: for ( i = k; i < row && !a[i][k]; i++ );
1927: if ( i == row )
1928: return 0;
1929: if ( i != k ) {
1930: j = perm[i]; perm[i] = perm[k]; perm[k] = j;
1931: t = a[i]; a[i] = a[k]; a[k] = t;
1932: }
1933: pivot = a[k];
1934: pivot[k] = inv = invm(pivot[k],md);
1935: for ( i = k+1; i < row; i++ ) {
1936: t = a[i];
1937: if ( m = t[k] ) {
1938: DMAR(inv,m,0,md,t[k])
1939: for ( j = k+1, m = md - t[k]; j < col; j++ )
1940: if ( pivot[j] ) {
1941: unsigned int tj;
1942:
1943: DMAR(m,pivot[j],t[j],md,tj)
1944: t[j] = tj;
1945: }
1946: }
1947: }
1948: }
1949: return 1;
1950: }
1951:
1952: /*
1953: Input
1954: a: a row x col matrix
1955: md : a modulus
1956:
1957: Output:
1958: return : d = the rank of mat
1959: a[0..(d-1)][0..(d-1)] : LU decomposition (a[i][i] = 1/U[i][i])
1960: rinfo: array of length row
1961: cinfo: array of length col
1962: i-th row in new a <-> rinfo[i]-th row in old a
1963: cinfo[j]=1 <=> j-th column is contained in the LU decomp.
1964: */
1965:
1966: int find_lhs_and_lu_mod(unsigned int **a,int row,int col,
1967: unsigned int md,int **rinfo,int **cinfo)
1968: {
1969: int i,j,k,d;
1970: int *rp,*cp;
1971: unsigned int *t,*pivot;
1972: unsigned int inv,m;
1973:
1974: *rinfo = rp = (int *)MALLOC_ATOMIC(row*sizeof(int));
1975: *cinfo = cp = (int *)MALLOC_ATOMIC(col*sizeof(int));
1976: for ( i = 0; i < row; i++ )
1977: rp[i] = i;
1978: for ( k = 0, d = 0; k < col; k++ ) {
1979: for ( i = d; i < row && !a[i][k]; i++ );
1980: if ( i == row ) {
1981: cp[k] = 0;
1982: continue;
1983: } else
1984: cp[k] = 1;
1985: if ( i != d ) {
1986: j = rp[i]; rp[i] = rp[d]; rp[d] = j;
1987: t = a[i]; a[i] = a[d]; a[d] = t;
1988: }
1989: pivot = a[d];
1990: pivot[k] = inv = invm(pivot[k],md);
1991: for ( i = d+1; i < row; i++ ) {
1992: t = a[i];
1993: if ( m = t[k] ) {
1994: DMAR(inv,m,0,md,t[k])
1995: for ( j = k+1, m = md - t[k]; j < col; j++ )
1996: if ( pivot[j] ) {
1997: unsigned int tj;
1998: DMAR(m,pivot[j],t[j],md,tj)
1999: t[j] = tj;
2000: }
2001: }
2002: }
2003: d++;
2004: }
2005: return d;
2006: }
2007:
2008: int lu_mod(unsigned int **a,int n,unsigned int md,int **rinfo)
2009: {
2010: int i,j,k;
2011: int *rp;
2012: unsigned int *t,*pivot;
2013: unsigned int inv,m;
2014:
2015: *rinfo = rp = (int *)MALLOC_ATOMIC(n*sizeof(int));
2016: for ( i = 0; i < n; i++ ) rp[i] = i;
2017: for ( k = 0; k < n; k++ ) {
2018: for ( i = k; i < n && !a[i][k]; i++ );
2019: if ( i == n ) return 0;
2020: if ( i != k ) {
2021: j = rp[i]; rp[i] = rp[k]; rp[k] = j;
2022: t = a[i]; a[i] = a[k]; a[k] = t;
2023: }
2024: pivot = a[k];
2025: inv = invm(pivot[k],md);
2026: for ( i = k+1; i < n; i++ ) {
2027: t = a[i];
2028: if ( m = t[k] ) {
2029: DMAR(inv,m,0,md,t[k])
2030: for ( j = k+1, m = md - t[k]; j < n; j++ )
2031: if ( pivot[j] ) {
2032: unsigned int tj;
2033: DMAR(m,pivot[j],t[j],md,tj)
2034: t[j] = tj;
2035: }
2036: }
2037: }
2038: }
2039: return 1;
2040: }
2041:
2042: /*
2043: Input
2044: a : n x n matrix; a result of LU-decomposition
2045: md : modulus
2046: b : n x l matrix
2047: Output
2048: b = a^(-1)b
2049: */
2050:
2051: void solve_by_lu_mod(int **a,int n,int md,int **b,int l,int normalize)
2052: {
2053: unsigned int *y,*c;
2054: int i,j,k;
2055: unsigned int t,m,m2;
2056:
2057: y = (int *)MALLOC_ATOMIC(n*sizeof(int));
2058: c = (int *)MALLOC_ATOMIC(n*sizeof(int));
2059: m2 = md>>1;
2060: for ( k = 0; k < l; k++ ) {
2061: /* copy b[.][k] to c */
2062: for ( i = 0; i < n; i++ )
2063: c[i] = (unsigned int)b[i][k];
2064: /* solve Ly=c */
2065: for ( i = 0; i < n; i++ ) {
2066: for ( t = c[i], j = 0; j < i; j++ )
2067: if ( a[i][j] ) {
2068: m = md - a[i][j];
2069: DMAR(m,y[j],t,md,t)
2070: }
2071: y[i] = t;
2072: }
2073: /* solve Uc=y */
2074: for ( i = n-1; i >= 0; i-- ) {
2075: for ( t = y[i], j =i+1; j < n; j++ )
2076: if ( a[i][j] ) {
2077: m = md - a[i][j];
2078: DMAR(m,c[j],t,md,t)
2079: }
2080: /* a[i][i] = 1/U[i][i] */
2081: DMAR(t,a[i][i],0,md,c[i])
2082: }
2083: /* copy c to b[.][k] with normalization */
2084: if ( normalize )
2085: for ( i = 0; i < n; i++ )
2086: b[i][k] = (int)(c[i]>m2 ? c[i]-md : c[i]);
2087: else
2088: for ( i = 0; i < n; i++ )
2089: b[i][k] = c[i];
2090: }
2091: }
2092:
2093: void Pleqm1(NODE arg,VECT *rp)
2094: {
2095: MAT m;
2096: VECT vect;
2097: pointer **mat;
2098: Z *v;
2099: Z q;
2100: int **wmat;
2101: int md,i,j,row,col,t,n,status;
2102:
2103: asir_assert(ARG0(arg),O_MAT,"leqm1");
2104: asir_assert(ARG1(arg),O_N,"leqm1");
1.2 ! noro 2105: m = (MAT)ARG0(arg); md = ZTOS((Z)ARG1(arg));
1.1 noro 2106: row = m->row; col = m->col; mat = m->body;
2107: wmat = (int **)almat(row,col);
2108: for ( i = 0; i < row; i++ )
2109: for ( j = 0; j < col; j++ )
2110: wmat[i][j] = remqi((Q)mat[i][j],md);
2111: status = gauss_elim_mod1(wmat,row,col,md);
2112: if ( status < 0 )
2113: *rp = 0;
2114: else if ( status > 0 )
2115: *rp = (VECT)ONE;
2116: else {
2117: n = col - 1;
2118: MKVECT(vect,n);
2119: for ( i = 0, v = (Z *)vect->body; i < n; i++ ) {
1.2 ! noro 2120: t = (md-wmat[i][n])%md; STOZ(t,v[i]);
1.1 noro 2121: }
2122: *rp = vect;
2123: }
2124: }
2125:
2126: int gauss_elim_mod1(int **mat,int row,int col,int md)
2127: {
2128: int i,j,k,inv,a,n;
2129: int *t,*pivot;
2130:
2131: n = col - 1;
2132: for ( j = 0; j < n; j++ ) {
2133: for ( i = j; i < row && !mat[i][j]; i++ );
2134: if ( i == row )
2135: return 1;
2136: if ( i != j ) {
2137: t = mat[i]; mat[i] = mat[j]; mat[j] = t;
2138: }
2139: pivot = mat[j];
2140: inv = invm(pivot[j],md);
2141: for ( k = j; k <= n; k++ )
2142: pivot[k] = dmar(pivot[k],inv,0,md);
2143: for ( i = j+1; i < row; i++ ) {
2144: t = mat[i];
2145: if ( i != j && (a = t[j]) )
2146: for ( k = j, a = md - a; k <= n; k++ )
2147: t[k] = dmar(pivot[k],a,t[k],md);
2148: }
2149: }
2150: for ( i = n; i < row && !mat[i][n]; i++ );
2151: if ( i == row ) {
2152: for ( j = n-1; j >= 0; j-- ) {
2153: for ( i = j-1, a = (md-mat[j][n])%md; i >= 0; i-- ) {
2154: mat[i][n] = dmar(mat[i][j],a,mat[i][n],md);
2155: mat[i][j] = 0;
2156: }
2157: }
2158: return 0;
2159: } else
2160: return -1;
2161: }
2162:
2163: void Pgeninvm(NODE arg,LIST *rp)
2164: {
2165: MAT m;
2166: pointer **mat;
2167: Z **tmat;
2168: Z q;
2169: unsigned int **wmat;
2170: int md,i,j,row,col,t,status;
2171: MAT mat1,mat2;
2172: NODE node1,node2;
2173:
2174: asir_assert(ARG0(arg),O_MAT,"leqm1");
2175: asir_assert(ARG1(arg),O_N,"leqm1");
1.2 ! noro 2176: m = (MAT)ARG0(arg); md = ZTOS((Q)ARG1(arg));
1.1 noro 2177: row = m->row; col = m->col; mat = m->body;
2178: wmat = (unsigned int **)almat(row,col+row);
2179: for ( i = 0; i < row; i++ ) {
2180: bzero((char *)wmat[i],(col+row)*sizeof(int));
2181: for ( j = 0; j < col; j++ )
2182: wmat[i][j] = remqi((Q)mat[i][j],md);
2183: }
2184: status = gauss_elim_geninv_mod(wmat,row,col,md);
2185: if ( status > 0 )
2186: *rp = 0;
2187: else {
2188: MKMAT(mat1,col,row); MKMAT(mat2,row-col,row);
2189: for ( i = 0, tmat = (Z **)mat1->body; i < col; i++ )
2190: for ( j = 0; j < row; j++ )
1.2 ! noro 2191: UTOZ(wmat[i][j+col],tmat[i][j]);
1.1 noro 2192: for ( tmat = (Z **)mat2->body; i < row; i++ )
2193: for ( j = 0; j < row; j++ )
1.2 ! noro 2194: UTOZ(wmat[i][j+col],tmat[i-col][j]);
1.1 noro 2195: MKNODE(node2,mat2,0); MKNODE(node1,mat1,node2); MKLIST(*rp,node1);
2196: }
2197: }
2198:
2199: int gauss_elim_geninv_mod(unsigned int **mat,int row,int col,int md)
2200: {
2201: int i,j,k,inv,a,n,m;
2202: unsigned int *t,*pivot;
2203:
2204: n = col; m = row+col;
2205: for ( j = 0; j < n; j++ ) {
2206: for ( i = j; i < row && !mat[i][j]; i++ );
2207: if ( i == row )
2208: return 1;
2209: if ( i != j ) {
2210: t = mat[i]; mat[i] = mat[j]; mat[j] = t;
2211: }
2212: pivot = mat[j];
2213: inv = invm(pivot[j],md);
2214: for ( k = j; k < m; k++ )
2215: pivot[k] = dmar(pivot[k],inv,0,md);
2216: for ( i = j+1; i < row; i++ ) {
2217: t = mat[i];
2218: if ( a = t[j] )
2219: for ( k = j, a = md - a; k < m; k++ )
2220: t[k] = dmar(pivot[k],a,t[k],md);
2221: }
2222: }
2223: for ( j = n-1; j >= 0; j-- ) {
2224: pivot = mat[j];
2225: for ( i = j-1; i >= 0; i-- ) {
2226: t = mat[i];
2227: if ( a = t[j] )
2228: for ( k = j, a = md - a; k < m; k++ )
2229: t[k] = dmar(pivot[k],a,t[k],md);
2230: }
2231: }
2232: return 0;
2233: }
2234:
2235: void Psolve_by_lu_gfmmat(NODE arg,VECT *rp)
2236: {
2237: GFMMAT lu;
2238: Z *perm,*rhs,*v;
2239: int n,i;
2240: unsigned int md;
2241: unsigned int *b,*sol;
2242: VECT r;
2243:
2244: lu = (GFMMAT)ARG0(arg);
2245: perm = (Z *)BDY((VECT)ARG1(arg));
2246: rhs = (Z *)BDY((VECT)ARG2(arg));
1.2 ! noro 2247: md = (unsigned int)ZTOS((Z)ARG3(arg));
1.1 noro 2248: n = lu->col;
2249: b = (unsigned int *)MALLOC_ATOMIC(n*sizeof(int));
2250: sol = (unsigned int *)MALLOC_ATOMIC(n*sizeof(int));
2251: for ( i = 0; i < n; i++ )
1.2 ! noro 2252: b[i] = ZTOS(rhs[ZTOS(perm[i])]);
1.1 noro 2253: solve_by_lu_gfmmat(lu,md,b,sol);
2254: MKVECT(r,n);
2255: for ( i = 0, v = (Z *)r->body; i < n; i++ )
1.2 ! noro 2256: UTOZ(sol[i],v[i]);
1.1 noro 2257: *rp = r;
2258: }
2259:
2260: void solve_by_lu_gfmmat(GFMMAT lu,unsigned int md,
2261: unsigned int *b,unsigned int *x)
2262: {
2263: int n;
2264: unsigned int **a;
2265: unsigned int *y;
2266: int i,j;
2267: unsigned int t,m;
2268:
2269: n = lu->col;
2270: a = lu->body;
2271: y = (unsigned int *)MALLOC_ATOMIC(n*sizeof(int));
2272: /* solve Ly=b */
2273: for ( i = 0; i < n; i++ ) {
2274: for ( t = b[i], j = 0; j < i; j++ )
2275: if ( a[i][j] ) {
2276: m = md - a[i][j];
2277: DMAR(m,y[j],t,md,t)
2278: }
2279: y[i] = t;
2280: }
2281: /* solve Ux=y */
2282: for ( i = n-1; i >= 0; i-- ) {
2283: for ( t = y[i], j =i+1; j < n; j++ )
2284: if ( a[i][j] ) {
2285: m = md - a[i][j];
2286: DMAR(m,x[j],t,md,t)
2287: }
2288: /* a[i][i] = 1/U[i][i] */
2289: DMAR(t,a[i][i],0,md,x[i])
2290: }
2291: }
2292:
2293: #if 0
2294: void Plu_mat(NODE arg,LIST *rp)
2295: {
2296: MAT m,lu;
2297: Q dn;
2298: Q *v;
2299: int n,i;
2300: int *iperm;
2301: VECT perm;
2302: NODE n0;
2303:
2304: asir_assert(ARG0(arg),O_MAT,"lu_mat");
2305: m = (MAT)ARG0(arg);
2306: n = m->row;
2307: MKMAT(lu,n,n);
2308: lu_dec_cr(m,lu,&dn,&iperm);
2309: MKVECT(perm,n);
2310: for ( i = 0, v = (Q *)perm->body; i < n; i++ )
1.2 ! noro 2311: STOZ(iperm[i],v[i]);
1.1 noro 2312: n0 = mknode(3,lu,dn,perm);
2313: MKLIST(*rp,n0);
2314: }
2315: #endif
2316:
2317: void Plu_gfmmat(NODE arg,LIST *rp)
2318: {
2319: MAT m;
2320: GFMMAT mm;
2321: unsigned int md;
2322: int i,row,col,status;
2323: int *iperm;
2324: Z *v;
2325: VECT perm;
2326: NODE n0;
2327:
2328: asir_assert(ARG0(arg),O_MAT,"lu_gfmmat");
2329: asir_assert(ARG1(arg),O_N,"lu_gfmmat");
1.2 ! noro 2330: m = (MAT)ARG0(arg); md = (unsigned int)ZTOS((Q)ARG1(arg));
1.1 noro 2331: mat_to_gfmmat(m,md,&mm);
2332: row = m->row;
2333: col = m->col;
2334: iperm = (int *)MALLOC_ATOMIC(row*sizeof(int));
2335: status = lu_gfmmat(mm,md,iperm);
2336: if ( !status )
2337: n0 = 0;
2338: else {
2339: MKVECT(perm,row);
2340: for ( i = 0, v = (Z *)perm->body; i < row; i++ )
1.2 ! noro 2341: STOZ(iperm[i],v[i]);
1.1 noro 2342: n0 = mknode(2,mm,perm);
2343: }
2344: MKLIST(*rp,n0);
2345: }
2346:
2347: void Pmat_to_gfmmat(NODE arg,GFMMAT *rp)
2348: {
2349: MAT m;
2350: unsigned int md;
2351:
2352: asir_assert(ARG0(arg),O_MAT,"mat_to_gfmmat");
2353: asir_assert(ARG1(arg),O_N,"mat_to_gfmmat");
1.2 ! noro 2354: m = (MAT)ARG0(arg); md = (unsigned int)ZTOS((Q)ARG1(arg));
1.1 noro 2355: mat_to_gfmmat(m,md,rp);
2356: }
2357:
2358: void mat_to_gfmmat(MAT m,unsigned int md,GFMMAT *rp)
2359: {
2360: unsigned int **wmat;
2361: unsigned int t;
2362: Z **mat;
2363: Z q;
2364: int i,j,row,col;
2365:
2366: row = m->row; col = m->col; mat = (Z **)m->body;
2367: wmat = (unsigned int **)almat(row,col);
2368: for ( i = 0; i < row; i++ ) {
2369: bzero((char *)wmat[i],col*sizeof(unsigned int));
2370: for ( j = 0; j < col; j++ )
2371: wmat[i][j] = remqi((Q)mat[i][j],md);
2372: }
2373: TOGFMMAT(row,col,wmat,*rp);
2374: }
2375:
2376: void Pgeninvm_swap(NODE arg,LIST *rp)
2377: {
2378: MAT m;
2379: pointer **mat;
2380: Z **tmat;
2381: Z *tvect;
2382: Z q;
2383: unsigned int **wmat,**invmat;
2384: int *index;
2385: unsigned int t,md;
2386: int i,j,row,col,status;
2387: MAT mat1;
2388: VECT vect1;
2389: NODE node1,node2;
2390:
2391: asir_assert(ARG0(arg),O_MAT,"geninvm_swap");
2392: asir_assert(ARG1(arg),O_N,"geninvm_swap");
1.2 ! noro 2393: m = (MAT)ARG0(arg); md = ZTOS((Z)ARG1(arg));
1.1 noro 2394: row = m->row; col = m->col; mat = m->body;
2395: wmat = (unsigned int **)almat(row,col+row);
2396: for ( i = 0; i < row; i++ ) {
2397: bzero((char *)wmat[i],(col+row)*sizeof(int));
2398: for ( j = 0; j < col; j++ )
2399: wmat[i][j] = remqi((Q)mat[i][j],md);
2400: wmat[i][col+i] = 1;
2401: }
2402: status = gauss_elim_geninv_mod_swap(wmat,row,col,md,&invmat,&index);
2403: if ( status > 0 )
2404: *rp = 0;
2405: else {
2406: MKMAT(mat1,col,col);
2407: for ( i = 0, tmat = (Z **)mat1->body; i < col; i++ )
2408: for ( j = 0; j < col; j++ )
1.2 ! noro 2409: UTOZ(invmat[i][j],tmat[i][j]);
1.1 noro 2410: MKVECT(vect1,row);
2411: for ( i = 0, tvect = (Z *)vect1->body; i < row; i++ )
1.2 ! noro 2412: STOZ(index[i],tvect[i]);
1.1 noro 2413: MKNODE(node2,vect1,0); MKNODE(node1,mat1,node2); MKLIST(*rp,node1);
2414: }
2415: }
2416:
2417: int gauss_elim_geninv_mod_swap(unsigned int **mat,int row,int col,unsigned int md,
2418: unsigned int ***invmatp,int **indexp)
2419: {
2420: int i,j,k,inv,a,n,m;
2421: unsigned int *t,*pivot,*s;
2422: int *index;
2423: unsigned int **invmat;
2424:
2425: n = col; m = row+col;
2426: *indexp = index = (int *)MALLOC_ATOMIC(row*sizeof(int));
2427: for ( i = 0; i < row; i++ )
2428: index[i] = i;
2429: for ( j = 0; j < n; j++ ) {
2430: for ( i = j; i < row && !mat[i][j]; i++ );
2431: if ( i == row ) {
2432: *indexp = 0; *invmatp = 0; return 1;
2433: }
2434: if ( i != j ) {
2435: t = mat[i]; mat[i] = mat[j]; mat[j] = t;
2436: k = index[i]; index[i] = index[j]; index[j] = k;
2437: }
2438: pivot = mat[j];
2439: inv = (unsigned int)invm(pivot[j],md);
2440: for ( k = j; k < m; k++ )
2441: if ( pivot[k] )
2442: pivot[k] = (unsigned int)dmar(pivot[k],inv,0,md);
2443: for ( i = j+1; i < row; i++ ) {
2444: t = mat[i];
2445: if ( a = t[j] )
2446: for ( k = j, a = md - a; k < m; k++ )
2447: if ( pivot[k] )
2448: t[k] = dmar(pivot[k],a,t[k],md);
2449: }
2450: }
2451: for ( j = n-1; j >= 0; j-- ) {
2452: pivot = mat[j];
2453: for ( i = j-1; i >= 0; i-- ) {
2454: t = mat[i];
2455: if ( a = t[j] )
2456: for ( k = j, a = md - a; k < m; k++ )
2457: if ( pivot[k] )
2458: t[k] = dmar(pivot[k],a,t[k],md);
2459: }
2460: }
2461: *invmatp = invmat = (unsigned int **)almat(col,col);
2462: for ( i = 0; i < col; i++ )
2463: for ( j = 0, s = invmat[i], t = mat[i]; j < col; j++ )
2464: s[j] = t[col+index[j]];
2465: return 0;
2466: }
2467:
2468: int gauss_elim_geninv_sf_swap(int **mat,int row,int col, int ***invmatp,int **indexp);
2469:
2470: void Pgeninv_sf_swap(NODE arg,LIST *rp)
2471: {
2472: MAT m;
2473: GFS **mat,**tmat;
2474: Z *tvect;
2475: GFS q;
2476: int **wmat,**invmat;
2477: int *index;
2478: unsigned int t;
2479: int i,j,row,col,status;
2480: MAT mat1;
2481: VECT vect1;
2482: NODE node1,node2;
2483:
2484: asir_assert(ARG0(arg),O_MAT,"geninv_sf_swap");
2485: m = (MAT)ARG0(arg);
2486: row = m->row; col = m->col; mat = (GFS **)m->body;
2487: wmat = (int **)almat(row,col+row);
2488: for ( i = 0; i < row; i++ ) {
2489: bzero((char *)wmat[i],(col+row)*sizeof(int));
2490: for ( j = 0; j < col; j++ )
2491: if ( q = (GFS)mat[i][j] )
2492: wmat[i][j] = FTOIF(CONT(q));
2493: wmat[i][col+i] = _onesf();
2494: }
2495: status = gauss_elim_geninv_sf_swap(wmat,row,col,&invmat,&index);
2496: if ( status > 0 )
2497: *rp = 0;
2498: else {
2499: MKMAT(mat1,col,col);
2500: for ( i = 0, tmat = (GFS **)mat1->body; i < col; i++ )
2501: for ( j = 0; j < col; j++ )
2502: if ( t = invmat[i][j] ) {
2503: MKGFS(IFTOF(t),tmat[i][j]);
2504: }
2505: MKVECT(vect1,row);
2506: for ( i = 0, tvect = (Z *)vect1->body; i < row; i++ )
1.2 ! noro 2507: STOZ(index[i],tvect[i]);
1.1 noro 2508: MKNODE(node2,vect1,0); MKNODE(node1,mat1,node2); MKLIST(*rp,node1);
2509: }
2510: }
2511:
2512: int gauss_elim_geninv_sf_swap(int **mat,int row,int col, int ***invmatp,int **indexp)
2513: {
2514: int i,j,k,inv,a,n,m,u;
2515: int *t,*pivot,*s;
2516: int *index;
2517: int **invmat;
2518:
2519: n = col; m = row+col;
2520: *indexp = index = (int *)MALLOC_ATOMIC(row*sizeof(int));
2521: for ( i = 0; i < row; i++ )
2522: index[i] = i;
2523: for ( j = 0; j < n; j++ ) {
2524: for ( i = j; i < row && !mat[i][j]; i++ );
2525: if ( i == row ) {
2526: *indexp = 0; *invmatp = 0; return 1;
2527: }
2528: if ( i != j ) {
2529: t = mat[i]; mat[i] = mat[j]; mat[j] = t;
2530: k = index[i]; index[i] = index[j]; index[j] = k;
2531: }
2532: pivot = mat[j];
2533: inv = _invsf(pivot[j]);
2534: for ( k = j; k < m; k++ )
2535: if ( pivot[k] )
2536: pivot[k] = _mulsf(pivot[k],inv);
2537: for ( i = j+1; i < row; i++ ) {
2538: t = mat[i];
2539: if ( a = t[j] )
2540: for ( k = j, a = _chsgnsf(a); k < m; k++ )
2541: if ( pivot[k] ) {
2542: u = _mulsf(pivot[k],a);
2543: t[k] = _addsf(u,t[k]);
2544: }
2545: }
2546: }
2547: for ( j = n-1; j >= 0; j-- ) {
2548: pivot = mat[j];
2549: for ( i = j-1; i >= 0; i-- ) {
2550: t = mat[i];
2551: if ( a = t[j] )
2552: for ( k = j, a = _chsgnsf(a); k < m; k++ )
2553: if ( pivot[k] ) {
2554: u = _mulsf(pivot[k],a);
2555: t[k] = _addsf(u,t[k]);
2556: }
2557: }
2558: }
2559: *invmatp = invmat = (int **)almat(col,col);
2560: for ( i = 0; i < col; i++ )
2561: for ( j = 0, s = invmat[i], t = mat[i]; j < col; j++ )
2562: s[j] = t[col+index[j]];
2563: return 0;
2564: }
2565:
2566: void inner_product_int(Z *a,Z *b,int n,Z *r)
2567: {
2568: int i;
2569: Z t;
2570:
2571: t = 0;
2572: for ( i = 0; i < n; i++ )
2573: muladdtoz(a[i],b[i],&t);
2574: *r = t;
2575: }
2576:
2577: void Pmul_mat_vect_int(NODE arg,VECT *rp)
2578: {
2579: MAT mat;
2580: VECT vect,r;
2581: int row,col,i;
2582:
2583: mat = (MAT)ARG0(arg);
2584: vect = (VECT)ARG1(arg);
2585: row = mat->row;
2586: col = mat->col;
2587: MKVECT(r,row);
2588: for ( i = 0; i < row; i++ ) {
2589: inner_product_int((Z *)mat->body[i],(Z *)vect->body,col,(Z *)&r->body[i]);
2590: }
2591: *rp = r;
2592: }
2593:
2594: void Pnbpoly_up2(NODE arg,GF2N *rp)
2595: {
2596: int m,type,ret;
2597: UP2 r;
2598:
1.2 ! noro 2599: m = ZTOS((Z)ARG0(arg));
! 2600: type = ZTOS((Z)ARG1(arg));
1.1 noro 2601: ret = generate_ONB_polynomial(&r,m,type);
2602: if ( ret == 0 )
2603: MKGF2N(r,*rp);
2604: else
2605: *rp = 0;
2606: }
2607:
2608: void Px962_irredpoly_up2(NODE arg,GF2N *rp)
2609: {
2610: int m,ret,w;
2611: GF2N prev;
2612: UP2 r;
2613:
1.2 ! noro 2614: m = ZTOS((Q)ARG0(arg));
1.1 noro 2615: prev = (GF2N)ARG1(arg);
2616: if ( !prev ) {
2617: w = (m>>5)+1; NEWUP2(r,w); r->w = 0;
2618: bzero((char *)r->b,w*sizeof(unsigned int));
2619: } else {
2620: r = prev->body;
2621: if ( degup2(r) != m ) {
2622: w = (m>>5)+1; NEWUP2(r,w); r->w = 0;
2623: bzero((char *)r->b,w*sizeof(unsigned int));
2624: }
2625: }
2626: ret = _generate_irreducible_polynomial(r,m);
2627: if ( ret == 0 )
2628: MKGF2N(r,*rp);
2629: else
2630: *rp = 0;
2631: }
2632:
2633: void Pirredpoly_up2(NODE arg,GF2N *rp)
2634: {
2635: int m,ret,w;
2636: GF2N prev;
2637: UP2 r;
2638:
1.2 ! noro 2639: m = ZTOS((Q)ARG0(arg));
1.1 noro 2640: prev = (GF2N)ARG1(arg);
2641: if ( !prev ) {
2642: w = (m>>5)+1; NEWUP2(r,w); r->w = 0;
2643: bzero((char *)r->b,w*sizeof(unsigned int));
2644: } else {
2645: r = prev->body;
2646: if ( degup2(r) != m ) {
2647: w = (m>>5)+1; NEWUP2(r,w); r->w = 0;
2648: bzero((char *)r->b,w*sizeof(unsigned int));
2649: }
2650: }
2651: ret = _generate_good_irreducible_polynomial(r,m);
2652: if ( ret == 0 )
2653: MKGF2N(r,*rp);
2654: else
2655: *rp = 0;
2656: }
2657:
2658: void Pmat_swap_row_destructive(NODE arg, MAT *m)
2659: {
2660: int i1,i2;
2661: pointer *t;
2662: MAT mat;
2663:
2664: asir_assert(ARG0(arg),O_MAT,"mat_swap_row_destructive");
2665: asir_assert(ARG1(arg),O_N,"mat_swap_row_destructive");
2666: asir_assert(ARG2(arg),O_N,"mat_swap_row_destructive");
2667: mat = (MAT)ARG0(arg);
1.2 ! noro 2668: i1 = ZTOS((Q)ARG1(arg));
! 2669: i2 = ZTOS((Q)ARG2(arg));
1.1 noro 2670: if ( i1 < 0 || i2 < 0 || i1 >= mat->row || i2 >= mat->row )
2671: error("mat_swap_row_destructive : Out of range");
2672: t = mat->body[i1];
2673: mat->body[i1] = mat->body[i2];
2674: mat->body[i2] = t;
2675: *m = mat;
2676: }
2677:
2678: void Pmat_swap_col_destructive(NODE arg, MAT *m)
2679: {
2680: int j1,j2,i,n;
2681: pointer *mi;
2682: pointer t;
2683: MAT mat;
2684:
2685: asir_assert(ARG0(arg),O_MAT,"mat_swap_col_destructive");
2686: asir_assert(ARG1(arg),O_N,"mat_swap_col_destructive");
2687: asir_assert(ARG2(arg),O_N,"mat_swap_col_destructive");
2688: mat = (MAT)ARG0(arg);
1.2 ! noro 2689: j1 = ZTOS((Q)ARG1(arg));
! 2690: j2 = ZTOS((Q)ARG2(arg));
1.1 noro 2691: if ( j1 < 0 || j2 < 0 || j1 >= mat->col || j2 >= mat->col )
2692: error("mat_swap_col_destructive : Out of range");
2693: n = mat->row;
2694: for ( i = 0; i < n; i++ ) {
2695: mi = mat->body[i];
2696: t = mi[j1]; mi[j1] = mi[j2]; mi[j2] = t;
2697: }
2698: *m = mat;
2699: }
2700: /*
2701: * f = type 'type' normal polynomial of degree m if exists
2702: * IEEE P1363 A.7.2
2703: *
2704: * return value : 0 --- exists
2705: * 1 --- does not exist
2706: * -1 --- failure (memory allocation error)
2707: */
2708:
2709: int generate_ONB_polynomial(UP2 *rp,int m,int type)
2710: {
2711: int i,r;
2712: int w;
2713: UP2 f,f0,f1,f2,t;
2714:
2715: w = (m>>5)+1;
2716: switch ( type ) {
2717: case 1:
2718: if ( !TypeT_NB_check(m,1) ) return 1;
2719: NEWUP2(f,w); *rp = f; f->w = w;
2720: /* set all the bits */
2721: for ( i = 0; i < w; i++ )
2722: f->b[i] = 0xffffffff;
2723: /* mask the top word if necessary */
2724: if ( r = (m+1)&31 )
2725: f->b[w-1] &= (1<<r)-1;
2726: return 0;
2727: break;
2728: case 2:
2729: if ( !TypeT_NB_check(m,2) ) return 1;
2730: NEWUP2(f,w); *rp = f;
2731: W_NEWUP2(f0,w);
2732: W_NEWUP2(f1,w);
2733: W_NEWUP2(f2,w);
2734:
2735: /* recursion for genrating Type II normal polynomial */
2736:
2737: /* f0 = 1, f1 = t+1 */
2738: f0->w = 1; f0->b[0] = 1;
2739: f1->w = 1; f1->b[0] = 3;
2740: for ( i = 2; i <= m; i++ ) {
2741: /* f2 = t*f1+f0 */
2742: _bshiftup2(f1,-1,f2);
2743: _addup2_destructive(f2,f0);
2744: /* cyclic change of the variables */
2745: t = f0; f0 = f1; f1 = f2; f2 = t;
2746: }
2747: _copyup2(f1,f);
2748: return 0;
2749: break;
2750: default:
2751: return -1;
2752: break;
2753: }
2754: }
2755:
2756: /*
2757: * f = an irreducible trinomial or pentanomial of degree d 'after' f
2758: * return value : 0 --- exists
2759: * 1 --- does not exist (exhaustion)
2760: */
2761:
2762: int _generate_irreducible_polynomial(UP2 f,int d)
2763: {
2764: int ret,i,j,k,nz,i0,j0,k0;
2765: int w;
2766: unsigned int *fd;
2767:
2768: /*
2769: * if f = x^d+x^i+1 then i0 <- i, j0 <- 0, k0 <-0.
2770: * if f = x^d+x^k+x^j+x^i+1 (k>j>i) then i0 <- i, j0 <- j, k0 <-k.
2771: * otherwise i0,j0,k0 is set to 0.
2772: */
2773:
2774: fd = f->b;
2775: w = (d>>5)+1;
2776: if ( f->w && (d==degup2(f)) ) {
2777: for ( nz = 0, i = d; i >= 0; i-- )
2778: if ( fd[i>>5]&(1<<(i&31)) ) nz++;
2779: switch ( nz ) {
2780: case 3:
2781: for ( i0 = 1; !(fd[i0>>5]&(1<<(i0&31))) ; i0++ );
2782: /* reset i0-th bit */
2783: fd[i0>>5] &= ~(1<<(i0&31));
2784: j0 = k0 = 0;
2785: break;
2786: case 5:
2787: for ( i0 = 1; !(fd[i0>>5]&(1<<(i0&31))) ; i0++ );
2788: /* reset i0-th bit */
2789: fd[i0>>5] &= ~(1<<(i0&31));
2790: for ( j0 = i0+1; !(fd[j0>>5]&(1<<(j0&31))) ; j0++ );
2791: /* reset j0-th bit */
2792: fd[j0>>5] &= ~(1<<(j0&31));
2793: for ( k0 = j0+1; !(fd[k0>>5]&(1<<(k0&31))) ; k0++ );
2794: /* reset k0-th bit */
2795: fd[k0>>5] &= ~(1<<(k0&31));
2796: break;
2797: default:
2798: f->w = 0; break;
2799: }
2800: } else
2801: f->w = 0;
2802:
2803: if ( !f->w ) {
2804: fd = f->b;
2805: f->w = w; fd[0] |= 1; fd[d>>5] |= (1<<(d&31));
2806: i0 = j0 = k0 = 0;
2807: }
2808: /* if j0 > 0 then f is already a pentanomial */
2809: if ( j0 > 0 ) goto PENTA;
2810:
2811: /* searching for an irreducible trinomial */
2812:
2813: for ( i = 1; 2*i <= d; i++ ) {
2814: /* skip the polynomials 'before' f */
2815: if ( i < i0 ) continue;
2816: if ( i == i0 ) { i0 = 0; continue; }
2817: /* set i-th bit */
2818: fd[i>>5] |= (1<<(i&31));
2819: ret = irredcheck_dddup2(f);
2820: if ( ret == 1 ) return 0;
2821: /* reset i-th bit */
2822: fd[i>>5] &= ~(1<<(i&31));
2823: }
2824:
2825: /* searching for an irreducible pentanomial */
2826: PENTA:
2827: for ( i = 1; i < d; i++ ) {
2828: /* skip the polynomials 'before' f */
2829: if ( i < i0 ) continue;
2830: if ( i == i0 ) i0 = 0;
2831: /* set i-th bit */
2832: fd[i>>5] |= (1<<(i&31));
2833: for ( j = i+1; j < d; j++ ) {
2834: /* skip the polynomials 'before' f */
2835: if ( j < j0 ) continue;
2836: if ( j == j0 ) j0 = 0;
2837: /* set j-th bit */
2838: fd[j>>5] |= (1<<(j&31));
2839: for ( k = j+1; k < d; k++ ) {
2840: /* skip the polynomials 'before' f */
2841: if ( k < k0 ) continue;
2842: else if ( k == k0 ) { k0 = 0; continue; }
2843: /* set k-th bit */
2844: fd[k>>5] |= (1<<(k&31));
2845: ret = irredcheck_dddup2(f);
2846: if ( ret == 1 ) return 0;
2847: /* reset k-th bit */
2848: fd[k>>5] &= ~(1<<(k&31));
2849: }
2850: /* reset j-th bit */
2851: fd[j>>5] &= ~(1<<(j&31));
2852: }
2853: /* reset i-th bit */
2854: fd[i>>5] &= ~(1<<(i&31));
2855: }
2856: /* exhausted */
2857: return 1;
2858: }
2859:
2860: /*
2861: * f = an irreducible trinomial or pentanomial of degree d 'after' f
2862: *
2863: * searching strategy:
2864: * trinomial x^d+x^i+1:
2865: * i is as small as possible.
2866: * trinomial x^d+x^i+x^j+x^k+1:
2867: * i is as small as possible.
2868: * For such i, j is as small as possible.
2869: * For such i and j, 'k' is as small as possible.
2870: *
2871: * return value : 0 --- exists
2872: * 1 --- does not exist (exhaustion)
2873: */
2874:
2875: int _generate_good_irreducible_polynomial(UP2 f,int d)
2876: {
2877: int ret,i,j,k,nz,i0,j0,k0;
2878: int w;
2879: unsigned int *fd;
2880:
2881: /*
2882: * if f = x^d+x^i+1 then i0 <- i, j0 <- 0, k0 <-0.
2883: * if f = x^d+x^k+x^j+x^i+1 (k>j>i) then i0 <- i, j0 <- j, k0 <-k.
2884: * otherwise i0,j0,k0 is set to 0.
2885: */
2886:
2887: fd = f->b;
2888: w = (d>>5)+1;
2889: if ( f->w && (d==degup2(f)) ) {
2890: for ( nz = 0, i = d; i >= 0; i-- )
2891: if ( fd[i>>5]&(1<<(i&31)) ) nz++;
2892: switch ( nz ) {
2893: case 3:
2894: for ( i0 = 1; !(fd[i0>>5]&(1<<(i0&31))) ; i0++ );
2895: /* reset i0-th bit */
2896: fd[i0>>5] &= ~(1<<(i0&31));
2897: j0 = k0 = 0;
2898: break;
2899: case 5:
2900: for ( i0 = 1; !(fd[i0>>5]&(1<<(i0&31))) ; i0++ );
2901: /* reset i0-th bit */
2902: fd[i0>>5] &= ~(1<<(i0&31));
2903: for ( j0 = i0+1; !(fd[j0>>5]&(1<<(j0&31))) ; j0++ );
2904: /* reset j0-th bit */
2905: fd[j0>>5] &= ~(1<<(j0&31));
2906: for ( k0 = j0+1; !(fd[k0>>5]&(1<<(k0&31))) ; k0++ );
2907: /* reset k0-th bit */
2908: fd[k0>>5] &= ~(1<<(k0&31));
2909: break;
2910: default:
2911: f->w = 0; break;
2912: }
2913: } else
2914: f->w = 0;
2915:
2916: if ( !f->w ) {
2917: fd = f->b;
2918: f->w = w; fd[0] |= 1; fd[d>>5] |= (1<<(d&31));
2919: i0 = j0 = k0 = 0;
2920: }
2921: /* if j0 > 0 then f is already a pentanomial */
2922: if ( j0 > 0 ) goto PENTA;
2923:
2924: /* searching for an irreducible trinomial */
2925:
2926: for ( i = 1; 2*i <= d; i++ ) {
2927: /* skip the polynomials 'before' f */
2928: if ( i < i0 ) continue;
2929: if ( i == i0 ) { i0 = 0; continue; }
2930: /* set i-th bit */
2931: fd[i>>5] |= (1<<(i&31));
2932: ret = irredcheck_dddup2(f);
2933: if ( ret == 1 ) return 0;
2934: /* reset i-th bit */
2935: fd[i>>5] &= ~(1<<(i&31));
2936: }
2937:
2938: /* searching for an irreducible pentanomial */
2939: PENTA:
2940: for ( i = 3; i < d; i++ ) {
2941: /* skip the polynomials 'before' f */
2942: if ( i < i0 ) continue;
2943: if ( i == i0 ) i0 = 0;
2944: /* set i-th bit */
2945: fd[i>>5] |= (1<<(i&31));
2946: for ( j = 2; j < i; j++ ) {
2947: /* skip the polynomials 'before' f */
2948: if ( j < j0 ) continue;
2949: if ( j == j0 ) j0 = 0;
2950: /* set j-th bit */
2951: fd[j>>5] |= (1<<(j&31));
2952: for ( k = 1; k < j; k++ ) {
2953: /* skip the polynomials 'before' f */
2954: if ( k < k0 ) continue;
2955: else if ( k == k0 ) { k0 = 0; continue; }
2956: /* set k-th bit */
2957: fd[k>>5] |= (1<<(k&31));
2958: ret = irredcheck_dddup2(f);
2959: if ( ret == 1 ) return 0;
2960: /* reset k-th bit */
2961: fd[k>>5] &= ~(1<<(k&31));
2962: }
2963: /* reset j-th bit */
2964: fd[j>>5] &= ~(1<<(j&31));
2965: }
2966: /* reset i-th bit */
2967: fd[i>>5] &= ~(1<<(i&31));
2968: }
2969: /* exhausted */
2970: return 1;
2971: }
2972:
2973: void printqmat(Q **mat,int row,int col)
2974: {
2975: int i,j;
2976:
2977: for ( i = 0; i < row; i++ ) {
2978: for ( j = 0; j < col; j++ ) {
2979: printnum((Num)mat[i][j]); printf(" ");
2980: }
2981: printf("\n");
2982: }
2983: }
2984:
2985: void printimat(int **mat,int row,int col)
2986: {
2987: int i,j;
2988:
2989: for ( i = 0; i < row; i++ ) {
2990: for ( j = 0; j < col; j++ ) {
2991: printf("%d ",mat[i][j]);
2992: }
2993: printf("\n");
2994: }
2995: }
2996:
2997: void Pnd_det(NODE arg,P *rp)
2998: {
2999: if ( argc(arg) == 1 )
3000: nd_det(0,ARG0(arg),rp);
3001: else
1.2 ! noro 3002: nd_det(ZTOS((Q)ARG1(arg)),ARG0(arg),rp);
1.1 noro 3003: }
3004:
3005: void Pmat_col(NODE arg,VECT *rp)
3006: {
3007: int i,j,n;
3008: MAT mat;
3009: VECT vect;
3010:
3011: asir_assert(ARG0(arg),O_MAT,"mat_col");
3012: asir_assert(ARG1(arg),O_N,"mat_col");
3013: mat = (MAT)ARG0(arg);
1.2 ! noro 3014: j = ZTOS((Q)ARG1(arg));
1.1 noro 3015: if ( j < 0 || j >= mat->col) {
3016: error("mat_col : Out of range");
3017: }
3018: n = mat->row;
3019: MKVECT(vect,n);
3020: for(i=0; i<n; i++) {
3021: BDY(vect)[i] = BDY(mat)[i][j];
3022: }
3023: *rp = vect;
3024: }
3025:
3026: NODE triangleq(NODE e)
3027: {
3028: int n,i,k;
3029: V v;
3030: VL vl;
3031: P *p;
3032: NODE r,r1;
3033:
3034: n = length(e);
3035: p = (P *)MALLOC(n*sizeof(P));
3036: for ( i = 0; i < n; i++, e = NEXT(e) ) p[i] = (P)BDY(e);
3037: i = 0;
3038: while ( 1 ) {
3039: for ( ; i < n && !p[i]; i++ );
3040: if ( i == n ) break;
3041: if ( OID(p[i]) == O_N ) return 0;
3042: v = p[i]->v;
3043: for ( k = i+1; k < n; k++ )
3044: if ( p[k] ) {
3045: if ( OID(p[k]) == O_N ) return 0;
3046: if ( p[k]->v == v ) p[k] = 0;
3047: }
3048: i++;
3049: }
3050: for ( r = 0, i = 0; i < n; i++ ) {
3051: if ( p[i] ) {
3052: MKNODE(r1,p[i],r); r = r1;
3053: }
3054: }
3055: return r;
3056: }
3057:
3058: void Ptriangleq(NODE arg,LIST *rp)
3059: {
3060: NODE ret;
3061:
3062: asir_assert(ARG0(arg),O_LIST,"sparseleq");
3063: ret = triangleq(BDY((LIST)ARG0(arg)));
3064: MKLIST(*rp,ret);
3065: }
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 / asir2018 / builtin