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