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