Annotation of OpenXM_contrib2/asir2018/builtin/array.c, Revision 1.3
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.3 ! noro 48: * $OpenXM: OpenXM_contrib2/asir2018/builtin/array.c,v 1.2 2018/09/28 08:20:27 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
1.3 ! noro 1167: rank = generic_gauss_elim64(m,&nm,&dn,&ri,&ci);
1.1 noro 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:
1.3 ! noro 1232: #if SIZEOF_LONG==8
! 1233: void Pgeneric_gauss_elim_mod64(NODE arg,LIST *rp)
! 1234: {
! 1235: NODE n0;
! 1236: MAT m,mat;
! 1237: VECT rind,cind,rnum;
! 1238: Z **tmat;
! 1239: mp_limb_t **wmat,**row0;
! 1240: Z *rib,*cib,*rnb;
! 1241: int *colstat;
! 1242: mp_limb_t *p;
! 1243: Z q;
! 1244: mp_limb_t md;
! 1245: int i,j,k,l,row,col,t,rank;
! 1246:
! 1247: asir_assert(ARG0(arg),O_MAT,"generic_gauss_elim_mod64");
! 1248: asir_assert(ARG1(arg),O_N,"generic_gauss_elim_mod64");
! 1249: m = (MAT)ARG0(arg); md = ZTOS((Z)ARG1(arg));
! 1250: row = m->row; col = m->col; tmat = (Z **)m->body;
! 1251: wmat = (mp_limb_t **)almat64(row,col);
! 1252:
! 1253: row0 = (mp_limb_t **)ALLOCA(row*sizeof(mp_limb_t *));
! 1254: for ( i = 0; i < row; i++ ) row0[i] = wmat[i];
! 1255:
! 1256: colstat = (int *)MALLOC_ATOMIC(col*sizeof(int));
! 1257: for ( i = 0; i < row; i++ )
! 1258: for ( j = 0; j < col; j++ )
! 1259: wmat[i][j] = remqi64((Q)tmat[i][j],md);
! 1260: rank = generic_gauss_elim_mod64(wmat,row,col,md,colstat);
! 1261:
! 1262: MKVECT(rnum,rank);
! 1263: rnb = (Z *)rnum->body;
! 1264: for ( i = 0; i < rank; i++ )
! 1265: for ( j = 0, p = wmat[i]; j < row; j++ )
! 1266: if ( p == row0[j] )
! 1267: STOZ(j,rnb[i]);
! 1268:
! 1269: MKMAT(mat,rank,col-rank);
! 1270: tmat = (Z **)mat->body;
! 1271: for ( i = 0; i < rank; i++ )
! 1272: for ( j = k = 0; j < col; j++ )
! 1273: if ( !colstat[j] ) {
! 1274: UTOZ(wmat[i][j],tmat[i][k]); k++;
! 1275: }
! 1276:
! 1277: MKVECT(rind,rank);
! 1278: MKVECT(cind,col-rank);
! 1279: rib = (Z *)rind->body; cib = (Z *)cind->body;
! 1280: for ( j = k = l = 0; j < col; j++ )
! 1281: if ( colstat[j] ) {
! 1282: STOZ(j,rib[k]); k++;
! 1283: } else {
! 1284: STOZ(j,cib[l]); l++;
! 1285: }
! 1286: n0 = mknode(4,mat,rind,cind,rnum);
! 1287: MKLIST(*rp,n0);
! 1288: }
! 1289: #endif
! 1290:
1.1 noro 1291: void Pgeneric_gauss_elim_mod(NODE arg,LIST *rp)
1292: {
1293: NODE n0;
1294: MAT m,mat;
1295: VECT rind,cind,rnum;
1296: Z **tmat;
1297: int **wmat,**row0;
1298: Z *rib,*cib,*rnb;
1299: int *colstat,*p;
1300: Z q;
1.3 ! noro 1301: long mdl;
1.1 noro 1302: int md,i,j,k,l,row,col,t,rank;
1303:
1304: asir_assert(ARG0(arg),O_MAT,"generic_gauss_elim_mod");
1305: asir_assert(ARG1(arg),O_N,"generic_gauss_elim_mod");
1.3 ! noro 1306: #if SIZEOF_LONG==8
! 1307: mdl = ZTOS((Z)ARG1(arg));
! 1308: if ( mdl >= ((mp_limb_t)1)<<32 ) {
! 1309: Pgeneric_gauss_elim_mod64(arg,rp);
! 1310: return;
! 1311: }
! 1312: #endif
! 1313: m = (MAT)ARG0(arg);
! 1314: md = ZTOS((Z)ARG1(arg));
1.1 noro 1315: row = m->row; col = m->col; tmat = (Z **)m->body;
1316: wmat = (int **)almat(row,col);
1317:
1318: row0 = (int **)ALLOCA(row*sizeof(int *));
1319: for ( i = 0; i < row; i++ ) row0[i] = wmat[i];
1320:
1321: colstat = (int *)MALLOC_ATOMIC(col*sizeof(int));
1322: for ( i = 0; i < row; i++ )
1323: for ( j = 0; j < col; j++ )
1324: wmat[i][j] = remqi((Q)tmat[i][j],md);
1325: rank = generic_gauss_elim_mod(wmat,row,col,md,colstat);
1326:
1327: MKVECT(rnum,rank);
1328: rnb = (Z *)rnum->body;
1329: for ( i = 0; i < rank; i++ )
1330: for ( j = 0, p = wmat[i]; j < row; j++ )
1331: if ( p == row0[j] )
1.2 noro 1332: STOZ(j,rnb[i]);
1.1 noro 1333:
1334: MKMAT(mat,rank,col-rank);
1335: tmat = (Z **)mat->body;
1336: for ( i = 0; i < rank; i++ )
1337: for ( j = k = 0; j < col; j++ )
1338: if ( !colstat[j] ) {
1.2 noro 1339: UTOZ(wmat[i][j],tmat[i][k]); k++;
1.1 noro 1340: }
1341:
1342: MKVECT(rind,rank);
1343: MKVECT(cind,col-rank);
1344: rib = (Z *)rind->body; cib = (Z *)cind->body;
1345: for ( j = k = l = 0; j < col; j++ )
1346: if ( colstat[j] ) {
1.2 noro 1347: STOZ(j,rib[k]); k++;
1.1 noro 1348: } else {
1.2 noro 1349: STOZ(j,cib[l]); l++;
1.1 noro 1350: }
1351: n0 = mknode(4,mat,rind,cind,rnum);
1352: MKLIST(*rp,n0);
1353: }
1354:
1.3 ! noro 1355:
1.1 noro 1356: void Pleqm(NODE arg,VECT *rp)
1357: {
1358: MAT m;
1359: VECT vect;
1360: pointer **mat;
1361: Z *v;
1362: Z q;
1363: int **wmat;
1364: int md,i,j,row,col,t,n,status;
1365:
1366: asir_assert(ARG0(arg),O_MAT,"leqm");
1367: asir_assert(ARG1(arg),O_N,"leqm");
1.2 noro 1368: m = (MAT)ARG0(arg); md = ZTOS((Z)ARG1(arg));
1.1 noro 1369: row = m->row; col = m->col; mat = m->body;
1370: wmat = (int **)almat(row,col);
1371: for ( i = 0; i < row; i++ )
1372: for ( j = 0; j < col; j++ )
1373: wmat[i][j] = remqi((Q)mat[i][j],md);
1374: status = gauss_elim_mod(wmat,row,col,md);
1375: if ( status < 0 )
1376: *rp = 0;
1377: else if ( status > 0 )
1378: *rp = (VECT)ONE;
1379: else {
1380: n = col - 1;
1381: MKVECT(vect,n);
1382: for ( i = 0, v = (Z *)vect->body; i < n; i++ ) {
1.2 noro 1383: t = (md-wmat[i][n])%md; STOZ(t,v[i]);
1.1 noro 1384: }
1385: *rp = vect;
1386: }
1387: }
1388:
1389: int gauss_elim_mod(int **mat,int row,int col,int md)
1390: {
1391: int i,j,k,inv,a,n;
1392: int *t,*pivot;
1393:
1394: n = col - 1;
1395: for ( j = 0; j < n; j++ ) {
1396: for ( i = j; i < row && !mat[i][j]; i++ );
1397: if ( i == row )
1398: return 1;
1399: if ( i != j ) {
1400: t = mat[i]; mat[i] = mat[j]; mat[j] = t;
1401: }
1402: pivot = mat[j];
1403: inv = invm(pivot[j],md);
1404: for ( k = j; k <= n; k++ ) {
1405: /* pivot[k] = dmar(pivot[k],inv,0,md); */
1406: DMAR(pivot[k],inv,0,md,pivot[k])
1407: }
1408: for ( i = 0; i < row; i++ ) {
1409: t = mat[i];
1410: if ( i != j && (a = t[j]) )
1411: for ( k = j, a = md - a; k <= n; k++ ) {
1412: unsigned int tk;
1413: /* t[k] = dmar(pivot[k],a,t[k],md); */
1414: DMAR(pivot[k],a,t[k],md,tk)
1415: t[k] = tk;
1416: }
1417: }
1418: }
1419: for ( i = n; i < row && !mat[i][n]; i++ );
1420: if ( i == row )
1421: return 0;
1422: else
1423: return -1;
1424: }
1425:
1426: struct oEGT eg_mod,eg_elim,eg_elim1,eg_elim2,eg_chrem,eg_gschk,eg_intrat,eg_symb;
1427: struct oEGT eg_conv;
1428:
1429: #if 0
1430: void lu_dec_cr(MAT mat,MAT lu,Q *dn,int **perm);
1431:
1432: /* XXX broken */
1433: void lu_dec_cr(MAT mat,MAT lu,Q *dn,int **perm)
1434: {
1435: Q **a0,**b;
1436: Q *aiq;
1437: N **a;
1438: N *ai;
1439: Q q,q1,dn2,a1,q0,bik;
1440: MAT m;
1441: unsigned int md;
1442: int n,ind,i,j,rank,t,inv,t1,ret,min,k;
1443: int **w;
1444: int *wi,*rinfo0,*rinfo;
1445: N m1,m2,m3,u,s;
1446:
1447: a0 = (Q **)mat->body;
1448: n = mat->row;
1449: if ( n != mat->col )
1450: error("lu_dec_cr : non-square matrix");
1451: w = (int **)almat(n,n);
1452: MKMAT(m,n,n);
1453: a = (N **)m->body;
1454: UTON(1,m1);
1455: rinfo0 = 0;
1456: ind = 0;
1457: while ( 1 ) {
1458: md = get_lprime(ind);
1459: /* mat mod md */
1460: for ( i = 0; i < n; i++ )
1461: for ( j = 0, aiq = a0[i], wi = w[i]; j < n; j++ )
1462: if ( q = aiq[j] ) {
1463: t = rem(NM(q),md);
1464: if ( t && SGN(q) < 0 )
1465: t = (md - t) % md;
1466: wi[j] = t;
1467: } else
1468: wi[j] = 0;
1469:
1470: if ( !lu_mod((unsigned int **)w,n,md,&rinfo) ) continue;
1471: printf("."); fflush(stdout);
1472: if ( !rinfo0 )
1473: *perm = rinfo0 = rinfo;
1474: else {
1475: for ( i = 0; i < n; i++ )
1476: if ( rinfo[i] != rinfo0[i] ) break;
1477: if ( i < n ) continue;
1478: }
1479: if ( UNIN(m1) ) {
1480: for ( i = 0; i < n; i++ )
1481: for ( j = 0, ai = a[i], wi = w[i]; j < n; j++ ) {
1482: UTON(wi[j],u); ai[j] = u;
1483: }
1484: UTON(md,m1);
1485: } else {
1486: inv = invm(rem(m1,md),md);
1487: UTON(md,m2); muln(m1,m2,&m3);
1488: for ( i = 0; i < n; i++ )
1489: for ( j = 0, ai = a[i], wi = w[i]; j < n; j++ )
1490: if ( ai[i] ) {
1491: /* f3 = f1+m1*(m1 mod m2)^(-1)*(f2 - f1 mod m2) */
1492: t = rem(ai[j],md);
1493: if ( wi[j] >= t )
1494: t = wi[j]-t;
1495: else
1496: t = md-(t-wi[j]);
1497: DMAR(t,inv,0,md,t1)
1498: UTON(t1,u);
1499: muln(m1,u,&s);
1500: addn(ai[j],s,&u); ai[j] = u;
1501: } else if ( wi[j] ) {
1502: /* f3 = m1*(m1 mod m2)^(-1)*f2 */
1503: DMAR(wi[j],inv,0,md,t)
1504: UTON(t,u);
1505: muln(m1,u,&s); ai[j] = s;
1506: }
1507: m1 = m3;
1508: }
1509: if ( (++ind%8) == 0 ) {
1510: ret = intmtoratm(m,m1,lu,dn);
1511: if ( ret ) {
1512: b = (Q **)lu->body;
1513: mulq(*dn,*dn,&dn2);
1514: for ( i = 0; i < n; i++ ) {
1515: for ( j = 0; j < n; j++ ) {
1516: q = 0;
1517: min = MIN(i,j);
1518: for ( k = 0; k <= min; k++ ) {
1519: bik = k==i ? *dn : b[i][k];
1520: mulq(bik,b[k][j],&q0);
1521: addq(q,q0,&q1); q = q1;
1522: }
1523: mulq(a0[rinfo0[i]][j],dn2,&q1);
1524: if ( cmpq(q,q1) ) break;
1525: }
1526: if ( j < n ) break;
1527: }
1528: if ( i == n )
1529: return;
1530: }
1531: }
1532: }
1533: }
1534: #endif
1535:
1536:
1537: int f4_nocheck;
1538:
1539: #define ONE_STEP1 if ( zzz = *s ) { DMAR(zzz,hc,*tj,md,*tj) } tj++; s++;
1540:
1541: void reduce_reducers_mod(int **mat,int row,int col,int md)
1542: {
1543: int i,j,k,l,hc,zzz;
1544: int *t,*s,*tj,*ind;
1545:
1546: /* reduce the reducers */
1547: ind = (int *)ALLOCA(row*sizeof(int));
1548: for ( i = 0; i < row; i++ ) {
1549: t = mat[i];
1550: for ( j = 0; j < col && !t[j]; j++ );
1551: /* register the position of the head term */
1552: ind[i] = j;
1553: for ( l = i-1; l >= 0; l-- ) {
1554: /* reduce mat[i] by mat[l] */
1555: if ( hc = t[ind[l]] ) {
1556: /* mat[i] = mat[i]-hc*mat[l] */
1557: j = ind[l];
1558: s = mat[l]+j;
1559: tj = t+j;
1560: hc = md-hc;
1561: k = col-j;
1562: for ( ; k >= 64; k -= 64 ) {
1563: ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
1564: ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
1565: ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
1566: ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
1567: ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
1568: ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
1569: ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
1570: ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
1571: ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
1572: ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
1573: ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
1574: ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
1575: ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
1576: ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
1577: ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
1578: ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
1579: }
1580: for ( ; k > 0; k-- ) {
1581: if ( zzz = *s ) { DMAR(zzz,hc,*tj,md,*tj) } tj++; s++;
1582: }
1583: }
1584: }
1585: }
1586: }
1587:
1588: /*
1589: mat[i] : reducers (i=0,...,nred-1)
1590: spolys (i=nred,...,row-1)
1591: mat[0] < mat[1] < ... < mat[nred-1] w.r.t the term order
1592: 1. reduce the reducers
1593: 2. reduce spolys by the reduced reducers
1594: */
1595:
1596: void pre_reduce_mod(int **mat,int row,int col,int nred,int md)
1597: {
1598: int i,j,k,l,hc,inv;
1599: int *t,*s,*tk,*ind;
1600:
1601: #if 1
1602: /* reduce the reducers */
1603: ind = (int *)ALLOCA(row*sizeof(int));
1604: for ( i = 0; i < nred; i++ ) {
1605: /* make mat[i] monic and mat[i] by mat[0],...,mat[i-1] */
1606: t = mat[i];
1607: for ( j = 0; j < col && !t[j]; j++ );
1608: /* register the position of the head term */
1609: ind[i] = j;
1610: inv = invm(t[j],md);
1611: for ( k = j; k < col; k++ )
1612: if ( t[k] )
1613: DMAR(t[k],inv,0,md,t[k])
1614: for ( l = i-1; l >= 0; l-- ) {
1615: /* reduce mat[i] by mat[l] */
1616: if ( hc = t[ind[l]] ) {
1617: /* mat[i] = mat[i]-hc*mat[l] */
1618: for ( k = ind[l], hc = md-hc, s = mat[l]+k, tk = t+k;
1619: k < col; k++, tk++, s++ )
1620: if ( *s )
1621: DMAR(*s,hc,*tk,md,*tk)
1622: }
1623: }
1624: }
1625: /* reduce the spolys */
1626: for ( i = nred; i < row; i++ ) {
1627: t = mat[i];
1628: for ( l = nred-1; l >= 0; l-- ) {
1629: /* reduce mat[i] by mat[l] */
1630: if ( hc = t[ind[l]] ) {
1631: /* mat[i] = mat[i]-hc*mat[l] */
1632: for ( k = ind[l], hc = md-hc, s = mat[l]+k, tk = t+k;
1633: k < col; k++, tk++, s++ )
1634: if ( *s )
1635: DMAR(*s,hc,*tk,md,*tk)
1636: }
1637: }
1638: }
1639: #endif
1640: }
1641: /*
1642: mat[i] : reducers (i=0,...,nred-1)
1643: mat[0] < mat[1] < ... < mat[nred-1] w.r.t the term order
1644: */
1645:
1646: void reduce_sp_by_red_mod(int *sp,int **redmat,int *ind,int nred,int col,int md)
1647: {
1648: int i,j,k,hc,zzz;
1649: int *s,*tj;
1650:
1651: /* reduce the spolys by redmat */
1652: for ( i = nred-1; i >= 0; i-- ) {
1653: /* reduce sp by redmat[i] */
1654: if ( hc = sp[ind[i]] ) {
1655: /* sp = sp-hc*redmat[i] */
1656: j = ind[i];
1657: hc = md-hc;
1658: s = redmat[i]+j;
1659: tj = sp+j;
1660: for ( k = col-j; k > 0; k-- ) {
1661: if ( zzz = *s ) { DMAR(zzz,hc,*tj,md,*tj) } tj++; s++;
1662: }
1663: }
1664: }
1665: }
1666:
1667: /*
1668: mat[i] : compressed reducers (i=0,...,nred-1)
1669: mat[0] < mat[1] < ... < mat[nred-1] w.r.t the term order
1670: */
1671:
1672: void red_by_compress(int m,unsigned int *p,unsigned int *r,
1673: unsigned int *ri,unsigned int hc,int len)
1674: {
1675: unsigned int up,lo;
1676: unsigned int dmy;
1677: unsigned int *pj;
1678:
1679: p[*ri] = 0; r++; ri++;
1680: for ( len--; len; len--, r++, ri++ ) {
1681: pj = p+ *ri;
1682: DMA(*r,hc,*pj,up,lo);
1683: if ( up ) {
1684: DSAB(m,up,lo,dmy,*pj);
1685: } else
1686: *pj = lo;
1687: }
1688: }
1689:
1690: /* p -= hc*r */
1691:
1692: void red_by_vect(int m,unsigned int *p,unsigned int *r,unsigned int hc,int len)
1693: {
1694: unsigned int up,lo,dmy;
1695:
1696: *p++ = 0; r++; len--;
1697: for ( ; len; len--, r++, p++ )
1698: if ( *r ) {
1699: DMA(*r,hc,*p,up,lo);
1700: if ( up ) {
1701: DSAB(m,up,lo,dmy,*p);
1702: } else
1703: *p = lo;
1704: }
1705: }
1706:
1.3 ! noro 1707: #if SIZEOF_LONG==8
! 1708: /* mp_limb_t vector += U32 vector(normalized)*U32 */
1.1 noro 1709:
1.3 ! noro 1710: void red_by_vect64(int m, mp_limb_t *p,unsigned int *c,mp_limb_t *r,unsigned int hc,int len)
1.1 noro 1711: {
1.3 ! noro 1712: mp_limb_t t;
1.1 noro 1713:
1714: /* (p[0],c[0]) is normalized */
1715: *p++ = 0; *c++ = 0; r++; len--;
1716: for ( ; len; len--, r++, p++, c++ )
1717: if ( *r ) {
1718: t = (*p)+(*r)*hc;
1719: if ( t < *p ) (*c)++;
1720: *p = t;
1721: }
1722: }
1.3 ! noro 1723:
! 1724: /* mp_limb_t vector = (mp_limb_t vector+mp_limb_t vector*mp_limb_t)%mp_limb_t */
! 1725:
! 1726: void red_by_vect64mod(mp_limb_t m, mp_limb_t *p,mp_limb_t *r,mp_limb_t hc,int len)
! 1727: {
! 1728: *p++ = 0; r++; len--;
! 1729: for ( ; len; len--, r++, p++ )
! 1730: if ( *r )
! 1731: *p = muladdmod64(*r,hc,*p,m);
! 1732: }
! 1733:
! 1734: int generic_gauss_elim_mod64(mp_limb_t **mat,int row,int col,mp_limb_t md,int *colstat)
! 1735: {
! 1736: int i,j,k,l,rank;
! 1737: mp_limb_t inv,a;
! 1738: mp_limb_t *t,*pivot,*pk;
! 1739:
! 1740: for ( rank = 0, j = 0; j < col; j++ ) {
! 1741: for ( i = rank; i < row; i++ )
! 1742: if ( mat[i][j] )
! 1743: break;
! 1744: if ( i == row ) {
! 1745: colstat[j] = 0;
! 1746: continue;
! 1747: } else
! 1748: colstat[j] = 1;
! 1749: if ( i != rank ) {
! 1750: t = mat[i]; mat[i] = mat[rank]; mat[rank] = t;
! 1751: }
! 1752: pivot = mat[rank];
! 1753: inv = invmod64(pivot[j],md);
! 1754: for ( k = j, pk = pivot+k; k < col; k++, pk++ )
! 1755: if ( *pk )
! 1756: *pk = mulmod64(*pk,inv,md);
! 1757: for ( i = rank+1; i < row; i++ ) {
! 1758: t = mat[i];
! 1759: if ( a = t[j] )
! 1760: red_by_vect64mod(md,t+j,pivot+j,md-a,col-j);
! 1761: }
! 1762: rank++;
! 1763: }
! 1764: for ( j = col-1, l = rank-1; j >= 0; j-- )
! 1765: if ( colstat[j] ) {
! 1766: pivot = mat[l];
! 1767: for ( i = 0; i < l; i++ ) {
! 1768: t = mat[i];
! 1769: if ( a = t[j] )
! 1770: red_by_vect64mod(md,t+j,pivot+j,md-a,col-j);
! 1771: }
! 1772: l--;
! 1773: }
! 1774: return rank;
! 1775: }
! 1776:
1.1 noro 1777: #endif
1778:
1779: void red_by_vect_sf(int m,unsigned int *p,unsigned int *r,unsigned int hc,int len)
1780: {
1781: *p++ = 0; r++; len--;
1782: for ( ; len; len--, r++, p++ )
1783: if ( *r )
1784: *p = _addsf(_mulsf(*r,hc),*p);
1785: }
1786:
1787: extern Z current_mod_lf;
1788: extern int current_mod_lf_size;
1789:
1790: void red_by_vect_lf(mpz_t *p,mpz_t *r,mpz_t hc,int len)
1791: {
1792: mpz_set_ui(*p++,0); r++; len--;
1793: for ( ; len; len--, r++, p++ ) {
1794: mpz_addmul(*p,*r,hc);
1795: #if 0
1796: if ( mpz_size(*p) > current_mod_lf_size )
1797: mpz_mod(*p,*p,BDY(current_mod_lf));
1798: #endif
1799: }
1800: }
1801:
1802:
1803: extern unsigned int **psca;
1804:
1805: void reduce_sp_by_red_mod_compress (int *sp,CDP *redmat,int *ind,
1806: int nred,int col,int md)
1807: {
1808: int i,len;
1809: CDP ri;
1810: unsigned int hc;
1811: unsigned int *usp;
1812:
1813: usp = (unsigned int *)sp;
1814: /* reduce the spolys by redmat */
1815: for ( i = nred-1; i >= 0; i-- ) {
1816: /* reduce sp by redmat[i] */
1817: usp[ind[i]] %= md;
1818: if ( hc = usp[ind[i]] ) {
1819: /* sp = sp-hc*redmat[i] */
1820: hc = md-hc;
1821: ri = redmat[i];
1822: len = ri->len;
1823: red_by_compress(md,usp,psca[ri->psindex],ri->body,hc,len);
1824: }
1825: }
1826: for ( i = 0; i < col; i++ )
1827: if ( usp[i] >= (unsigned int)md )
1828: usp[i] %= md;
1829: }
1830:
1831: #define ONE_STEP2 if ( zzz = *pk ) { DMAR(zzz,a,*tk,md,*tk) } pk++; tk++;
1832:
1833: int generic_gauss_elim_mod(int **mat0,int row,int col,int md,int *colstat)
1834: {
1835: int i,j,k,l,inv,a,rank;
1836: unsigned int *t,*pivot,*pk;
1837: unsigned int **mat;
1838:
1839: mat = (unsigned int **)mat0;
1840: for ( rank = 0, j = 0; j < col; j++ ) {
1841: for ( i = rank; i < row; i++ )
1842: mat[i][j] %= md;
1843: for ( i = rank; i < row; i++ )
1844: if ( mat[i][j] )
1845: break;
1846: if ( i == row ) {
1847: colstat[j] = 0;
1848: continue;
1849: } else
1850: colstat[j] = 1;
1851: if ( i != rank ) {
1852: t = mat[i]; mat[i] = mat[rank]; mat[rank] = t;
1853: }
1854: pivot = mat[rank];
1855: inv = invm(pivot[j],md);
1856: for ( k = j, pk = pivot+k; k < col; k++, pk++ )
1857: if ( *pk ) {
1858: if ( *pk >= (unsigned int)md )
1859: *pk %= md;
1860: DMAR(*pk,inv,0,md,*pk)
1861: }
1862: for ( i = rank+1; i < row; i++ ) {
1863: t = mat[i];
1864: if ( a = t[j] )
1865: red_by_vect(md,t+j,pivot+j,md-a,col-j);
1866: }
1867: rank++;
1868: }
1869: for ( j = col-1, l = rank-1; j >= 0; j-- )
1870: if ( colstat[j] ) {
1871: pivot = mat[l];
1872: for ( i = 0; i < l; i++ ) {
1873: t = mat[i];
1874: t[j] %= md;
1875: if ( a = t[j] )
1876: red_by_vect(md,t+j,pivot+j,md-a,col-j);
1877: }
1878: l--;
1879: }
1880: for ( j = 0, l = 0; l < rank; j++ )
1881: if ( colstat[j] ) {
1882: t = mat[l];
1883: for ( k = j; k < col; k++ )
1884: if ( t[k] >= (unsigned int)md )
1885: t[k] %= md;
1886: l++;
1887: }
1888: return rank;
1889: }
1890:
1891: int generic_gauss_elim_mod2(int **mat0,int row,int col,int md,int *colstat,int *rowstat)
1892: {
1893: int i,j,k,l,inv,a,rank;
1894: unsigned int *t,*pivot,*pk;
1895: unsigned int **mat;
1896:
1897: for ( i = 0; i < row; i++ ) rowstat[i] = i;
1898: mat = (unsigned int **)mat0;
1899: for ( rank = 0, j = 0; j < col; j++ ) {
1900: for ( i = rank; i < row; i++ )
1901: mat[i][j] %= md;
1902: for ( i = rank; i < row; i++ )
1903: if ( mat[i][j] )
1904: break;
1905: if ( i == row ) {
1906: colstat[j] = 0;
1907: continue;
1908: } else
1909: colstat[j] = 1;
1910: if ( i != rank ) {
1911: t = mat[i]; mat[i] = mat[rank]; mat[rank] = t;
1912: k = rowstat[i]; rowstat[i] = rowstat[rank]; rowstat[rank] = k;
1913: }
1914: pivot = mat[rank];
1915: inv = invm(pivot[j],md);
1916: for ( k = j, pk = pivot+k; k < col; k++, pk++ )
1917: if ( *pk ) {
1918: if ( *pk >= (unsigned int)md )
1919: *pk %= md;
1920: DMAR(*pk,inv,0,md,*pk)
1921: }
1922: for ( i = rank+1; i < row; i++ ) {
1923: t = mat[i];
1924: if ( a = t[j] )
1925: red_by_vect(md,t+j,pivot+j,md-a,col-j);
1926: }
1927: rank++;
1928: }
1929: for ( j = col-1, l = rank-1; j >= 0; j-- )
1930: if ( colstat[j] ) {
1931: pivot = mat[l];
1932: for ( i = 0; i < l; i++ ) {
1933: t = mat[i];
1934: t[j] %= md;
1935: if ( a = t[j] )
1936: red_by_vect(md,t+j,pivot+j,md-a,col-j);
1937: }
1938: l--;
1939: }
1940: for ( j = 0, l = 0; l < rank; j++ )
1941: if ( colstat[j] ) {
1942: t = mat[l];
1943: for ( k = j; k < col; k++ )
1944: if ( t[k] >= (unsigned int)md )
1945: t[k] %= md;
1946: l++;
1947: }
1948: return rank;
1949: }
1950:
1951: int indep_rows_mod(int **mat0,int row,int col,int md,int *rowstat)
1952: {
1953: int i,j,k,l,inv,a,rank;
1954: unsigned int *t,*pivot,*pk;
1955: unsigned int **mat;
1956:
1957: for ( i = 0; i < row; i++ ) rowstat[i] = i;
1958: mat = (unsigned int **)mat0;
1959: for ( rank = 0, j = 0; j < col; j++ ) {
1960: for ( i = rank; i < row; i++ )
1961: mat[i][j] %= md;
1962: for ( i = rank; i < row; i++ )
1963: if ( mat[i][j] )
1964: break;
1965: if ( i == row ) continue;
1966: if ( i != rank ) {
1967: t = mat[i]; mat[i] = mat[rank]; mat[rank] = t;
1968: k = rowstat[i]; rowstat[i] = rowstat[rank]; rowstat[rank] = k;
1969: }
1970: pivot = mat[rank];
1971: inv = invm(pivot[j],md);
1972: for ( k = j, pk = pivot+k; k < col; k++, pk++ )
1973: if ( *pk ) {
1974: if ( *pk >= (unsigned int)md )
1975: *pk %= md;
1976: DMAR(*pk,inv,0,md,*pk)
1977: }
1978: for ( i = rank+1; i < row; i++ ) {
1979: t = mat[i];
1980: if ( a = t[j] )
1981: red_by_vect(md,t+j,pivot+j,md-a,col-j);
1982: }
1983: rank++;
1984: }
1985: return rank;
1986: }
1987:
1988: int generic_gauss_elim_sf(int **mat0,int row,int col,int md,int *colstat)
1989: {
1990: int i,j,k,l,inv,a,rank;
1991: unsigned int *t,*pivot,*pk;
1992: unsigned int **mat;
1993:
1994: mat = (unsigned int **)mat0;
1995: for ( rank = 0, j = 0; j < col; j++ ) {
1996: for ( i = rank; i < row; i++ )
1997: if ( mat[i][j] )
1998: break;
1999: if ( i == row ) {
2000: colstat[j] = 0;
2001: continue;
2002: } else
2003: colstat[j] = 1;
2004: if ( i != rank ) {
2005: t = mat[i]; mat[i] = mat[rank]; mat[rank] = t;
2006: }
2007: pivot = mat[rank];
2008: inv = _invsf(pivot[j]);
2009: for ( k = j, pk = pivot+k; k < col; k++, pk++ )
2010: if ( *pk )
2011: *pk = _mulsf(*pk,inv);
2012: for ( i = rank+1; i < row; i++ ) {
2013: t = mat[i];
2014: if ( a = t[j] )
2015: red_by_vect_sf(md,t+j,pivot+j,_chsgnsf(a),col-j);
2016: }
2017: rank++;
2018: }
2019: for ( j = col-1, l = rank-1; j >= 0; j-- )
2020: if ( colstat[j] ) {
2021: pivot = mat[l];
2022: for ( i = 0; i < l; i++ ) {
2023: t = mat[i];
2024: if ( a = t[j] )
2025: red_by_vect_sf(md,t+j,pivot+j,_chsgnsf(a),col-j);
2026: }
2027: l--;
2028: }
2029: return rank;
2030: }
2031:
2032: /* LU decomposition; a[i][i] = 1/U[i][i] */
2033:
2034: int lu_gfmmat(GFMMAT mat,unsigned int md,int *perm)
2035: {
2036: int row,col;
2037: int i,j,k;
2038: unsigned int *t,*pivot;
2039: unsigned int **a;
2040: unsigned int inv,m;
2041:
2042: row = mat->row; col = mat->col;
2043: a = mat->body;
2044: bzero(perm,row*sizeof(int));
2045:
2046: for ( i = 0; i < row; i++ )
2047: perm[i] = i;
2048: for ( k = 0; k < col; k++ ) {
2049: for ( i = k; i < row && !a[i][k]; i++ );
2050: if ( i == row )
2051: return 0;
2052: if ( i != k ) {
2053: j = perm[i]; perm[i] = perm[k]; perm[k] = j;
2054: t = a[i]; a[i] = a[k]; a[k] = t;
2055: }
2056: pivot = a[k];
2057: pivot[k] = inv = invm(pivot[k],md);
2058: for ( i = k+1; i < row; i++ ) {
2059: t = a[i];
2060: if ( m = t[k] ) {
2061: DMAR(inv,m,0,md,t[k])
2062: for ( j = k+1, m = md - t[k]; j < col; j++ )
2063: if ( pivot[j] ) {
2064: unsigned int tj;
2065:
2066: DMAR(m,pivot[j],t[j],md,tj)
2067: t[j] = tj;
2068: }
2069: }
2070: }
2071: }
2072: return 1;
2073: }
2074:
2075: /*
2076: Input
2077: a: a row x col matrix
2078: md : a modulus
2079:
2080: Output:
2081: return : d = the rank of mat
2082: a[0..(d-1)][0..(d-1)] : LU decomposition (a[i][i] = 1/U[i][i])
2083: rinfo: array of length row
2084: cinfo: array of length col
2085: i-th row in new a <-> rinfo[i]-th row in old a
2086: cinfo[j]=1 <=> j-th column is contained in the LU decomp.
2087: */
2088:
2089: int find_lhs_and_lu_mod(unsigned int **a,int row,int col,
2090: unsigned int md,int **rinfo,int **cinfo)
2091: {
2092: int i,j,k,d;
2093: int *rp,*cp;
2094: unsigned int *t,*pivot;
2095: unsigned int inv,m;
2096:
2097: *rinfo = rp = (int *)MALLOC_ATOMIC(row*sizeof(int));
2098: *cinfo = cp = (int *)MALLOC_ATOMIC(col*sizeof(int));
2099: for ( i = 0; i < row; i++ )
2100: rp[i] = i;
2101: for ( k = 0, d = 0; k < col; k++ ) {
2102: for ( i = d; i < row && !a[i][k]; i++ );
2103: if ( i == row ) {
2104: cp[k] = 0;
2105: continue;
2106: } else
2107: cp[k] = 1;
2108: if ( i != d ) {
2109: j = rp[i]; rp[i] = rp[d]; rp[d] = j;
2110: t = a[i]; a[i] = a[d]; a[d] = t;
2111: }
2112: pivot = a[d];
2113: pivot[k] = inv = invm(pivot[k],md);
2114: for ( i = d+1; i < row; i++ ) {
2115: t = a[i];
2116: if ( m = t[k] ) {
2117: DMAR(inv,m,0,md,t[k])
2118: for ( j = k+1, m = md - t[k]; j < col; j++ )
2119: if ( pivot[j] ) {
2120: unsigned int tj;
2121: DMAR(m,pivot[j],t[j],md,tj)
2122: t[j] = tj;
2123: }
2124: }
2125: }
2126: d++;
2127: }
2128: return d;
2129: }
2130:
2131: int lu_mod(unsigned int **a,int n,unsigned int md,int **rinfo)
2132: {
2133: int i,j,k;
2134: int *rp;
2135: unsigned int *t,*pivot;
2136: unsigned int inv,m;
2137:
2138: *rinfo = rp = (int *)MALLOC_ATOMIC(n*sizeof(int));
2139: for ( i = 0; i < n; i++ ) rp[i] = i;
2140: for ( k = 0; k < n; k++ ) {
2141: for ( i = k; i < n && !a[i][k]; i++ );
2142: if ( i == n ) return 0;
2143: if ( i != k ) {
2144: j = rp[i]; rp[i] = rp[k]; rp[k] = j;
2145: t = a[i]; a[i] = a[k]; a[k] = t;
2146: }
2147: pivot = a[k];
2148: inv = invm(pivot[k],md);
2149: for ( i = k+1; i < n; i++ ) {
2150: t = a[i];
2151: if ( m = t[k] ) {
2152: DMAR(inv,m,0,md,t[k])
2153: for ( j = k+1, m = md - t[k]; j < n; j++ )
2154: if ( pivot[j] ) {
2155: unsigned int tj;
2156: DMAR(m,pivot[j],t[j],md,tj)
2157: t[j] = tj;
2158: }
2159: }
2160: }
2161: }
2162: return 1;
2163: }
2164:
2165: /*
2166: Input
2167: a : n x n matrix; a result of LU-decomposition
2168: md : modulus
2169: b : n x l matrix
2170: Output
2171: b = a^(-1)b
2172: */
2173:
2174: void solve_by_lu_mod(int **a,int n,int md,int **b,int l,int normalize)
2175: {
2176: unsigned int *y,*c;
2177: int i,j,k;
2178: unsigned int t,m,m2;
2179:
2180: y = (int *)MALLOC_ATOMIC(n*sizeof(int));
2181: c = (int *)MALLOC_ATOMIC(n*sizeof(int));
2182: m2 = md>>1;
2183: for ( k = 0; k < l; k++ ) {
2184: /* copy b[.][k] to c */
2185: for ( i = 0; i < n; i++ )
2186: c[i] = (unsigned int)b[i][k];
2187: /* solve Ly=c */
2188: for ( i = 0; i < n; i++ ) {
2189: for ( t = c[i], j = 0; j < i; j++ )
2190: if ( a[i][j] ) {
2191: m = md - a[i][j];
2192: DMAR(m,y[j],t,md,t)
2193: }
2194: y[i] = t;
2195: }
2196: /* solve Uc=y */
2197: for ( i = n-1; i >= 0; i-- ) {
2198: for ( t = y[i], j =i+1; j < n; j++ )
2199: if ( a[i][j] ) {
2200: m = md - a[i][j];
2201: DMAR(m,c[j],t,md,t)
2202: }
2203: /* a[i][i] = 1/U[i][i] */
2204: DMAR(t,a[i][i],0,md,c[i])
2205: }
2206: /* copy c to b[.][k] with normalization */
2207: if ( normalize )
2208: for ( i = 0; i < n; i++ )
2209: b[i][k] = (int)(c[i]>m2 ? c[i]-md : c[i]);
2210: else
2211: for ( i = 0; i < n; i++ )
2212: b[i][k] = c[i];
2213: }
2214: }
2215:
2216: void Pleqm1(NODE arg,VECT *rp)
2217: {
2218: MAT m;
2219: VECT vect;
2220: pointer **mat;
2221: Z *v;
2222: Z q;
2223: int **wmat;
2224: int md,i,j,row,col,t,n,status;
2225:
2226: asir_assert(ARG0(arg),O_MAT,"leqm1");
2227: asir_assert(ARG1(arg),O_N,"leqm1");
1.2 noro 2228: m = (MAT)ARG0(arg); md = ZTOS((Z)ARG1(arg));
1.1 noro 2229: row = m->row; col = m->col; mat = m->body;
2230: wmat = (int **)almat(row,col);
2231: for ( i = 0; i < row; i++ )
2232: for ( j = 0; j < col; j++ )
2233: wmat[i][j] = remqi((Q)mat[i][j],md);
2234: status = gauss_elim_mod1(wmat,row,col,md);
2235: if ( status < 0 )
2236: *rp = 0;
2237: else if ( status > 0 )
2238: *rp = (VECT)ONE;
2239: else {
2240: n = col - 1;
2241: MKVECT(vect,n);
2242: for ( i = 0, v = (Z *)vect->body; i < n; i++ ) {
1.2 noro 2243: t = (md-wmat[i][n])%md; STOZ(t,v[i]);
1.1 noro 2244: }
2245: *rp = vect;
2246: }
2247: }
2248:
2249: int gauss_elim_mod1(int **mat,int row,int col,int md)
2250: {
2251: int i,j,k,inv,a,n;
2252: int *t,*pivot;
2253:
2254: n = col - 1;
2255: for ( j = 0; j < n; j++ ) {
2256: for ( i = j; i < row && !mat[i][j]; i++ );
2257: if ( i == row )
2258: return 1;
2259: if ( i != j ) {
2260: t = mat[i]; mat[i] = mat[j]; mat[j] = t;
2261: }
2262: pivot = mat[j];
2263: inv = invm(pivot[j],md);
2264: for ( k = j; k <= n; k++ )
2265: pivot[k] = dmar(pivot[k],inv,0,md);
2266: for ( i = j+1; i < row; i++ ) {
2267: t = mat[i];
2268: if ( i != j && (a = t[j]) )
2269: for ( k = j, a = md - a; k <= n; k++ )
2270: t[k] = dmar(pivot[k],a,t[k],md);
2271: }
2272: }
2273: for ( i = n; i < row && !mat[i][n]; i++ );
2274: if ( i == row ) {
2275: for ( j = n-1; j >= 0; j-- ) {
2276: for ( i = j-1, a = (md-mat[j][n])%md; i >= 0; i-- ) {
2277: mat[i][n] = dmar(mat[i][j],a,mat[i][n],md);
2278: mat[i][j] = 0;
2279: }
2280: }
2281: return 0;
2282: } else
2283: return -1;
2284: }
2285:
2286: void Pgeninvm(NODE arg,LIST *rp)
2287: {
2288: MAT m;
2289: pointer **mat;
2290: Z **tmat;
2291: Z q;
2292: unsigned int **wmat;
2293: int md,i,j,row,col,t,status;
2294: MAT mat1,mat2;
2295: NODE node1,node2;
2296:
2297: asir_assert(ARG0(arg),O_MAT,"leqm1");
2298: asir_assert(ARG1(arg),O_N,"leqm1");
1.2 noro 2299: m = (MAT)ARG0(arg); md = ZTOS((Q)ARG1(arg));
1.1 noro 2300: row = m->row; col = m->col; mat = m->body;
2301: wmat = (unsigned int **)almat(row,col+row);
2302: for ( i = 0; i < row; i++ ) {
2303: bzero((char *)wmat[i],(col+row)*sizeof(int));
2304: for ( j = 0; j < col; j++ )
2305: wmat[i][j] = remqi((Q)mat[i][j],md);
2306: }
2307: status = gauss_elim_geninv_mod(wmat,row,col,md);
2308: if ( status > 0 )
2309: *rp = 0;
2310: else {
2311: MKMAT(mat1,col,row); MKMAT(mat2,row-col,row);
2312: for ( i = 0, tmat = (Z **)mat1->body; i < col; i++ )
2313: for ( j = 0; j < row; j++ )
1.2 noro 2314: UTOZ(wmat[i][j+col],tmat[i][j]);
1.1 noro 2315: for ( tmat = (Z **)mat2->body; i < row; i++ )
2316: for ( j = 0; j < row; j++ )
1.2 noro 2317: UTOZ(wmat[i][j+col],tmat[i-col][j]);
1.1 noro 2318: MKNODE(node2,mat2,0); MKNODE(node1,mat1,node2); MKLIST(*rp,node1);
2319: }
2320: }
2321:
2322: int gauss_elim_geninv_mod(unsigned int **mat,int row,int col,int md)
2323: {
2324: int i,j,k,inv,a,n,m;
2325: unsigned int *t,*pivot;
2326:
2327: n = col; m = row+col;
2328: for ( j = 0; j < n; j++ ) {
2329: for ( i = j; i < row && !mat[i][j]; i++ );
2330: if ( i == row )
2331: return 1;
2332: if ( i != j ) {
2333: t = mat[i]; mat[i] = mat[j]; mat[j] = t;
2334: }
2335: pivot = mat[j];
2336: inv = invm(pivot[j],md);
2337: for ( k = j; k < m; k++ )
2338: pivot[k] = dmar(pivot[k],inv,0,md);
2339: for ( i = j+1; i < row; i++ ) {
2340: t = mat[i];
2341: if ( a = t[j] )
2342: for ( k = j, a = md - a; k < m; k++ )
2343: t[k] = dmar(pivot[k],a,t[k],md);
2344: }
2345: }
2346: for ( j = n-1; j >= 0; j-- ) {
2347: pivot = mat[j];
2348: for ( i = j-1; i >= 0; i-- ) {
2349: t = mat[i];
2350: if ( a = t[j] )
2351: for ( k = j, a = md - a; k < m; k++ )
2352: t[k] = dmar(pivot[k],a,t[k],md);
2353: }
2354: }
2355: return 0;
2356: }
2357:
2358: void Psolve_by_lu_gfmmat(NODE arg,VECT *rp)
2359: {
2360: GFMMAT lu;
2361: Z *perm,*rhs,*v;
2362: int n,i;
2363: unsigned int md;
2364: unsigned int *b,*sol;
2365: VECT r;
2366:
2367: lu = (GFMMAT)ARG0(arg);
2368: perm = (Z *)BDY((VECT)ARG1(arg));
2369: rhs = (Z *)BDY((VECT)ARG2(arg));
1.2 noro 2370: md = (unsigned int)ZTOS((Z)ARG3(arg));
1.1 noro 2371: n = lu->col;
2372: b = (unsigned int *)MALLOC_ATOMIC(n*sizeof(int));
2373: sol = (unsigned int *)MALLOC_ATOMIC(n*sizeof(int));
2374: for ( i = 0; i < n; i++ )
1.2 noro 2375: b[i] = ZTOS(rhs[ZTOS(perm[i])]);
1.1 noro 2376: solve_by_lu_gfmmat(lu,md,b,sol);
2377: MKVECT(r,n);
2378: for ( i = 0, v = (Z *)r->body; i < n; i++ )
1.2 noro 2379: UTOZ(sol[i],v[i]);
1.1 noro 2380: *rp = r;
2381: }
2382:
2383: void solve_by_lu_gfmmat(GFMMAT lu,unsigned int md,
2384: unsigned int *b,unsigned int *x)
2385: {
2386: int n;
2387: unsigned int **a;
2388: unsigned int *y;
2389: int i,j;
2390: unsigned int t,m;
2391:
2392: n = lu->col;
2393: a = lu->body;
2394: y = (unsigned int *)MALLOC_ATOMIC(n*sizeof(int));
2395: /* solve Ly=b */
2396: for ( i = 0; i < n; i++ ) {
2397: for ( t = b[i], j = 0; j < i; j++ )
2398: if ( a[i][j] ) {
2399: m = md - a[i][j];
2400: DMAR(m,y[j],t,md,t)
2401: }
2402: y[i] = t;
2403: }
2404: /* solve Ux=y */
2405: for ( i = n-1; i >= 0; i-- ) {
2406: for ( t = y[i], j =i+1; j < n; j++ )
2407: if ( a[i][j] ) {
2408: m = md - a[i][j];
2409: DMAR(m,x[j],t,md,t)
2410: }
2411: /* a[i][i] = 1/U[i][i] */
2412: DMAR(t,a[i][i],0,md,x[i])
2413: }
2414: }
2415:
2416: #if 0
2417: void Plu_mat(NODE arg,LIST *rp)
2418: {
2419: MAT m,lu;
2420: Q dn;
2421: Q *v;
2422: int n,i;
2423: int *iperm;
2424: VECT perm;
2425: NODE n0;
2426:
2427: asir_assert(ARG0(arg),O_MAT,"lu_mat");
2428: m = (MAT)ARG0(arg);
2429: n = m->row;
2430: MKMAT(lu,n,n);
2431: lu_dec_cr(m,lu,&dn,&iperm);
2432: MKVECT(perm,n);
2433: for ( i = 0, v = (Q *)perm->body; i < n; i++ )
1.2 noro 2434: STOZ(iperm[i],v[i]);
1.1 noro 2435: n0 = mknode(3,lu,dn,perm);
2436: MKLIST(*rp,n0);
2437: }
2438: #endif
2439:
2440: void Plu_gfmmat(NODE arg,LIST *rp)
2441: {
2442: MAT m;
2443: GFMMAT mm;
2444: unsigned int md;
2445: int i,row,col,status;
2446: int *iperm;
2447: Z *v;
2448: VECT perm;
2449: NODE n0;
2450:
2451: asir_assert(ARG0(arg),O_MAT,"lu_gfmmat");
2452: asir_assert(ARG1(arg),O_N,"lu_gfmmat");
1.2 noro 2453: m = (MAT)ARG0(arg); md = (unsigned int)ZTOS((Q)ARG1(arg));
1.1 noro 2454: mat_to_gfmmat(m,md,&mm);
2455: row = m->row;
2456: col = m->col;
2457: iperm = (int *)MALLOC_ATOMIC(row*sizeof(int));
2458: status = lu_gfmmat(mm,md,iperm);
2459: if ( !status )
2460: n0 = 0;
2461: else {
2462: MKVECT(perm,row);
2463: for ( i = 0, v = (Z *)perm->body; i < row; i++ )
1.2 noro 2464: STOZ(iperm[i],v[i]);
1.1 noro 2465: n0 = mknode(2,mm,perm);
2466: }
2467: MKLIST(*rp,n0);
2468: }
2469:
2470: void Pmat_to_gfmmat(NODE arg,GFMMAT *rp)
2471: {
2472: MAT m;
2473: unsigned int md;
2474:
2475: asir_assert(ARG0(arg),O_MAT,"mat_to_gfmmat");
2476: asir_assert(ARG1(arg),O_N,"mat_to_gfmmat");
1.2 noro 2477: m = (MAT)ARG0(arg); md = (unsigned int)ZTOS((Q)ARG1(arg));
1.1 noro 2478: mat_to_gfmmat(m,md,rp);
2479: }
2480:
2481: void mat_to_gfmmat(MAT m,unsigned int md,GFMMAT *rp)
2482: {
2483: unsigned int **wmat;
2484: unsigned int t;
2485: Z **mat;
2486: Z q;
2487: int i,j,row,col;
2488:
2489: row = m->row; col = m->col; mat = (Z **)m->body;
2490: wmat = (unsigned int **)almat(row,col);
2491: for ( i = 0; i < row; i++ ) {
2492: bzero((char *)wmat[i],col*sizeof(unsigned int));
2493: for ( j = 0; j < col; j++ )
2494: wmat[i][j] = remqi((Q)mat[i][j],md);
2495: }
2496: TOGFMMAT(row,col,wmat,*rp);
2497: }
2498:
2499: void Pgeninvm_swap(NODE arg,LIST *rp)
2500: {
2501: MAT m;
2502: pointer **mat;
2503: Z **tmat;
2504: Z *tvect;
2505: Z q;
2506: unsigned int **wmat,**invmat;
2507: int *index;
2508: unsigned int t,md;
2509: int i,j,row,col,status;
2510: MAT mat1;
2511: VECT vect1;
2512: NODE node1,node2;
2513:
2514: asir_assert(ARG0(arg),O_MAT,"geninvm_swap");
2515: asir_assert(ARG1(arg),O_N,"geninvm_swap");
1.2 noro 2516: m = (MAT)ARG0(arg); md = ZTOS((Z)ARG1(arg));
1.1 noro 2517: row = m->row; col = m->col; mat = m->body;
2518: wmat = (unsigned int **)almat(row,col+row);
2519: for ( i = 0; i < row; i++ ) {
2520: bzero((char *)wmat[i],(col+row)*sizeof(int));
2521: for ( j = 0; j < col; j++ )
2522: wmat[i][j] = remqi((Q)mat[i][j],md);
2523: wmat[i][col+i] = 1;
2524: }
2525: status = gauss_elim_geninv_mod_swap(wmat,row,col,md,&invmat,&index);
2526: if ( status > 0 )
2527: *rp = 0;
2528: else {
2529: MKMAT(mat1,col,col);
2530: for ( i = 0, tmat = (Z **)mat1->body; i < col; i++ )
2531: for ( j = 0; j < col; j++ )
1.2 noro 2532: UTOZ(invmat[i][j],tmat[i][j]);
1.1 noro 2533: MKVECT(vect1,row);
2534: for ( i = 0, tvect = (Z *)vect1->body; i < row; i++ )
1.2 noro 2535: STOZ(index[i],tvect[i]);
1.1 noro 2536: MKNODE(node2,vect1,0); MKNODE(node1,mat1,node2); MKLIST(*rp,node1);
2537: }
2538: }
2539:
2540: int gauss_elim_geninv_mod_swap(unsigned int **mat,int row,int col,unsigned int md,
2541: unsigned int ***invmatp,int **indexp)
2542: {
2543: int i,j,k,inv,a,n,m;
2544: unsigned int *t,*pivot,*s;
2545: int *index;
2546: unsigned int **invmat;
2547:
2548: n = col; m = row+col;
2549: *indexp = index = (int *)MALLOC_ATOMIC(row*sizeof(int));
2550: for ( i = 0; i < row; i++ )
2551: index[i] = i;
2552: for ( j = 0; j < n; j++ ) {
2553: for ( i = j; i < row && !mat[i][j]; i++ );
2554: if ( i == row ) {
2555: *indexp = 0; *invmatp = 0; return 1;
2556: }
2557: if ( i != j ) {
2558: t = mat[i]; mat[i] = mat[j]; mat[j] = t;
2559: k = index[i]; index[i] = index[j]; index[j] = k;
2560: }
2561: pivot = mat[j];
2562: inv = (unsigned int)invm(pivot[j],md);
2563: for ( k = j; k < m; k++ )
2564: if ( pivot[k] )
2565: pivot[k] = (unsigned int)dmar(pivot[k],inv,0,md);
2566: for ( i = j+1; i < row; i++ ) {
2567: t = mat[i];
2568: if ( a = t[j] )
2569: for ( k = j, a = md - a; k < m; k++ )
2570: if ( pivot[k] )
2571: t[k] = dmar(pivot[k],a,t[k],md);
2572: }
2573: }
2574: for ( j = n-1; j >= 0; j-- ) {
2575: pivot = mat[j];
2576: for ( i = j-1; i >= 0; i-- ) {
2577: t = mat[i];
2578: if ( a = t[j] )
2579: for ( k = j, a = md - a; k < m; k++ )
2580: if ( pivot[k] )
2581: t[k] = dmar(pivot[k],a,t[k],md);
2582: }
2583: }
2584: *invmatp = invmat = (unsigned int **)almat(col,col);
2585: for ( i = 0; i < col; i++ )
2586: for ( j = 0, s = invmat[i], t = mat[i]; j < col; j++ )
2587: s[j] = t[col+index[j]];
2588: return 0;
2589: }
2590:
2591: int gauss_elim_geninv_sf_swap(int **mat,int row,int col, int ***invmatp,int **indexp);
2592:
2593: void Pgeninv_sf_swap(NODE arg,LIST *rp)
2594: {
2595: MAT m;
2596: GFS **mat,**tmat;
2597: Z *tvect;
2598: GFS q;
2599: int **wmat,**invmat;
2600: int *index;
2601: unsigned int t;
2602: int i,j,row,col,status;
2603: MAT mat1;
2604: VECT vect1;
2605: NODE node1,node2;
2606:
2607: asir_assert(ARG0(arg),O_MAT,"geninv_sf_swap");
2608: m = (MAT)ARG0(arg);
2609: row = m->row; col = m->col; mat = (GFS **)m->body;
2610: wmat = (int **)almat(row,col+row);
2611: for ( i = 0; i < row; i++ ) {
2612: bzero((char *)wmat[i],(col+row)*sizeof(int));
2613: for ( j = 0; j < col; j++ )
2614: if ( q = (GFS)mat[i][j] )
2615: wmat[i][j] = FTOIF(CONT(q));
2616: wmat[i][col+i] = _onesf();
2617: }
2618: status = gauss_elim_geninv_sf_swap(wmat,row,col,&invmat,&index);
2619: if ( status > 0 )
2620: *rp = 0;
2621: else {
2622: MKMAT(mat1,col,col);
2623: for ( i = 0, tmat = (GFS **)mat1->body; i < col; i++ )
2624: for ( j = 0; j < col; j++ )
2625: if ( t = invmat[i][j] ) {
2626: MKGFS(IFTOF(t),tmat[i][j]);
2627: }
2628: MKVECT(vect1,row);
2629: for ( i = 0, tvect = (Z *)vect1->body; i < row; i++ )
1.2 noro 2630: STOZ(index[i],tvect[i]);
1.1 noro 2631: MKNODE(node2,vect1,0); MKNODE(node1,mat1,node2); MKLIST(*rp,node1);
2632: }
2633: }
2634:
2635: int gauss_elim_geninv_sf_swap(int **mat,int row,int col, int ***invmatp,int **indexp)
2636: {
2637: int i,j,k,inv,a,n,m,u;
2638: int *t,*pivot,*s;
2639: int *index;
2640: int **invmat;
2641:
2642: n = col; m = row+col;
2643: *indexp = index = (int *)MALLOC_ATOMIC(row*sizeof(int));
2644: for ( i = 0; i < row; i++ )
2645: index[i] = i;
2646: for ( j = 0; j < n; j++ ) {
2647: for ( i = j; i < row && !mat[i][j]; i++ );
2648: if ( i == row ) {
2649: *indexp = 0; *invmatp = 0; return 1;
2650: }
2651: if ( i != j ) {
2652: t = mat[i]; mat[i] = mat[j]; mat[j] = t;
2653: k = index[i]; index[i] = index[j]; index[j] = k;
2654: }
2655: pivot = mat[j];
2656: inv = _invsf(pivot[j]);
2657: for ( k = j; k < m; k++ )
2658: if ( pivot[k] )
2659: pivot[k] = _mulsf(pivot[k],inv);
2660: for ( i = j+1; i < row; i++ ) {
2661: t = mat[i];
2662: if ( a = t[j] )
2663: for ( k = j, a = _chsgnsf(a); k < m; k++ )
2664: if ( pivot[k] ) {
2665: u = _mulsf(pivot[k],a);
2666: t[k] = _addsf(u,t[k]);
2667: }
2668: }
2669: }
2670: for ( j = n-1; j >= 0; j-- ) {
2671: pivot = mat[j];
2672: for ( i = j-1; i >= 0; i-- ) {
2673: t = mat[i];
2674: if ( a = t[j] )
2675: for ( k = j, a = _chsgnsf(a); k < m; k++ )
2676: if ( pivot[k] ) {
2677: u = _mulsf(pivot[k],a);
2678: t[k] = _addsf(u,t[k]);
2679: }
2680: }
2681: }
2682: *invmatp = invmat = (int **)almat(col,col);
2683: for ( i = 0; i < col; i++ )
2684: for ( j = 0, s = invmat[i], t = mat[i]; j < col; j++ )
2685: s[j] = t[col+index[j]];
2686: return 0;
2687: }
2688:
2689: void inner_product_int(Z *a,Z *b,int n,Z *r)
2690: {
2691: int i;
2692: Z t;
2693:
2694: t = 0;
2695: for ( i = 0; i < n; i++ )
2696: muladdtoz(a[i],b[i],&t);
2697: *r = t;
2698: }
2699:
2700: void Pmul_mat_vect_int(NODE arg,VECT *rp)
2701: {
2702: MAT mat;
2703: VECT vect,r;
2704: int row,col,i;
2705:
2706: mat = (MAT)ARG0(arg);
2707: vect = (VECT)ARG1(arg);
2708: row = mat->row;
2709: col = mat->col;
2710: MKVECT(r,row);
2711: for ( i = 0; i < row; i++ ) {
2712: inner_product_int((Z *)mat->body[i],(Z *)vect->body,col,(Z *)&r->body[i]);
2713: }
2714: *rp = r;
2715: }
2716:
2717: void Pnbpoly_up2(NODE arg,GF2N *rp)
2718: {
2719: int m,type,ret;
2720: UP2 r;
2721:
1.2 noro 2722: m = ZTOS((Z)ARG0(arg));
2723: type = ZTOS((Z)ARG1(arg));
1.1 noro 2724: ret = generate_ONB_polynomial(&r,m,type);
2725: if ( ret == 0 )
2726: MKGF2N(r,*rp);
2727: else
2728: *rp = 0;
2729: }
2730:
2731: void Px962_irredpoly_up2(NODE arg,GF2N *rp)
2732: {
2733: int m,ret,w;
2734: GF2N prev;
2735: UP2 r;
2736:
1.2 noro 2737: m = ZTOS((Q)ARG0(arg));
1.1 noro 2738: prev = (GF2N)ARG1(arg);
2739: if ( !prev ) {
2740: w = (m>>5)+1; NEWUP2(r,w); r->w = 0;
2741: bzero((char *)r->b,w*sizeof(unsigned int));
2742: } else {
2743: r = prev->body;
2744: if ( degup2(r) != m ) {
2745: w = (m>>5)+1; NEWUP2(r,w); r->w = 0;
2746: bzero((char *)r->b,w*sizeof(unsigned int));
2747: }
2748: }
2749: ret = _generate_irreducible_polynomial(r,m);
2750: if ( ret == 0 )
2751: MKGF2N(r,*rp);
2752: else
2753: *rp = 0;
2754: }
2755:
2756: void Pirredpoly_up2(NODE arg,GF2N *rp)
2757: {
2758: int m,ret,w;
2759: GF2N prev;
2760: UP2 r;
2761:
1.2 noro 2762: m = ZTOS((Q)ARG0(arg));
1.1 noro 2763: prev = (GF2N)ARG1(arg);
2764: if ( !prev ) {
2765: w = (m>>5)+1; NEWUP2(r,w); r->w = 0;
2766: bzero((char *)r->b,w*sizeof(unsigned int));
2767: } else {
2768: r = prev->body;
2769: if ( degup2(r) != m ) {
2770: w = (m>>5)+1; NEWUP2(r,w); r->w = 0;
2771: bzero((char *)r->b,w*sizeof(unsigned int));
2772: }
2773: }
2774: ret = _generate_good_irreducible_polynomial(r,m);
2775: if ( ret == 0 )
2776: MKGF2N(r,*rp);
2777: else
2778: *rp = 0;
2779: }
2780:
2781: void Pmat_swap_row_destructive(NODE arg, MAT *m)
2782: {
2783: int i1,i2;
2784: pointer *t;
2785: MAT mat;
2786:
2787: asir_assert(ARG0(arg),O_MAT,"mat_swap_row_destructive");
2788: asir_assert(ARG1(arg),O_N,"mat_swap_row_destructive");
2789: asir_assert(ARG2(arg),O_N,"mat_swap_row_destructive");
2790: mat = (MAT)ARG0(arg);
1.2 noro 2791: i1 = ZTOS((Q)ARG1(arg));
2792: i2 = ZTOS((Q)ARG2(arg));
1.1 noro 2793: if ( i1 < 0 || i2 < 0 || i1 >= mat->row || i2 >= mat->row )
2794: error("mat_swap_row_destructive : Out of range");
2795: t = mat->body[i1];
2796: mat->body[i1] = mat->body[i2];
2797: mat->body[i2] = t;
2798: *m = mat;
2799: }
2800:
2801: void Pmat_swap_col_destructive(NODE arg, MAT *m)
2802: {
2803: int j1,j2,i,n;
2804: pointer *mi;
2805: pointer t;
2806: MAT mat;
2807:
2808: asir_assert(ARG0(arg),O_MAT,"mat_swap_col_destructive");
2809: asir_assert(ARG1(arg),O_N,"mat_swap_col_destructive");
2810: asir_assert(ARG2(arg),O_N,"mat_swap_col_destructive");
2811: mat = (MAT)ARG0(arg);
1.2 noro 2812: j1 = ZTOS((Q)ARG1(arg));
2813: j2 = ZTOS((Q)ARG2(arg));
1.1 noro 2814: if ( j1 < 0 || j2 < 0 || j1 >= mat->col || j2 >= mat->col )
2815: error("mat_swap_col_destructive : Out of range");
2816: n = mat->row;
2817: for ( i = 0; i < n; i++ ) {
2818: mi = mat->body[i];
2819: t = mi[j1]; mi[j1] = mi[j2]; mi[j2] = t;
2820: }
2821: *m = mat;
2822: }
2823: /*
2824: * f = type 'type' normal polynomial of degree m if exists
2825: * IEEE P1363 A.7.2
2826: *
2827: * return value : 0 --- exists
2828: * 1 --- does not exist
2829: * -1 --- failure (memory allocation error)
2830: */
2831:
2832: int generate_ONB_polynomial(UP2 *rp,int m,int type)
2833: {
2834: int i,r;
2835: int w;
2836: UP2 f,f0,f1,f2,t;
2837:
2838: w = (m>>5)+1;
2839: switch ( type ) {
2840: case 1:
2841: if ( !TypeT_NB_check(m,1) ) return 1;
2842: NEWUP2(f,w); *rp = f; f->w = w;
2843: /* set all the bits */
2844: for ( i = 0; i < w; i++ )
2845: f->b[i] = 0xffffffff;
2846: /* mask the top word if necessary */
2847: if ( r = (m+1)&31 )
2848: f->b[w-1] &= (1<<r)-1;
2849: return 0;
2850: break;
2851: case 2:
2852: if ( !TypeT_NB_check(m,2) ) return 1;
2853: NEWUP2(f,w); *rp = f;
2854: W_NEWUP2(f0,w);
2855: W_NEWUP2(f1,w);
2856: W_NEWUP2(f2,w);
2857:
2858: /* recursion for genrating Type II normal polynomial */
2859:
2860: /* f0 = 1, f1 = t+1 */
2861: f0->w = 1; f0->b[0] = 1;
2862: f1->w = 1; f1->b[0] = 3;
2863: for ( i = 2; i <= m; i++ ) {
2864: /* f2 = t*f1+f0 */
2865: _bshiftup2(f1,-1,f2);
2866: _addup2_destructive(f2,f0);
2867: /* cyclic change of the variables */
2868: t = f0; f0 = f1; f1 = f2; f2 = t;
2869: }
2870: _copyup2(f1,f);
2871: return 0;
2872: break;
2873: default:
2874: return -1;
2875: break;
2876: }
2877: }
2878:
2879: /*
2880: * f = an irreducible trinomial or pentanomial of degree d 'after' f
2881: * return value : 0 --- exists
2882: * 1 --- does not exist (exhaustion)
2883: */
2884:
2885: int _generate_irreducible_polynomial(UP2 f,int d)
2886: {
2887: int ret,i,j,k,nz,i0,j0,k0;
2888: int w;
2889: unsigned int *fd;
2890:
2891: /*
2892: * if f = x^d+x^i+1 then i0 <- i, j0 <- 0, k0 <-0.
2893: * if f = x^d+x^k+x^j+x^i+1 (k>j>i) then i0 <- i, j0 <- j, k0 <-k.
2894: * otherwise i0,j0,k0 is set to 0.
2895: */
2896:
2897: fd = f->b;
2898: w = (d>>5)+1;
2899: if ( f->w && (d==degup2(f)) ) {
2900: for ( nz = 0, i = d; i >= 0; i-- )
2901: if ( fd[i>>5]&(1<<(i&31)) ) nz++;
2902: switch ( nz ) {
2903: case 3:
2904: for ( i0 = 1; !(fd[i0>>5]&(1<<(i0&31))) ; i0++ );
2905: /* reset i0-th bit */
2906: fd[i0>>5] &= ~(1<<(i0&31));
2907: j0 = k0 = 0;
2908: break;
2909: case 5:
2910: for ( i0 = 1; !(fd[i0>>5]&(1<<(i0&31))) ; i0++ );
2911: /* reset i0-th bit */
2912: fd[i0>>5] &= ~(1<<(i0&31));
2913: for ( j0 = i0+1; !(fd[j0>>5]&(1<<(j0&31))) ; j0++ );
2914: /* reset j0-th bit */
2915: fd[j0>>5] &= ~(1<<(j0&31));
2916: for ( k0 = j0+1; !(fd[k0>>5]&(1<<(k0&31))) ; k0++ );
2917: /* reset k0-th bit */
2918: fd[k0>>5] &= ~(1<<(k0&31));
2919: break;
2920: default:
2921: f->w = 0; break;
2922: }
2923: } else
2924: f->w = 0;
2925:
2926: if ( !f->w ) {
2927: fd = f->b;
2928: f->w = w; fd[0] |= 1; fd[d>>5] |= (1<<(d&31));
2929: i0 = j0 = k0 = 0;
2930: }
2931: /* if j0 > 0 then f is already a pentanomial */
2932: if ( j0 > 0 ) goto PENTA;
2933:
2934: /* searching for an irreducible trinomial */
2935:
2936: for ( i = 1; 2*i <= d; i++ ) {
2937: /* skip the polynomials 'before' f */
2938: if ( i < i0 ) continue;
2939: if ( i == i0 ) { i0 = 0; continue; }
2940: /* set i-th bit */
2941: fd[i>>5] |= (1<<(i&31));
2942: ret = irredcheck_dddup2(f);
2943: if ( ret == 1 ) return 0;
2944: /* reset i-th bit */
2945: fd[i>>5] &= ~(1<<(i&31));
2946: }
2947:
2948: /* searching for an irreducible pentanomial */
2949: PENTA:
2950: for ( i = 1; i < d; i++ ) {
2951: /* skip the polynomials 'before' f */
2952: if ( i < i0 ) continue;
2953: if ( i == i0 ) i0 = 0;
2954: /* set i-th bit */
2955: fd[i>>5] |= (1<<(i&31));
2956: for ( j = i+1; j < d; j++ ) {
2957: /* skip the polynomials 'before' f */
2958: if ( j < j0 ) continue;
2959: if ( j == j0 ) j0 = 0;
2960: /* set j-th bit */
2961: fd[j>>5] |= (1<<(j&31));
2962: for ( k = j+1; k < d; k++ ) {
2963: /* skip the polynomials 'before' f */
2964: if ( k < k0 ) continue;
2965: else if ( k == k0 ) { k0 = 0; continue; }
2966: /* set k-th bit */
2967: fd[k>>5] |= (1<<(k&31));
2968: ret = irredcheck_dddup2(f);
2969: if ( ret == 1 ) return 0;
2970: /* reset k-th bit */
2971: fd[k>>5] &= ~(1<<(k&31));
2972: }
2973: /* reset j-th bit */
2974: fd[j>>5] &= ~(1<<(j&31));
2975: }
2976: /* reset i-th bit */
2977: fd[i>>5] &= ~(1<<(i&31));
2978: }
2979: /* exhausted */
2980: return 1;
2981: }
2982:
2983: /*
2984: * f = an irreducible trinomial or pentanomial of degree d 'after' f
2985: *
2986: * searching strategy:
2987: * trinomial x^d+x^i+1:
2988: * i is as small as possible.
2989: * trinomial x^d+x^i+x^j+x^k+1:
2990: * i is as small as possible.
2991: * For such i, j is as small as possible.
2992: * For such i and j, 'k' is as small as possible.
2993: *
2994: * return value : 0 --- exists
2995: * 1 --- does not exist (exhaustion)
2996: */
2997:
2998: int _generate_good_irreducible_polynomial(UP2 f,int d)
2999: {
3000: int ret,i,j,k,nz,i0,j0,k0;
3001: int w;
3002: unsigned int *fd;
3003:
3004: /*
3005: * if f = x^d+x^i+1 then i0 <- i, j0 <- 0, k0 <-0.
3006: * if f = x^d+x^k+x^j+x^i+1 (k>j>i) then i0 <- i, j0 <- j, k0 <-k.
3007: * otherwise i0,j0,k0 is set to 0.
3008: */
3009:
3010: fd = f->b;
3011: w = (d>>5)+1;
3012: if ( f->w && (d==degup2(f)) ) {
3013: for ( nz = 0, i = d; i >= 0; i-- )
3014: if ( fd[i>>5]&(1<<(i&31)) ) nz++;
3015: switch ( nz ) {
3016: case 3:
3017: for ( i0 = 1; !(fd[i0>>5]&(1<<(i0&31))) ; i0++ );
3018: /* reset i0-th bit */
3019: fd[i0>>5] &= ~(1<<(i0&31));
3020: j0 = k0 = 0;
3021: break;
3022: case 5:
3023: for ( i0 = 1; !(fd[i0>>5]&(1<<(i0&31))) ; i0++ );
3024: /* reset i0-th bit */
3025: fd[i0>>5] &= ~(1<<(i0&31));
3026: for ( j0 = i0+1; !(fd[j0>>5]&(1<<(j0&31))) ; j0++ );
3027: /* reset j0-th bit */
3028: fd[j0>>5] &= ~(1<<(j0&31));
3029: for ( k0 = j0+1; !(fd[k0>>5]&(1<<(k0&31))) ; k0++ );
3030: /* reset k0-th bit */
3031: fd[k0>>5] &= ~(1<<(k0&31));
3032: break;
3033: default:
3034: f->w = 0; break;
3035: }
3036: } else
3037: f->w = 0;
3038:
3039: if ( !f->w ) {
3040: fd = f->b;
3041: f->w = w; fd[0] |= 1; fd[d>>5] |= (1<<(d&31));
3042: i0 = j0 = k0 = 0;
3043: }
3044: /* if j0 > 0 then f is already a pentanomial */
3045: if ( j0 > 0 ) goto PENTA;
3046:
3047: /* searching for an irreducible trinomial */
3048:
3049: for ( i = 1; 2*i <= d; i++ ) {
3050: /* skip the polynomials 'before' f */
3051: if ( i < i0 ) continue;
3052: if ( i == i0 ) { i0 = 0; continue; }
3053: /* set i-th bit */
3054: fd[i>>5] |= (1<<(i&31));
3055: ret = irredcheck_dddup2(f);
3056: if ( ret == 1 ) return 0;
3057: /* reset i-th bit */
3058: fd[i>>5] &= ~(1<<(i&31));
3059: }
3060:
3061: /* searching for an irreducible pentanomial */
3062: PENTA:
3063: for ( i = 3; i < d; i++ ) {
3064: /* skip the polynomials 'before' f */
3065: if ( i < i0 ) continue;
3066: if ( i == i0 ) i0 = 0;
3067: /* set i-th bit */
3068: fd[i>>5] |= (1<<(i&31));
3069: for ( j = 2; j < i; j++ ) {
3070: /* skip the polynomials 'before' f */
3071: if ( j < j0 ) continue;
3072: if ( j == j0 ) j0 = 0;
3073: /* set j-th bit */
3074: fd[j>>5] |= (1<<(j&31));
3075: for ( k = 1; k < j; k++ ) {
3076: /* skip the polynomials 'before' f */
3077: if ( k < k0 ) continue;
3078: else if ( k == k0 ) { k0 = 0; continue; }
3079: /* set k-th bit */
3080: fd[k>>5] |= (1<<(k&31));
3081: ret = irredcheck_dddup2(f);
3082: if ( ret == 1 ) return 0;
3083: /* reset k-th bit */
3084: fd[k>>5] &= ~(1<<(k&31));
3085: }
3086: /* reset j-th bit */
3087: fd[j>>5] &= ~(1<<(j&31));
3088: }
3089: /* reset i-th bit */
3090: fd[i>>5] &= ~(1<<(i&31));
3091: }
3092: /* exhausted */
3093: return 1;
3094: }
3095:
3096: void printqmat(Q **mat,int row,int col)
3097: {
3098: int i,j;
3099:
3100: for ( i = 0; i < row; i++ ) {
3101: for ( j = 0; j < col; j++ ) {
3102: printnum((Num)mat[i][j]); printf(" ");
3103: }
3104: printf("\n");
3105: }
3106: }
3107:
3108: void printimat(int **mat,int row,int col)
3109: {
3110: int i,j;
3111:
3112: for ( i = 0; i < row; i++ ) {
3113: for ( j = 0; j < col; j++ ) {
3114: printf("%d ",mat[i][j]);
3115: }
3116: printf("\n");
3117: }
3118: }
3119:
3120: void Pnd_det(NODE arg,P *rp)
3121: {
3122: if ( argc(arg) == 1 )
3123: nd_det(0,ARG0(arg),rp);
3124: else
1.2 noro 3125: nd_det(ZTOS((Q)ARG1(arg)),ARG0(arg),rp);
1.1 noro 3126: }
3127:
3128: void Pmat_col(NODE arg,VECT *rp)
3129: {
3130: int i,j,n;
3131: MAT mat;
3132: VECT vect;
3133:
3134: asir_assert(ARG0(arg),O_MAT,"mat_col");
3135: asir_assert(ARG1(arg),O_N,"mat_col");
3136: mat = (MAT)ARG0(arg);
1.2 noro 3137: j = ZTOS((Q)ARG1(arg));
1.1 noro 3138: if ( j < 0 || j >= mat->col) {
3139: error("mat_col : Out of range");
3140: }
3141: n = mat->row;
3142: MKVECT(vect,n);
3143: for(i=0; i<n; i++) {
3144: BDY(vect)[i] = BDY(mat)[i][j];
3145: }
3146: *rp = vect;
3147: }
3148:
3149: NODE triangleq(NODE e)
3150: {
3151: int n,i,k;
3152: V v;
3153: VL vl;
3154: P *p;
3155: NODE r,r1;
3156:
3157: n = length(e);
3158: p = (P *)MALLOC(n*sizeof(P));
3159: for ( i = 0; i < n; i++, e = NEXT(e) ) p[i] = (P)BDY(e);
3160: i = 0;
3161: while ( 1 ) {
3162: for ( ; i < n && !p[i]; i++ );
3163: if ( i == n ) break;
3164: if ( OID(p[i]) == O_N ) return 0;
3165: v = p[i]->v;
3166: for ( k = i+1; k < n; k++ )
3167: if ( p[k] ) {
3168: if ( OID(p[k]) == O_N ) return 0;
3169: if ( p[k]->v == v ) p[k] = 0;
3170: }
3171: i++;
3172: }
3173: for ( r = 0, i = 0; i < n; i++ ) {
3174: if ( p[i] ) {
3175: MKNODE(r1,p[i],r); r = r1;
3176: }
3177: }
3178: return r;
3179: }
3180:
3181: void Ptriangleq(NODE arg,LIST *rp)
3182: {
3183: NODE ret;
3184:
3185: asir_assert(ARG0(arg),O_LIST,"sparseleq");
3186: ret = triangleq(BDY((LIST)ARG0(arg)));
3187: MKLIST(*rp,ret);
3188: }
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