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