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