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