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