![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to array.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib2 / asir2000 / builtin|
Return to array.c CVS log | Up to [local] / OpenXM_contrib2 / asir2000 / builtin |
1.1 ! noro 1: /* $OpenXM: OpenXM/src/asir99/builtin/array.c,v 1.2 1999/11/23 07:14:14 noro Exp $ */
! 2: #include "ca.h"
! 3: #include "base.h"
! 4: #include "parse.h"
! 5: #include "inline.h"
! 6: /*
! 7: #undef DMAR
! 8: #define DMAR(a1,a2,a3,d,r) (r)=dmar(a1,a2,a3,d);
! 9: */
! 10:
! 11: extern int Print; /* XXX */
! 12:
! 13: void solve_by_lu_gfmmat(GFMMAT,unsigned int,unsigned int *,unsigned int *);
! 14: int lu_gfmmat(GFMMAT,unsigned int,int *);
! 15: void mat_to_gfmmat(MAT,unsigned int,GFMMAT *);
! 16:
! 17: int generic_gauss_elim_mod(int **,int,int,int,int *);
! 18: int generic_gauss_elim(MAT ,MAT *,Q *,int **,int **);
! 19:
! 20: int gauss_elim_mod(int **,int,int,int);
! 21: int gauss_elim_mod1(int **,int,int,int);
! 22: int gauss_elim_geninv_mod(unsigned int **,int,int,int);
! 23: int gauss_elim_geninv_mod_swap(unsigned int **,int,int,unsigned int,unsigned int ***,int **);
! 24: void Pnewvect(), Pnewmat(), Psepvect(), Psize(), Pdet(), Pleqm(), Pleqm1(), Pgeninvm();
! 25:
! 26: void Pgeneric_gauss_elim_mod();
! 27:
! 28: void Pmat_to_gfmmat(),Plu_gfmmat(),Psolve_by_lu_gfmmat();
! 29: void Pgeninvm_swap(), Premainder(), Psremainder(), Pvtol();
! 30: void sepvect();
! 31: void Pmulmat_gf2n();
! 32: void Pbconvmat_gf2n();
! 33: void Pmul_vect_mat_gf2n();
! 34: void PNBmul_gf2n();
! 35: void Pmul_mat_vect_int();
! 36: void Psepmat_destructive();
! 37: void Px962_irredpoly_up2();
! 38: void Pirredpoly_up2();
! 39: void Pnbpoly_up2();
! 40: void Pqsort();
! 41:
! 42: struct ftab array_tab[] = {
! 43: {"solve_by_lu_gfmmat",Psolve_by_lu_gfmmat,4},
! 44: {"lu_gfmmat",Plu_gfmmat,2},
! 45: {"mat_to_gfmmat",Pmat_to_gfmmat,2},
! 46: {"generic_gauss_elim_mod",Pgeneric_gauss_elim_mod,2},
! 47: {"newvect",Pnewvect,-2},
! 48: {"newmat",Pnewmat,-3},
! 49: {"sepmat_destructive",Psepmat_destructive,2},
! 50: {"sepvect",Psepvect,2},
! 51: {"qsort",Pqsort,-2},
! 52: {"vtol",Pvtol,1},
! 53: {"size",Psize,1},
! 54: {"det",Pdet,-2},
! 55: {"leqm",Pleqm,2},
! 56: {"leqm1",Pleqm1,2},
! 57: {"geninvm",Pgeninvm,2},
! 58: {"geninvm_swap",Pgeninvm_swap,2},
! 59: {"remainder",Premainder,2},
! 60: {"sremainder",Psremainder,2},
! 61: {"mulmat_gf2n",Pmulmat_gf2n,1},
! 62: {"bconvmat_gf2n",Pbconvmat_gf2n,-4},
! 63: {"mul_vect_mat_gf2n",Pmul_vect_mat_gf2n,2},
! 64: {"mul_mat_vect_int",Pmul_mat_vect_int,2},
! 65: {"nbmul_gf2n",PNBmul_gf2n,3},
! 66: {"x962_irredpoly_up2",Px962_irredpoly_up2,2},
! 67: {"irredpoly_up2",Pirredpoly_up2,2},
! 68: {"nbpoly_up2",Pnbpoly_up2,2},
! 69: {0,0,0},
! 70: };
! 71:
! 72: int comp_obj(a,b)
! 73: Obj *a,*b;
! 74: {
! 75: return arf_comp(CO,*a,*b);
! 76: }
! 77:
! 78: static FUNC generic_comp_obj_func;
! 79: static NODE generic_comp_obj_arg;
! 80:
! 81: int generic_comp_obj(a,b)
! 82: Obj *a,*b;
! 83: {
! 84: Q r;
! 85:
! 86: BDY(generic_comp_obj_arg)=(pointer)(*a);
! 87: BDY(NEXT(generic_comp_obj_arg))=(pointer)(*b);
! 88: r = (Q)bevalf(generic_comp_obj_func,generic_comp_obj_arg);
! 89: if ( !r )
! 90: return 0;
! 91: else
! 92: return SGN(r)>0?1:-1;
! 93: }
! 94:
! 95:
! 96: void Pqsort(arg,rp)
! 97: NODE arg;
! 98: VECT *rp;
! 99: {
! 100: VECT vect;
! 101: char buf[BUFSIZ];
! 102: char *fname;
! 103: NODE n;
! 104: P p;
! 105: V v;
! 106:
! 107: asir_assert(ARG0(arg),O_VECT,"qsort");
! 108: vect = (VECT)ARG0(arg);
! 109: if ( argc(arg) == 1 )
! 110: qsort(BDY(vect),vect->len,sizeof(Obj),(int (*)(const void *,const void *))comp_obj);
! 111: else {
! 112: p = (P)ARG1(arg);
! 113: if ( !p || OID(p)!=2 )
! 114: error("qsort : invalid argument");
! 115: v = VR(p);
! 116: if ( (int)v->attr != V_SR )
! 117: error("qsort : no such function");
! 118: generic_comp_obj_func = (FUNC)v->priv;
! 119: MKNODE(n,0,0); MKNODE(generic_comp_obj_arg,0,n);
! 120: qsort(BDY(vect),vect->len,sizeof(Obj),(int (*)(const void *,const void *))generic_comp_obj);
! 121: }
! 122: *rp = vect;
! 123: }
! 124:
! 125: void PNBmul_gf2n(arg,rp)
! 126: NODE arg;
! 127: GF2N *rp;
! 128: {
! 129: GF2N a,b;
! 130: GF2MAT mat;
! 131: int n,w;
! 132: unsigned int *ab,*bb;
! 133: UP2 r;
! 134:
! 135: a = (GF2N)ARG0(arg);
! 136: b = (GF2N)ARG1(arg);
! 137: mat = (GF2MAT)ARG2(arg);
! 138: if ( !a || !b )
! 139: *rp = 0;
! 140: else {
! 141: n = mat->row;
! 142: w = (n+BSH-1)/BSH;
! 143:
! 144: ab = (unsigned int *)ALLOCA(w*sizeof(unsigned int));
! 145: bzero((char *)ab,w*sizeof(unsigned int));
! 146: bcopy(a->body->b,ab,(a->body->w)*sizeof(unsigned int));
! 147:
! 148: bb = (unsigned int *)ALLOCA(w*sizeof(unsigned int));
! 149: bzero((char *)bb,w*sizeof(unsigned int));
! 150: bcopy(b->body->b,bb,(b->body->w)*sizeof(unsigned int));
! 151:
! 152: NEWUP2(r,w);
! 153: bzero((char *)r->b,w*sizeof(unsigned int));
! 154: mul_nb(mat,ab,bb,r->b);
! 155: r->w = w;
! 156: _adjup2(r);
! 157: if ( !r->w )
! 158: *rp = 0;
! 159: else
! 160: MKGF2N(r,*rp);
! 161: }
! 162: }
! 163:
! 164: void Pmul_vect_mat_gf2n(arg,rp)
! 165: NODE arg;
! 166: GF2N *rp;
! 167: {
! 168: GF2N a;
! 169: GF2MAT mat;
! 170: int n,w;
! 171: unsigned int *b;
! 172: UP2 r;
! 173:
! 174: a = (GF2N)ARG0(arg);
! 175: mat = (GF2MAT)ARG1(arg);
! 176: if ( !a )
! 177: *rp = 0;
! 178: else {
! 179: n = mat->row;
! 180: w = (n+BSH-1)/BSH;
! 181: b = (unsigned int *)ALLOCA(w*sizeof(unsigned int));
! 182: bzero((char *)b,w*sizeof(unsigned int));
! 183: bcopy(a->body->b,b,(a->body->w)*sizeof(unsigned int));
! 184: NEWUP2(r,w);
! 185: bzero((char *)r->b,w*sizeof(unsigned int));
! 186: mulgf2vectmat(mat->row,b,mat->body,r->b);
! 187: r->w = w;
! 188: _adjup2(r);
! 189: if ( !r->w )
! 190: *rp = 0;
! 191: else {
! 192: MKGF2N(r,*rp);
! 193: }
! 194: }
! 195: }
! 196:
! 197: void Pbconvmat_gf2n(arg,rp)
! 198: NODE arg;
! 199: LIST *rp;
! 200: {
! 201: P p0,p1;
! 202: int to;
! 203: GF2MAT p01,p10;
! 204: GF2N root;
! 205: NODE n0,n1;
! 206:
! 207: p0 = (P)ARG0(arg);
! 208: p1 = (P)ARG1(arg);
! 209: to = ARG2(arg)?1:0;
! 210: if ( argc(arg) == 4 ) {
! 211: root = (GF2N)ARG3(arg);
! 212: compute_change_of_basis_matrix_with_root(p0,p1,to,root,&p01,&p10);
! 213: } else
! 214: compute_change_of_basis_matrix(p0,p1,to,&p01,&p10);
! 215: MKNODE(n1,p10,0); MKNODE(n0,p01,n1);
! 216: MKLIST(*rp,n0);
! 217: }
! 218:
! 219: void Pmulmat_gf2n(arg,rp)
! 220: NODE arg;
! 221: GF2MAT *rp;
! 222: {
! 223: GF2MAT m;
! 224:
! 225: if ( !compute_multiplication_matrix((P)ARG0(arg),&m) )
! 226: error("mulmat_gf2n : input is not a normal polynomial");
! 227: *rp = m;
! 228: }
! 229:
! 230: void Psepmat_destructive(arg,rp)
! 231: NODE arg;
! 232: LIST *rp;
! 233: {
! 234: MAT mat,mat1;
! 235: int i,j,row,col;
! 236: Q **a,**a1;
! 237: Q ent;
! 238: N nm,mod,rem,quo;
! 239: int sgn;
! 240: NODE n0,n1;
! 241:
! 242: mat = (MAT)ARG0(arg); mod = NM((Q)ARG1(arg));
! 243: row = mat->row; col = mat->col;
! 244: MKMAT(mat1,row,col);
! 245: a = (Q **)mat->body; a1 = (Q **)mat1->body;
! 246: for ( i = 0; i < row; i++ )
! 247: for ( j = 0; j < col; j++ ) {
! 248: ent = a[i][j];
! 249: if ( !ent )
! 250: continue;
! 251: nm = NM(ent);
! 252: sgn = SGN(ent);
! 253: divn(nm,mod,&quo,&rem);
! 254: /* if ( quo != nm && rem != nm ) */
! 255: /* GC_free(nm); */
! 256: /* GC_free(ent); */
! 257: NTOQ(rem,sgn,a[i][j]); NTOQ(quo,sgn,a1[i][j]);
! 258: }
! 259: MKNODE(n1,mat1,0); MKNODE(n0,mat,n1);
! 260: MKLIST(*rp,n0);
! 261: }
! 262:
! 263: void Psepvect(arg,rp)
! 264: NODE arg;
! 265: VECT *rp;
! 266: {
! 267: sepvect((VECT)ARG0(arg),QTOS((Q)ARG1(arg)),rp);
! 268: }
! 269:
! 270: void sepvect(v,d,rp)
! 271: VECT v;
! 272: int d;
! 273: VECT *rp;
! 274: {
! 275: int i,j,k,n,q,q1,r;
! 276: pointer *pv,*pw,*pu;
! 277: VECT w,u;
! 278:
! 279: n = v->len;
! 280: if ( d > n )
! 281: d = n;
! 282: q = n/d; r = n%d; q1 = q+1;
! 283: MKVECT(w,d); *rp = w;
! 284: pv = BDY(v); pw = BDY(w); k = 0;
! 285: for ( i = 0; i < r; i++ ) {
! 286: MKVECT(u,q1); pw[i] = (pointer)u;
! 287: for ( pu = BDY(u), j = 0; j < q1; j++, k++ )
! 288: pu[j] = pv[k];
! 289: }
! 290: for ( ; i < d; i++ ) {
! 291: MKVECT(u,q); pw[i] = (pointer)u;
! 292: for ( pu = BDY(u), j = 0; j < q; j++, k++ )
! 293: pu[j] = pv[k];
! 294: }
! 295: }
! 296:
! 297: void Pnewvect(arg,rp)
! 298: NODE arg;
! 299: VECT *rp;
! 300: {
! 301: int len,i,r;
! 302: VECT vect;
! 303: pointer *vb;
! 304: LIST list;
! 305: NODE tn;
! 306:
! 307: asir_assert(ARG0(arg),O_N,"newvect");
! 308: len = QTOS((Q)ARG0(arg));
! 309: if ( len <= 0 )
! 310: error("newvect : invalid size");
! 311: MKVECT(vect,len);
! 312: if ( argc(arg) == 2 ) {
! 313: list = (LIST)ARG1(arg);
! 314: asir_assert(list,O_LIST,"newvect");
! 315: for ( r = 0, tn = BDY(list); tn; r++, tn = NEXT(tn) );
! 316: if ( r > len ) {
! 317: *rp = vect;
! 318: return;
! 319: }
! 320: for ( i = 0, tn = BDY(list), vb = BDY(vect); tn; i++, tn = NEXT(tn) )
! 321: vb[i] = (pointer)BDY(tn);
! 322: }
! 323: *rp = vect;
! 324: }
! 325:
! 326: void Pnewmat(arg,rp)
! 327: NODE arg;
! 328: MAT *rp;
! 329: {
! 330: int row,col;
! 331: int i,j,r,c;
! 332: NODE tn,sn;
! 333: MAT m;
! 334: pointer **mb;
! 335: LIST list;
! 336:
! 337: asir_assert(ARG0(arg),O_N,"newmat");
! 338: asir_assert(ARG1(arg),O_N,"newmat");
! 339: row = QTOS((Q)ARG0(arg)); col = QTOS((Q)ARG1(arg));
! 340: if ( row <= 0 || col <= 0 )
! 341: error("newmat : invalid size");
! 342: MKMAT(m,row,col);
! 343: if ( argc(arg) == 3 ) {
! 344: list = (LIST)ARG2(arg);
! 345: asir_assert(list,O_LIST,"newmat");
! 346: for ( r = 0, c = 0, tn = BDY(list); tn; r++, tn = NEXT(tn) ) {
! 347: for ( j = 0, sn = BDY((LIST)BDY(tn)); sn; j++, sn = NEXT(sn) );
! 348: c = MAX(c,j);
! 349: }
! 350: if ( (r > row) || (c > col) ) {
! 351: *rp = m;
! 352: return;
! 353: }
! 354: for ( i = 0, tn = BDY(list), mb = BDY(m); tn; i++, tn = NEXT(tn) ) {
! 355: asir_assert(BDY(tn),O_LIST,"newmat");
! 356: for ( j = 0, sn = BDY((LIST)BDY(tn)); sn; j++, sn = NEXT(sn) )
! 357: mb[i][j] = (pointer)BDY(sn);
! 358: }
! 359: }
! 360: *rp = m;
! 361: }
! 362:
! 363: void Pvtol(arg,rp)
! 364: NODE arg;
! 365: LIST *rp;
! 366: {
! 367: NODE n,n1;
! 368: VECT v;
! 369: pointer *a;
! 370: int len,i;
! 371:
! 372: asir_assert(ARG0(arg),O_VECT,"vtol");
! 373: v = (VECT)ARG0(arg); len = v->len; a = BDY(v);
! 374: for ( i = len - 1, n = 0; i >= 0; i-- ) {
! 375: MKNODE(n1,a[i],n); n = n1;
! 376: }
! 377: MKLIST(*rp,n);
! 378: }
! 379:
! 380: void Premainder(arg,rp)
! 381: NODE arg;
! 382: Obj *rp;
! 383: {
! 384: Obj a;
! 385: VECT v,w;
! 386: MAT m,l;
! 387: pointer *vb,*wb;
! 388: pointer **mb,**lb;
! 389: int id,i,j,n,row,col,t,smd,sgn;
! 390: Q md,q;
! 391:
! 392: a = (Obj)ARG0(arg); md = (Q)ARG1(arg);
! 393: if ( !a )
! 394: *rp = 0;
! 395: else {
! 396: id = OID(a);
! 397: switch ( id ) {
! 398: case O_N:
! 399: case O_P:
! 400: cmp(md,(P)a,(P *)rp); break;
! 401: case O_VECT:
! 402: smd = QTOS(md);
! 403: v = (VECT)a; n = v->len; vb = v->body;
! 404: MKVECT(w,n); wb = w->body;
! 405: for ( i = 0; i < n; i++ ) {
! 406: if ( q = (Q)vb[i] ) {
! 407: sgn = SGN(q); t = rem(NM(q),smd);
! 408: STOQ(t,q);
! 409: if ( q )
! 410: SGN(q) = sgn;
! 411: }
! 412: wb[i] = (pointer)q;
! 413: }
! 414: *rp = (Obj)w;
! 415: break;
! 416: case O_MAT:
! 417: m = (MAT)a; row = m->row; col = m->col; mb = m->body;
! 418: MKMAT(l,row,col); lb = l->body;
! 419: for ( i = 0; i < row; i++ )
! 420: for ( j = 0, vb = mb[i], wb = lb[i]; j < col; j++ )
! 421: cmp(md,(P)vb[j],(P *)&wb[j]);
! 422: *rp = (Obj)l;
! 423: break;
! 424: default:
! 425: error("remainder : invalid argument");
! 426: }
! 427: }
! 428: }
! 429:
! 430: void Psremainder(arg,rp)
! 431: NODE arg;
! 432: Obj *rp;
! 433: {
! 434: Obj a;
! 435: VECT v,w;
! 436: MAT m,l;
! 437: pointer *vb,*wb;
! 438: pointer **mb,**lb;
! 439: unsigned int t,smd;
! 440: int id,i,j,n,row,col;
! 441: Q md,q;
! 442:
! 443: a = (Obj)ARG0(arg); md = (Q)ARG1(arg);
! 444: if ( !a )
! 445: *rp = 0;
! 446: else {
! 447: id = OID(a);
! 448: switch ( id ) {
! 449: case O_N:
! 450: case O_P:
! 451: cmp(md,(P)a,(P *)rp); break;
! 452: case O_VECT:
! 453: smd = QTOS(md);
! 454: v = (VECT)a; n = v->len; vb = v->body;
! 455: MKVECT(w,n); wb = w->body;
! 456: for ( i = 0; i < n; i++ ) {
! 457: if ( q = (Q)vb[i] ) {
! 458: t = (unsigned int)rem(NM(q),smd);
! 459: if ( SGN(q) < 0 )
! 460: t = (smd - t) % smd;
! 461: UTOQ(t,q);
! 462: }
! 463: wb[i] = (pointer)q;
! 464: }
! 465: *rp = (Obj)w;
! 466: break;
! 467: case O_MAT:
! 468: m = (MAT)a; row = m->row; col = m->col; mb = m->body;
! 469: MKMAT(l,row,col); lb = l->body;
! 470: for ( i = 0; i < row; i++ )
! 471: for ( j = 0, vb = mb[i], wb = lb[i]; j < col; j++ )
! 472: cmp(md,(P)vb[j],(P *)&wb[j]);
! 473: *rp = (Obj)l;
! 474: break;
! 475: default:
! 476: error("remainder : invalid argument");
! 477: }
! 478: }
! 479: }
! 480:
! 481: void Psize(arg,rp)
! 482: NODE arg;
! 483: LIST *rp;
! 484: {
! 485:
! 486: int n,m;
! 487: Q q;
! 488: NODE t,s;
! 489:
! 490: if ( !ARG0(arg) )
! 491: t = 0;
! 492: else {
! 493: switch (OID(ARG0(arg))) {
! 494: case O_VECT:
! 495: n = ((VECT)ARG0(arg))->len;
! 496: STOQ(n,q); MKNODE(t,q,0);
! 497: break;
! 498: case O_MAT:
! 499: n = ((MAT)ARG0(arg))->row; m = ((MAT)ARG0(arg))->col;
! 500: STOQ(m,q); MKNODE(s,q,0); STOQ(n,q); MKNODE(t,q,s);
! 501: break;
! 502: default:
! 503: error("size : invalid argument"); break;
! 504: }
! 505: }
! 506: MKLIST(*rp,t);
! 507: }
! 508:
! 509: void Pdet(arg,rp)
! 510: NODE arg;
! 511: P *rp;
! 512: {
! 513: MAT m;
! 514: int n,i,j,mod;
! 515: P d;
! 516: P **mat,**w;
! 517:
! 518: m = (MAT)ARG0(arg);
! 519: asir_assert(m,O_MAT,"det");
! 520: if ( m->row != m->col )
! 521: error("det : non-square matrix");
! 522: else if ( argc(arg) == 1 )
! 523: detp(CO,(P **)BDY(m),m->row,rp);
! 524: else {
! 525: n = m->row; mod = QTOS((Q)ARG1(arg)); mat = (P **)BDY(m);
! 526: w = (P **)almat_pointer(n,n);
! 527: for ( i = 0; i < n; i++ )
! 528: for ( j = 0; j < n; j++ )
! 529: ptomp(mod,mat[i][j],&w[i][j]);
! 530: detmp(CO,mod,w,n,&d);
! 531: mptop(d,rp);
! 532: }
! 533: }
! 534:
! 535: /*
! 536: input : a row x col matrix A
! 537: A[I] <-> A[I][0]*x_0+A[I][1]*x_1+...
! 538:
! 539: output : [B,R,C]
! 540: B : a rank(A) x col-rank(A) matrix
! 541: R : a vector of length rank(A)
! 542: C : a vector of length col-rank(A)
! 543: B[I] <-> x_{R[I]}+B[I][0]x_{C[0]}+B[I][1]x_{C[1]}+...
! 544: */
! 545:
! 546: void Pgeneric_gauss_elim_mod(arg,rp)
! 547: NODE arg;
! 548: LIST *rp;
! 549: {
! 550: NODE n0;
! 551: MAT m,mat;
! 552: VECT rind,cind;
! 553: Q **tmat;
! 554: int **wmat;
! 555: Q *rib,*cib;
! 556: int *colstat;
! 557: Q q;
! 558: int md,i,j,k,l,row,col,t,n,rank;
! 559:
! 560: asir_assert(ARG0(arg),O_MAT,"generic_gauss_elim_mod");
! 561: asir_assert(ARG1(arg),O_N,"generic_gauss_elim_mod");
! 562: m = (MAT)ARG0(arg); md = QTOS((Q)ARG1(arg));
! 563: row = m->row; col = m->col; tmat = (Q **)m->body;
! 564: wmat = (int **)almat(row,col);
! 565: colstat = (int *)MALLOC_ATOMIC(col*sizeof(int));
! 566: for ( i = 0; i < row; i++ )
! 567: for ( j = 0; j < col; j++ )
! 568: if ( q = (Q)tmat[i][j] ) {
! 569: t = rem(NM(q),md);
! 570: if ( t && SGN(q) < 0 )
! 571: t = (md - t) % md;
! 572: wmat[i][j] = t;
! 573: } else
! 574: wmat[i][j] = 0;
! 575: rank = generic_gauss_elim_mod(wmat,row,col,md,colstat);
! 576:
! 577: MKMAT(mat,rank,col-rank);
! 578: tmat = (Q **)mat->body;
! 579: for ( i = 0; i < rank; i++ )
! 580: for ( j = k = 0; j < col; j++ )
! 581: if ( !colstat[j] ) {
! 582: UTOQ(wmat[i][j],tmat[i][k]); k++;
! 583: }
! 584:
! 585: MKVECT(rind,rank);
! 586: MKVECT(cind,col-rank);
! 587: rib = (Q *)rind->body; cib = (Q *)cind->body;
! 588: for ( j = k = l = 0; j < col; j++ )
! 589: if ( colstat[j] ) {
! 590: STOQ(j,rib[k]); k++;
! 591: } else {
! 592: STOQ(j,cib[l]); l++;
! 593: }
! 594: n0 = mknode(3,mat,rind,cind);
! 595: MKLIST(*rp,n0);
! 596: }
! 597:
! 598: void Pleqm(arg,rp)
! 599: NODE arg;
! 600: VECT *rp;
! 601: {
! 602: MAT m;
! 603: VECT vect;
! 604: pointer **mat;
! 605: Q *v;
! 606: Q q;
! 607: int **wmat;
! 608: int md,i,j,row,col,t,n,status;
! 609:
! 610: asir_assert(ARG0(arg),O_MAT,"leqm");
! 611: asir_assert(ARG1(arg),O_N,"leqm");
! 612: m = (MAT)ARG0(arg); md = QTOS((Q)ARG1(arg));
! 613: row = m->row; col = m->col; mat = m->body;
! 614: wmat = (int **)almat(row,col);
! 615: for ( i = 0; i < row; i++ )
! 616: for ( j = 0; j < col; j++ )
! 617: if ( q = (Q)mat[i][j] ) {
! 618: t = rem(NM(q),md);
! 619: if ( SGN(q) < 0 )
! 620: t = (md - t) % md;
! 621: wmat[i][j] = t;
! 622: } else
! 623: wmat[i][j] = 0;
! 624: status = gauss_elim_mod(wmat,row,col,md);
! 625: if ( status < 0 )
! 626: *rp = 0;
! 627: else if ( status > 0 )
! 628: *rp = (VECT)ONE;
! 629: else {
! 630: n = col - 1;
! 631: MKVECT(vect,n);
! 632: for ( i = 0, v = (Q *)vect->body; i < n; i++ ) {
! 633: t = (md-wmat[i][n])%md; STOQ(t,v[i]);
! 634: }
! 635: *rp = vect;
! 636: }
! 637: }
! 638:
! 639: int gauss_elim_mod(mat,row,col,md)
! 640: int **mat;
! 641: int row,col,md;
! 642: {
! 643: int i,j,k,inv,a,n;
! 644: int *t,*pivot;
! 645:
! 646: n = col - 1;
! 647: for ( j = 0; j < n; j++ ) {
! 648: for ( i = j; i < row && !mat[i][j]; i++ );
! 649: if ( i == row )
! 650: return 1;
! 651: if ( i != j ) {
! 652: t = mat[i]; mat[i] = mat[j]; mat[j] = t;
! 653: }
! 654: pivot = mat[j];
! 655: inv = invm(pivot[j],md);
! 656: for ( k = j; k <= n; k++ ) {
! 657: /* pivot[k] = dmar(pivot[k],inv,0,md); */
! 658: DMAR(pivot[k],inv,0,md,pivot[k])
! 659: }
! 660: for ( i = 0; i < row; i++ ) {
! 661: t = mat[i];
! 662: if ( i != j && (a = t[j]) )
! 663: for ( k = j, a = md - a; k <= n; k++ ) {
! 664: /* t[k] = dmar(pivot[k],a,t[k],md); */
! 665: DMAR(pivot[k],a,t[k],md,t[k])
! 666: }
! 667: }
! 668: }
! 669: for ( i = n; i < row && !mat[i][n]; i++ );
! 670: if ( i == row )
! 671: return 0;
! 672: else
! 673: return -1;
! 674: }
! 675:
! 676: struct oEGT eg_mod,eg_elim,eg_chrem,eg_gschk,eg_intrat,eg_symb;
! 677:
! 678: int generic_gauss_elim(mat,nm,dn,rindp,cindp)
! 679: MAT mat;
! 680: MAT *nm;
! 681: Q *dn;
! 682: int **rindp,**cindp;
! 683: {
! 684: int **wmat;
! 685: Q **bmat;
! 686: N **tmat;
! 687: Q *bmi;
! 688: N *tmi;
! 689: Q q;
! 690: int *wmi;
! 691: int *colstat,*wcolstat,*rind,*cind;
! 692: int row,col,ind,md,i,j,k,l,t,t1,rank,rank0,inv;
! 693: N m1,m2,m3,s,u;
! 694: MAT r,crmat;
! 695: struct oEGT tmp0,tmp1;
! 696: struct oEGT eg_mod_split,eg_elim_split,eg_chrem_split;
! 697: struct oEGT eg_intrat_split,eg_gschk_split;
! 698: int ret;
! 699:
! 700: init_eg(&eg_mod_split); init_eg(&eg_chrem_split);
! 701: init_eg(&eg_elim_split); init_eg(&eg_intrat_split);
! 702: init_eg(&eg_gschk_split);
! 703: bmat = (Q **)mat->body;
! 704: row = mat->row; col = mat->col;
! 705: wmat = (int **)almat(row,col);
! 706: colstat = (int *)MALLOC_ATOMIC(col*sizeof(int));
! 707: wcolstat = (int *)MALLOC_ATOMIC(col*sizeof(int));
! 708: for ( ind = 0; ; ind++ ) {
! 709: if ( Print )
! 710: fprintf(asir_out,".");
! 711: md = lprime[ind];
! 712: get_eg(&tmp0);
! 713: for ( i = 0; i < row; i++ )
! 714: for ( j = 0, bmi = bmat[i], wmi = wmat[i]; j < col; j++ )
! 715: if ( q = (Q)bmi[j] ) {
! 716: t = rem(NM(q),md);
! 717: if ( t && SGN(q) < 0 )
! 718: t = (md - t) % md;
! 719: wmi[j] = t;
! 720: } else
! 721: wmi[j] = 0;
! 722: get_eg(&tmp1);
! 723: add_eg(&eg_mod,&tmp0,&tmp1);
! 724: add_eg(&eg_mod_split,&tmp0,&tmp1);
! 725: get_eg(&tmp0);
! 726: rank = generic_gauss_elim_mod(wmat,row,col,md,wcolstat);
! 727: get_eg(&tmp1);
! 728: add_eg(&eg_elim,&tmp0,&tmp1);
! 729: add_eg(&eg_elim_split,&tmp0,&tmp1);
! 730: if ( !ind ) {
! 731: RESET:
! 732: UTON(md,m1);
! 733: rank0 = rank;
! 734: bcopy(wcolstat,colstat,col*sizeof(int));
! 735: MKMAT(crmat,rank,col-rank);
! 736: MKMAT(r,rank,col-rank); *nm = r;
! 737: tmat = (N **)crmat->body;
! 738: for ( i = 0; i < rank; i++ )
! 739: for ( j = k = 0, tmi = tmat[i], wmi = wmat[i]; j < col; j++ )
! 740: if ( !colstat[j] ) {
! 741: UTON(wmi[j],tmi[k]); k++;
! 742: }
! 743: } else {
! 744: if ( rank < rank0 ) {
! 745: if ( Print )
! 746: fprintf(asir_out,"lower rank matrix; continuing...\n");
! 747: continue;
! 748: } else if ( rank > rank0 ) {
! 749: if ( Print )
! 750: fprintf(asir_out,"higher rank matrix; resetting...\n");
! 751: goto RESET;
! 752: } else {
! 753: for ( j = 0; (j<col) && (colstat[j]==wcolstat[j]); j++ );
! 754: if ( j < col ) {
! 755: if ( Print )
! 756: fprintf(asir_out,"inconsitent colstat; resetting...\n");
! 757: goto RESET;
! 758: }
! 759: }
! 760:
! 761: get_eg(&tmp0);
! 762: inv = invm(rem(m1,md),md);
! 763: UTON(md,m2); muln(m1,m2,&m3);
! 764: for ( i = 0; i < rank; i++ )
! 765: for ( j = k = 0, tmi = tmat[i], wmi = wmat[i]; j < col; j++ )
! 766: if ( !colstat[j] ) {
! 767: if ( tmi[k] ) {
! 768: /* f3 = f1+m1*(m1 mod m2)^(-1)*(f2 - f1 mod m2) */
! 769: t = rem(tmi[k],md);
! 770: if ( wmi[j] >= t )
! 771: t = wmi[j]-t;
! 772: else
! 773: t = md-(t-wmi[j]);
! 774: DMAR(t,inv,0,md,t1)
! 775: UTON(t1,u);
! 776: muln(m1,u,&s);
! 777: addn(tmi[k],s,&u); tmi[k] = u;
! 778: } else if ( wmi[j] ) {
! 779: /* f3 = m1*(m1 mod m2)^(-1)*f2 */
! 780: DMAR(wmi[j],inv,0,md,t)
! 781: UTON(t,u);
! 782: muln(m1,u,&s); tmi[k] = s;
! 783: }
! 784: k++;
! 785: }
! 786: m1 = m3;
! 787: get_eg(&tmp1);
! 788: add_eg(&eg_chrem,&tmp0,&tmp1);
! 789: add_eg(&eg_chrem_split,&tmp0,&tmp1);
! 790:
! 791: get_eg(&tmp0);
! 792: ret = intmtoratm(crmat,m1,*nm,dn);
! 793: get_eg(&tmp1);
! 794: add_eg(&eg_intrat,&tmp0,&tmp1);
! 795: add_eg(&eg_intrat_split,&tmp0,&tmp1);
! 796: if ( ret ) {
! 797: *rindp = rind = (int *)MALLOC_ATOMIC(rank*sizeof(int));
! 798: *cindp = cind = (int *)MALLOC_ATOMIC((col-rank)*sizeof(int));
! 799: for ( j = k = l = 0; j < col; j++ )
! 800: if ( colstat[j] )
! 801: rind[k++] = j;
! 802: else
! 803: cind[l++] = j;
! 804: get_eg(&tmp0);
! 805: if ( gensolve_check(mat,*nm,*dn,rind,cind) )
! 806: get_eg(&tmp1);
! 807: add_eg(&eg_gschk,&tmp0,&tmp1);
! 808: add_eg(&eg_gschk_split,&tmp0,&tmp1);
! 809: if ( Print ) {
! 810: print_eg("Mod",&eg_mod_split);
! 811: print_eg("Elim",&eg_elim_split);
! 812: print_eg("ChRem",&eg_chrem_split);
! 813: print_eg("IntRat",&eg_intrat_split);
! 814: print_eg("Check",&eg_gschk_split);
! 815: fflush(asir_out);
! 816: }
! 817: return rank;
! 818: }
! 819: }
! 820: }
! 821: }
! 822:
! 823: int f4_nocheck;
! 824:
! 825: int gensolve_check(mat,nm,dn,rind,cind)
! 826: MAT mat,nm;
! 827: Q dn;
! 828: int *rind,*cind;
! 829: {
! 830: int row,col,rank,clen,i,j,k,l;
! 831: Q s,t,u;
! 832: Q *w;
! 833: Q *mati,*nmk;
! 834:
! 835: if ( f4_nocheck )
! 836: return 1;
! 837: row = mat->row; col = mat->col;
! 838: rank = nm->row; clen = nm->col;
! 839: w = (Q *)MALLOC(clen*sizeof(Q));
! 840: for ( i = 0; i < row; i++ ) {
! 841: mati = (Q *)mat->body[i];
! 842: #if 1
! 843: bzero(w,clen*sizeof(Q));
! 844: for ( k = 0; k < rank; k++ )
! 845: for ( l = 0, nmk = (Q *)nm->body[k]; l < clen; l++ ) {
! 846: mulq(mati[rind[k]],nmk[l],&t);
! 847: addq(w[l],t,&s); w[l] = s;
! 848: }
! 849: for ( j = 0; j < clen; j++ ) {
! 850: mulq(dn,mati[cind[j]],&t);
! 851: if ( cmpq(w[j],t) )
! 852: break;
! 853: }
! 854: #else
! 855: for ( j = 0; j < clen; j++ ) {
! 856: for ( k = 0, s = 0; k < rank; k++ ) {
! 857: mulq(mati[rind[k]],nm->body[k][j],&t);
! 858: addq(s,t,&u); s = u;
! 859: }
! 860: mulq(dn,mati[cind[j]],&t);
! 861: if ( cmpq(s,t) )
! 862: break;
! 863: }
! 864: #endif
! 865: if ( j != clen )
! 866: break;
! 867: }
! 868: if ( i != row )
! 869: return 0;
! 870: else
! 871: return 1;
! 872: }
! 873:
! 874: /* assuming 0 < c < m */
! 875:
! 876: int inttorat(c,m,b,sgnp,nmp,dnp)
! 877: N c,m,b;
! 878: int *sgnp;
! 879: N *nmp,*dnp;
! 880: {
! 881: Q qq,t,u1,v1,r1,nm;
! 882: N q,r,u2,v2,r2;
! 883:
! 884: u1 = 0; v1 = ONE; u2 = m; v2 = c;
! 885: while ( cmpn(v2,b) >= 0 ) {
! 886: divn(u2,v2,&q,&r2); u2 = v2; v2 = r2;
! 887: NTOQ(q,1,qq); mulq(qq,v1,&t); subq(u1,t,&r1); u1 = v1; v1 = r1;
! 888: }
! 889: if ( cmpn(NM(v1),b) >= 0 )
! 890: return 0;
! 891: else {
! 892: *nmp = v2;
! 893: *dnp = NM(v1);
! 894: *sgnp = SGN(v1);
! 895: return 1;
! 896: }
! 897: }
! 898:
! 899: /* mat->body = N ** */
! 900:
! 901: int intmtoratm(mat,md,nm,dn)
! 902: MAT mat;
! 903: N md;
! 904: MAT nm;
! 905: Q *dn;
! 906: {
! 907: N t,s,b;
! 908: Q bound,dn0,dn1,nm1,q,tq;
! 909: int i,j,k,l,row,col;
! 910: Q **rmat;
! 911: N **tmat;
! 912: N *tmi;
! 913: Q *nmk;
! 914: N u,unm,udn;
! 915: int sgn,ret;
! 916:
! 917: row = mat->row; col = mat->col;
! 918: bshiftn(md,1,&t);
! 919: isqrt(t,&s);
! 920: bshiftn(s,64,&b);
! 921: if ( !b )
! 922: b = ONEN;
! 923: dn0 = ONE;
! 924: tmat = (N **)mat->body;
! 925: rmat = (Q **)nm->body;
! 926: for ( i = 0; i < row; i++ )
! 927: for ( j = 0, tmi = tmat[i]; j < col; j++ )
! 928: if ( tmi[j] ) {
! 929: muln(tmi[j],NM(dn0),&s);
! 930: remn(s,md,&u);
! 931: ret = inttorat(u,md,b,&sgn,&unm,&udn);
! 932: if ( !ret )
! 933: return 0;
! 934: else {
! 935: NTOQ(unm,sgn,nm1);
! 936: NTOQ(udn,1,dn1);
! 937: if ( !UNIQ(dn1) ) {
! 938: for ( k = 0; k < i; k++ )
! 939: for ( l = 0, nmk = rmat[k]; l < col; l++ ) {
! 940: mulq(nmk[l],dn1,&q); nmk[l] = q;
! 941: }
! 942: for ( l = 0, nmk = rmat[i]; l < j; l++ ) {
! 943: mulq(nmk[l],dn1,&q); nmk[l] = q;
! 944: }
! 945: }
! 946: rmat[i][j] = nm1;
! 947: mulq(dn0,dn1,&q); dn0 = q;
! 948: }
! 949: }
! 950: *dn = dn0;
! 951: return 1;
! 952: }
! 953:
! 954: int generic_gauss_elim_mod(mat,row,col,md,colstat)
! 955: int **mat;
! 956: int row,col,md;
! 957: int *colstat;
! 958: {
! 959: int i,j,k,l,inv,a,rank;
! 960: int *t,*pivot;
! 961:
! 962: for ( rank = 0, j = 0; j < col; j++ ) {
! 963: for ( i = rank; i < row && !mat[i][j]; i++ );
! 964: if ( i == row ) {
! 965: colstat[j] = 0;
! 966: continue;
! 967: } else
! 968: colstat[j] = 1;
! 969: if ( i != rank ) {
! 970: t = mat[i]; mat[i] = mat[rank]; mat[rank] = t;
! 971: }
! 972: pivot = mat[rank];
! 973: inv = invm(pivot[j],md);
! 974: for ( k = j; k < col; k++ )
! 975: if ( pivot[k] ) {
! 976: DMAR(pivot[k],inv,0,md,pivot[k])
! 977: }
! 978: for ( i = rank+1; i < row; i++ ) {
! 979: t = mat[i];
! 980: if ( a = t[j] )
! 981: for ( k = j, a = md - a; k < col; k++ )
! 982: if ( pivot[k] ) {
! 983: DMAR(pivot[k],a,t[k],md,t[k])
! 984: }
! 985: }
! 986: rank++;
! 987: }
! 988: for ( j = col-1, l = rank-1; j >= 0; j-- )
! 989: if ( colstat[j] ) {
! 990: pivot = mat[l];
! 991: for ( i = 0; i < l; i++ ) {
! 992: t = mat[i];
! 993: if ( a = t[j] )
! 994: for ( k = j, a = md-a; k < col; k++ )
! 995: if ( pivot[k] ) {
! 996: DMAR(pivot[k],a,t[k],md,t[k])
! 997: }
! 998: }
! 999: l--;
! 1000: }
! 1001: return rank;
! 1002: }
! 1003:
! 1004: /* LU decomposition; a[i][i] = 1/U[i][i] */
! 1005:
! 1006: int lu_gfmmat(mat,md,perm)
! 1007: GFMMAT mat;
! 1008: unsigned int md;
! 1009: int *perm;
! 1010: {
! 1011: int row,col;
! 1012: int i,j,k,l;
! 1013: unsigned int *t,*pivot;
! 1014: unsigned int **a;
! 1015: unsigned int inv,m;
! 1016:
! 1017: row = mat->row; col = mat->col;
! 1018: a = mat->body;
! 1019: bzero(perm,row*sizeof(int));
! 1020:
! 1021: for ( i = 0; i < row; i++ )
! 1022: perm[i] = i;
! 1023: for ( k = 0; k < col; k++ ) {
! 1024: for ( i = k; i < row && !a[i][k]; i++ );
! 1025: if ( i == row )
! 1026: return 0;
! 1027: if ( i != k ) {
! 1028: j = perm[i]; perm[i] = perm[k]; perm[k] = j;
! 1029: t = a[i]; a[i] = a[k]; a[k] = t;
! 1030: }
! 1031: pivot = a[k];
! 1032: pivot[k] = inv = invm(pivot[k],md);
! 1033: for ( i = k+1; i < row; i++ ) {
! 1034: t = a[i];
! 1035: if ( m = t[k] ) {
! 1036: DMAR(inv,m,0,md,t[k])
! 1037: for ( j = k+1, m = md - t[k]; j < col; j++ )
! 1038: if ( pivot[j] ) {
! 1039: DMAR(m,pivot[j],t[j],md,t[j])
! 1040: }
! 1041: }
! 1042: }
! 1043: }
! 1044: return 1;
! 1045: }
! 1046:
! 1047: void Pleqm1(arg,rp)
! 1048: NODE arg;
! 1049: VECT *rp;
! 1050: {
! 1051: MAT m;
! 1052: VECT vect;
! 1053: pointer **mat;
! 1054: Q *v;
! 1055: Q q;
! 1056: int **wmat;
! 1057: int md,i,j,row,col,t,n,status;
! 1058:
! 1059: asir_assert(ARG0(arg),O_MAT,"leqm1");
! 1060: asir_assert(ARG1(arg),O_N,"leqm1");
! 1061: m = (MAT)ARG0(arg); md = QTOS((Q)ARG1(arg));
! 1062: row = m->row; col = m->col; mat = m->body;
! 1063: wmat = (int **)almat(row,col);
! 1064: for ( i = 0; i < row; i++ )
! 1065: for ( j = 0; j < col; j++ )
! 1066: if ( q = (Q)mat[i][j] ) {
! 1067: t = rem(NM(q),md);
! 1068: if ( SGN(q) < 0 )
! 1069: t = (md - t) % md;
! 1070: wmat[i][j] = t;
! 1071: } else
! 1072: wmat[i][j] = 0;
! 1073: status = gauss_elim_mod1(wmat,row,col,md);
! 1074: if ( status < 0 )
! 1075: *rp = 0;
! 1076: else if ( status > 0 )
! 1077: *rp = (VECT)ONE;
! 1078: else {
! 1079: n = col - 1;
! 1080: MKVECT(vect,n);
! 1081: for ( i = 0, v = (Q *)vect->body; i < n; i++ ) {
! 1082: t = (md-wmat[i][n])%md; STOQ(t,v[i]);
! 1083: }
! 1084: *rp = vect;
! 1085: }
! 1086: }
! 1087:
! 1088: gauss_elim_mod1(mat,row,col,md)
! 1089: int **mat;
! 1090: int row,col,md;
! 1091: {
! 1092: int i,j,k,inv,a,n;
! 1093: int *t,*pivot;
! 1094:
! 1095: n = col - 1;
! 1096: for ( j = 0; j < n; j++ ) {
! 1097: for ( i = j; i < row && !mat[i][j]; i++ );
! 1098: if ( i == row )
! 1099: return 1;
! 1100: if ( i != j ) {
! 1101: t = mat[i]; mat[i] = mat[j]; mat[j] = t;
! 1102: }
! 1103: pivot = mat[j];
! 1104: inv = invm(pivot[j],md);
! 1105: for ( k = j; k <= n; k++ )
! 1106: pivot[k] = dmar(pivot[k],inv,0,md);
! 1107: for ( i = j+1; i < row; i++ ) {
! 1108: t = mat[i];
! 1109: if ( i != j && (a = t[j]) )
! 1110: for ( k = j, a = md - a; k <= n; k++ )
! 1111: t[k] = dmar(pivot[k],a,t[k],md);
! 1112: }
! 1113: }
! 1114: for ( i = n; i < row && !mat[i][n]; i++ );
! 1115: if ( i == row ) {
! 1116: for ( j = n-1; j >= 0; j-- ) {
! 1117: for ( i = j-1, a = (md-mat[j][n])%md; i >= 0; i-- ) {
! 1118: mat[i][n] = dmar(mat[i][j],a,mat[i][n],md);
! 1119: mat[i][j] = 0;
! 1120: }
! 1121: }
! 1122: return 0;
! 1123: } else
! 1124: return -1;
! 1125: }
! 1126:
! 1127: void Pgeninvm(arg,rp)
! 1128: NODE arg;
! 1129: LIST *rp;
! 1130: {
! 1131: MAT m;
! 1132: pointer **mat;
! 1133: Q **tmat;
! 1134: Q q;
! 1135: unsigned int **wmat;
! 1136: int md,i,j,row,col,t,status;
! 1137: MAT mat1,mat2;
! 1138: NODE node1,node2;
! 1139:
! 1140: asir_assert(ARG0(arg),O_MAT,"leqm1");
! 1141: asir_assert(ARG1(arg),O_N,"leqm1");
! 1142: m = (MAT)ARG0(arg); md = QTOS((Q)ARG1(arg));
! 1143: row = m->row; col = m->col; mat = m->body;
! 1144: wmat = (unsigned int **)almat(row,col+row);
! 1145: for ( i = 0; i < row; i++ ) {
! 1146: bzero((char *)wmat[i],(col+row)*sizeof(int));
! 1147: for ( j = 0; j < col; j++ )
! 1148: if ( q = (Q)mat[i][j] ) {
! 1149: t = rem(NM(q),md);
! 1150: if ( SGN(q) < 0 )
! 1151: t = (md - t) % md;
! 1152: wmat[i][j] = t;
! 1153: }
! 1154: wmat[i][col+i] = 1;
! 1155: }
! 1156: status = gauss_elim_geninv_mod(wmat,row,col,md);
! 1157: if ( status > 0 )
! 1158: *rp = 0;
! 1159: else {
! 1160: MKMAT(mat1,col,row); MKMAT(mat2,row-col,row);
! 1161: for ( i = 0, tmat = (Q **)mat1->body; i < col; i++ )
! 1162: for ( j = 0; j < row; j++ )
! 1163: STOQ(wmat[i][j+col],tmat[i][j]);
! 1164: for ( tmat = (Q **)mat2->body; i < row; i++ )
! 1165: for ( j = 0; j < row; j++ )
! 1166: STOQ(wmat[i][j+col],tmat[i-col][j]);
! 1167: MKNODE(node2,mat2,0); MKNODE(node1,mat1,node2); MKLIST(*rp,node1);
! 1168: }
! 1169: }
! 1170:
! 1171: int gauss_elim_geninv_mod(mat,row,col,md)
! 1172: unsigned int **mat;
! 1173: int row,col,md;
! 1174: {
! 1175: int i,j,k,inv,a,n,m;
! 1176: unsigned int *t,*pivot;
! 1177:
! 1178: n = col; m = row+col;
! 1179: for ( j = 0; j < n; j++ ) {
! 1180: for ( i = j; i < row && !mat[i][j]; i++ );
! 1181: if ( i == row )
! 1182: return 1;
! 1183: if ( i != j ) {
! 1184: t = mat[i]; mat[i] = mat[j]; mat[j] = t;
! 1185: }
! 1186: pivot = mat[j];
! 1187: inv = invm(pivot[j],md);
! 1188: for ( k = j; k < m; k++ )
! 1189: pivot[k] = dmar(pivot[k],inv,0,md);
! 1190: for ( i = j+1; i < row; i++ ) {
! 1191: t = mat[i];
! 1192: if ( a = t[j] )
! 1193: for ( k = j, a = md - a; k < m; k++ )
! 1194: t[k] = dmar(pivot[k],a,t[k],md);
! 1195: }
! 1196: }
! 1197: for ( j = n-1; j >= 0; j-- ) {
! 1198: pivot = mat[j];
! 1199: for ( i = j-1; i >= 0; i-- ) {
! 1200: t = mat[i];
! 1201: if ( a = t[j] )
! 1202: for ( k = j, a = md - a; k < m; k++ )
! 1203: t[k] = dmar(pivot[k],a,t[k],md);
! 1204: }
! 1205: }
! 1206: return 0;
! 1207: }
! 1208:
! 1209: void Psolve_by_lu_gfmmat(arg,rp)
! 1210: NODE arg;
! 1211: VECT *rp;
! 1212: {
! 1213: GFMMAT lu;
! 1214: Q *perm,*rhs,*v;
! 1215: int n,i;
! 1216: unsigned int md;
! 1217: unsigned int *b,*sol;
! 1218: VECT r;
! 1219:
! 1220: lu = (GFMMAT)ARG0(arg);
! 1221: perm = (Q *)BDY((VECT)ARG1(arg));
! 1222: rhs = (Q *)BDY((VECT)ARG2(arg));
! 1223: md = (unsigned int)QTOS((Q)ARG3(arg));
! 1224: n = lu->col;
! 1225: b = (unsigned int *)MALLOC_ATOMIC(n*sizeof(int));
! 1226: sol = (unsigned int *)MALLOC_ATOMIC(n*sizeof(int));
! 1227: for ( i = 0; i < n; i++ )
! 1228: b[i] = QTOS(rhs[QTOS(perm[i])]);
! 1229: solve_by_lu_gfmmat(lu,md,b,sol);
! 1230: MKVECT(r,n);
! 1231: for ( i = 0, v = (Q *)r->body; i < n; i++ )
! 1232: STOQ(sol[i],v[i]);
! 1233: *rp = r;
! 1234: }
! 1235:
! 1236: void solve_by_lu_gfmmat(lu,md,b,x)
! 1237: GFMMAT lu;
! 1238: unsigned int md;
! 1239: unsigned int *b;
! 1240: unsigned int *x;
! 1241: {
! 1242: int n;
! 1243: unsigned int **a;
! 1244: unsigned int *y;
! 1245: int i,j;
! 1246: unsigned int t,m;
! 1247:
! 1248: n = lu->col;
! 1249: a = lu->body;
! 1250: y = (unsigned int *)MALLOC_ATOMIC(n*sizeof(int));
! 1251: /* solve Ly=b */
! 1252: for ( i = 0; i < n; i++ ) {
! 1253: for ( t = b[i], j = 0; j < i; j++ )
! 1254: if ( a[i][j] ) {
! 1255: m = md - a[i][j];
! 1256: DMAR(m,y[j],t,md,t)
! 1257: }
! 1258: y[i] = t;
! 1259: }
! 1260: /* solve Ux=y */
! 1261: for ( i = n-1; i >= 0; i-- ) {
! 1262: for ( t = y[i], j =i+1; j < n; j++ )
! 1263: if ( a[i][j] ) {
! 1264: m = md - a[i][j];
! 1265: DMAR(m,x[j],t,md,t)
! 1266: }
! 1267: /* a[i][i] = 1/U[i][i] */
! 1268: DMAR(t,a[i][i],0,md,x[i])
! 1269: }
! 1270: }
! 1271:
! 1272: void Plu_gfmmat(arg,rp)
! 1273: NODE arg;
! 1274: LIST *rp;
! 1275: {
! 1276: MAT m;
! 1277: GFMMAT mm;
! 1278: unsigned int md;
! 1279: int i,row,col,status;
! 1280: int *iperm;
! 1281: Q *v;
! 1282: VECT perm;
! 1283: NODE n0;
! 1284:
! 1285: asir_assert(ARG0(arg),O_MAT,"mat_to_gfmmat");
! 1286: asir_assert(ARG1(arg),O_N,"mat_to_gfmmat");
! 1287: m = (MAT)ARG0(arg); md = (unsigned int)QTOS((Q)ARG1(arg));
! 1288: mat_to_gfmmat(m,md,&mm);
! 1289: row = m->row;
! 1290: col = m->col;
! 1291: iperm = (int *)MALLOC_ATOMIC(row*sizeof(int));
! 1292: status = lu_gfmmat(mm,md,iperm);
! 1293: if ( !status )
! 1294: n0 = 0;
! 1295: else {
! 1296: MKVECT(perm,row);
! 1297: for ( i = 0, v = (Q *)perm->body; i < row; i++ )
! 1298: STOQ(iperm[i],v[i]);
! 1299: n0 = mknode(2,mm,perm);
! 1300: }
! 1301: MKLIST(*rp,n0);
! 1302: }
! 1303:
! 1304: void Pmat_to_gfmmat(arg,rp)
! 1305: NODE arg;
! 1306: GFMMAT *rp;
! 1307: {
! 1308: MAT m;
! 1309: unsigned int md;
! 1310:
! 1311: asir_assert(ARG0(arg),O_MAT,"mat_to_gfmmat");
! 1312: asir_assert(ARG1(arg),O_N,"mat_to_gfmmat");
! 1313: m = (MAT)ARG0(arg); md = (unsigned int)QTOS((Q)ARG1(arg));
! 1314: mat_to_gfmmat(m,md,rp);
! 1315: }
! 1316:
! 1317: void mat_to_gfmmat(m,md,rp)
! 1318: MAT m;
! 1319: unsigned int md;
! 1320: GFMMAT *rp;
! 1321: {
! 1322: unsigned int **wmat;
! 1323: unsigned int t;
! 1324: Q **mat;
! 1325: Q q;
! 1326: int i,j,row,col;
! 1327:
! 1328: row = m->row; col = m->col; mat = (Q **)m->body;
! 1329: wmat = (unsigned int **)almat(row,col);
! 1330: for ( i = 0; i < row; i++ ) {
! 1331: bzero((char *)wmat[i],col*sizeof(unsigned int));
! 1332: for ( j = 0; j < col; j++ )
! 1333: if ( q = mat[i][j] ) {
! 1334: t = (unsigned int)rem(NM(q),md);
! 1335: if ( SGN(q) < 0 )
! 1336: t = (md - t) % md;
! 1337: wmat[i][j] = t;
! 1338: }
! 1339: }
! 1340: TOGFMMAT(row,col,wmat,*rp);
! 1341: }
! 1342:
! 1343: void Pgeninvm_swap(arg,rp)
! 1344: NODE arg;
! 1345: LIST *rp;
! 1346: {
! 1347: MAT m;
! 1348: pointer **mat;
! 1349: Q **tmat;
! 1350: Q *tvect;
! 1351: Q q;
! 1352: unsigned int **wmat,**invmat;
! 1353: int *index;
! 1354: unsigned int t,md;
! 1355: int i,j,row,col,status;
! 1356: MAT mat1;
! 1357: VECT vect1;
! 1358: NODE node1,node2;
! 1359:
! 1360: asir_assert(ARG0(arg),O_MAT,"geninvm_swap");
! 1361: asir_assert(ARG1(arg),O_N,"geninvm_swap");
! 1362: m = (MAT)ARG0(arg); md = QTOS((Q)ARG1(arg));
! 1363: row = m->row; col = m->col; mat = m->body;
! 1364: wmat = (unsigned int **)almat(row,col+row);
! 1365: for ( i = 0; i < row; i++ ) {
! 1366: bzero((char *)wmat[i],(col+row)*sizeof(int));
! 1367: for ( j = 0; j < col; j++ )
! 1368: if ( q = (Q)mat[i][j] ) {
! 1369: t = (unsigned int)rem(NM(q),md);
! 1370: if ( SGN(q) < 0 )
! 1371: t = (md - t) % md;
! 1372: wmat[i][j] = t;
! 1373: }
! 1374: wmat[i][col+i] = 1;
! 1375: }
! 1376: status = gauss_elim_geninv_mod_swap(wmat,row,col,md,&invmat,&index);
! 1377: if ( status > 0 )
! 1378: *rp = 0;
! 1379: else {
! 1380: MKMAT(mat1,col,col);
! 1381: for ( i = 0, tmat = (Q **)mat1->body; i < col; i++ )
! 1382: for ( j = 0; j < col; j++ )
! 1383: UTOQ(invmat[i][j],tmat[i][j]);
! 1384: MKVECT(vect1,row);
! 1385: for ( i = 0, tvect = (Q *)vect1->body; i < row; i++ )
! 1386: STOQ(index[i],tvect[i]);
! 1387: MKNODE(node2,vect1,0); MKNODE(node1,mat1,node2); MKLIST(*rp,node1);
! 1388: }
! 1389: }
! 1390:
! 1391: gauss_elim_geninv_mod_swap(mat,row,col,md,invmatp,indexp)
! 1392: unsigned int **mat;
! 1393: int row,col;
! 1394: unsigned int md;
! 1395: unsigned int ***invmatp;
! 1396: int **indexp;
! 1397: {
! 1398: int i,j,k,inv,a,n,m;
! 1399: unsigned int *t,*pivot,*s;
! 1400: int *index;
! 1401: unsigned int **invmat;
! 1402:
! 1403: n = col; m = row+col;
! 1404: *indexp = index = (int *)MALLOC_ATOMIC(row*sizeof(int));
! 1405: for ( i = 0; i < row; i++ )
! 1406: index[i] = i;
! 1407: for ( j = 0; j < n; j++ ) {
! 1408: for ( i = j; i < row && !mat[i][j]; i++ );
! 1409: if ( i == row ) {
! 1410: *indexp = 0; *invmatp = 0; return 1;
! 1411: }
! 1412: if ( i != j ) {
! 1413: t = mat[i]; mat[i] = mat[j]; mat[j] = t;
! 1414: k = index[i]; index[i] = index[j]; index[j] = k;
! 1415: }
! 1416: pivot = mat[j];
! 1417: inv = (unsigned int)invm(pivot[j],md);
! 1418: for ( k = j; k < m; k++ )
! 1419: if ( pivot[k] )
! 1420: pivot[k] = (unsigned int)dmar(pivot[k],inv,0,md);
! 1421: for ( i = j+1; i < row; i++ ) {
! 1422: t = mat[i];
! 1423: if ( a = t[j] )
! 1424: for ( k = j, a = md - a; k < m; k++ )
! 1425: if ( pivot[k] )
! 1426: t[k] = dmar(pivot[k],a,t[k],md);
! 1427: }
! 1428: }
! 1429: for ( j = n-1; j >= 0; j-- ) {
! 1430: pivot = mat[j];
! 1431: for ( i = j-1; i >= 0; i-- ) {
! 1432: t = mat[i];
! 1433: if ( a = t[j] )
! 1434: for ( k = j, a = md - a; k < m; k++ )
! 1435: if ( pivot[k] )
! 1436: t[k] = dmar(pivot[k],a,t[k],md);
! 1437: }
! 1438: }
! 1439: *invmatp = invmat = (unsigned int **)almat(col,col);
! 1440: for ( i = 0; i < col; i++ )
! 1441: for ( j = 0, s = invmat[i], t = mat[i]; j < col; j++ )
! 1442: s[j] = t[col+index[j]];
! 1443: return 0;
! 1444: }
! 1445:
! 1446: void _addn(N,N,N);
! 1447: int _subn(N,N,N);
! 1448: void _muln(N,N,N);
! 1449:
! 1450: void inner_product_int(a,b,n,r)
! 1451: Q *a,*b;
! 1452: int n;
! 1453: Q *r;
! 1454: {
! 1455: int la,lb,i;
! 1456: int sgn,sgn1;
! 1457: N wm,wma,sum,t;
! 1458:
! 1459: for ( la = lb = 0, i = 0; i < n; i++ ) {
! 1460: if ( a[i] )
! 1461: if ( DN(a[i]) )
! 1462: error("inner_product_int : invalid argument");
! 1463: else
! 1464: la = MAX(PL(NM(a[i])),la);
! 1465: if ( b[i] )
! 1466: if ( DN(b[i]) )
! 1467: error("inner_product_int : invalid argument");
! 1468: else
! 1469: lb = MAX(PL(NM(b[i])),lb);
! 1470: }
! 1471: sgn = 0;
! 1472: sum= NALLOC(la+lb+2);
! 1473: bzero((char *)sum,(la+lb+3)*sizeof(unsigned int));
! 1474: wm = NALLOC(la+lb+2);
! 1475: wma = NALLOC(la+lb+2);
! 1476: for ( i = 0; i < n; i++ ) {
! 1477: if ( !a[i] || !b[i] )
! 1478: continue;
! 1479: _muln(NM(a[i]),NM(b[i]),wm);
! 1480: sgn1 = SGN(a[i])*SGN(b[i]);
! 1481: if ( !sgn ) {
! 1482: sgn = sgn1;
! 1483: t = wm; wm = sum; sum = t;
! 1484: } else if ( sgn == sgn1 ) {
! 1485: _addn(sum,wm,wma);
! 1486: if ( !PL(wma) )
! 1487: sgn = 0;
! 1488: t = wma; wma = sum; sum = t;
! 1489: } else {
! 1490: /* sgn*sum+sgn1*wm = sgn*(sum-wm) */
! 1491: sgn *= _subn(sum,wm,wma);
! 1492: t = wma; wma = sum; sum = t;
! 1493: }
! 1494: }
! 1495: GC_free(wm);
! 1496: GC_free(wma);
! 1497: if ( !sgn ) {
! 1498: GC_free(sum);
! 1499: *r = 0;
! 1500: } else
! 1501: NTOQ(sum,sgn,*r);
! 1502: }
! 1503:
! 1504: void Pmul_mat_vect_int(arg,rp)
! 1505: NODE arg;
! 1506: VECT *rp;
! 1507: {
! 1508: MAT mat;
! 1509: VECT vect,r;
! 1510: int row,col,i;
! 1511:
! 1512: mat = (MAT)ARG0(arg);
! 1513: vect = (VECT)ARG1(arg);
! 1514: row = mat->row;
! 1515: col = mat->col;
! 1516: MKVECT(r,row);
! 1517: for ( i = 0; i < row; i++ )
! 1518: inner_product_int(mat->body[i],vect->body,col,&r->body[i]);
! 1519: *rp = r;
! 1520: }
! 1521:
! 1522: void Pnbpoly_up2(arg,rp)
! 1523: NODE arg;
! 1524: GF2N *rp;
! 1525: {
! 1526: int m,type,ret;
! 1527: UP2 r;
! 1528:
! 1529: m = QTOS((Q)ARG0(arg));
! 1530: type = QTOS((Q)ARG1(arg));
! 1531: ret = generate_ONB_polynomial(&r,m,type);
! 1532: if ( ret == 0 )
! 1533: MKGF2N(r,*rp);
! 1534: else
! 1535: *rp = 0;
! 1536: }
! 1537:
! 1538: void Px962_irredpoly_up2(arg,rp)
! 1539: NODE arg;
! 1540: GF2N *rp;
! 1541: {
! 1542: int m,type,ret,w;
! 1543: GF2N prev;
! 1544: UP2 r;
! 1545:
! 1546: m = QTOS((Q)ARG0(arg));
! 1547: prev = (GF2N)ARG1(arg);
! 1548: if ( !prev ) {
! 1549: w = (m>>5)+1; NEWUP2(r,w); r->w = 0;
! 1550: bzero((char *)r->b,w*sizeof(unsigned int));
! 1551: } else {
! 1552: r = prev->body;
! 1553: if ( degup2(r) != m ) {
! 1554: w = (m>>5)+1; NEWUP2(r,w); r->w = 0;
! 1555: bzero((char *)r->b,w*sizeof(unsigned int));
! 1556: }
! 1557: }
! 1558: ret = _generate_irreducible_polynomial(r,m,type);
! 1559: if ( ret == 0 )
! 1560: MKGF2N(r,*rp);
! 1561: else
! 1562: *rp = 0;
! 1563: }
! 1564:
! 1565: void Pirredpoly_up2(arg,rp)
! 1566: NODE arg;
! 1567: GF2N *rp;
! 1568: {
! 1569: int m,type,ret,w;
! 1570: GF2N prev;
! 1571: UP2 r;
! 1572:
! 1573: m = QTOS((Q)ARG0(arg));
! 1574: prev = (GF2N)ARG1(arg);
! 1575: if ( !prev ) {
! 1576: w = (m>>5)+1; NEWUP2(r,w); r->w = 0;
! 1577: bzero((char *)r->b,w*sizeof(unsigned int));
! 1578: } else {
! 1579: r = prev->body;
! 1580: if ( degup2(r) != m ) {
! 1581: w = (m>>5)+1; NEWUP2(r,w); r->w = 0;
! 1582: bzero((char *)r->b,w*sizeof(unsigned int));
! 1583: }
! 1584: }
! 1585: ret = _generate_good_irreducible_polynomial(r,m,type);
! 1586: if ( ret == 0 )
! 1587: MKGF2N(r,*rp);
! 1588: else
! 1589: *rp = 0;
! 1590: }
! 1591:
! 1592: /*
! 1593: * f = type 'type' normal polynomial of degree m if exists
! 1594: * IEEE P1363 A.7.2
! 1595: *
! 1596: * return value : 0 --- exists
! 1597: * 1 --- does not exist
! 1598: * -1 --- failure (memory allocation error)
! 1599: */
! 1600:
! 1601: int generate_ONB_polynomial(UP2 *rp,int m,int type)
! 1602: {
! 1603: int i,r;
! 1604: int w;
! 1605: UP2 f,f0,f1,f2,t;
! 1606:
! 1607: w = (m>>5)+1;
! 1608: switch ( type ) {
! 1609: case 1:
! 1610: if ( !TypeT_NB_check(m,1) ) return 1;
! 1611: NEWUP2(f,w); *rp = f; f->w = w;
! 1612: /* set all the bits */
! 1613: for ( i = 0; i < w; i++ )
! 1614: f->b[i] = 0xffffffff;
! 1615: /* mask the top word if necessary */
! 1616: if ( r = (m+1)&31 )
! 1617: f->b[w-1] &= (1<<r)-1;
! 1618: return 0;
! 1619: break;
! 1620: case 2:
! 1621: if ( !TypeT_NB_check(m,2) ) return 1;
! 1622: NEWUP2(f,w); *rp = f;
! 1623: W_NEWUP2(f0,w);
! 1624: W_NEWUP2(f1,w);
! 1625: W_NEWUP2(f2,w);
! 1626:
! 1627: /* recursion for genrating Type II normal polynomial */
! 1628:
! 1629: /* f0 = 1, f1 = t+1 */
! 1630: f0->w = 1; f0->b[0] = 1;
! 1631: f1->w = 1; f1->b[0] = 3;
! 1632: for ( i = 2; i <= m; i++ ) {
! 1633: /* f2 = t*f1+f0 */
! 1634: _bshiftup2(f1,-1,f2);
! 1635: _addup2_destructive(f2,f0);
! 1636: /* cyclic change of the variables */
! 1637: t = f0; f0 = f1; f1 = f2; f2 = t;
! 1638: }
! 1639: _copyup2(f1,f);
! 1640: return 0;
! 1641: break;
! 1642: default:
! 1643: return -1;
! 1644: break;
! 1645: }
! 1646: }
! 1647:
! 1648: /*
! 1649: * f = an irreducible trinomial or pentanomial of degree d 'after' f
! 1650: * return value : 0 --- exists
! 1651: * 1 --- does not exist (exhaustion)
! 1652: */
! 1653:
! 1654: int _generate_irreducible_polynomial(UP2 f,int d)
! 1655: {
! 1656: int ret,i,j,k,nz,i0,j0,k0;
! 1657: int w;
! 1658: unsigned int *fd;
! 1659:
! 1660: /*
! 1661: * if f = x^d+x^i+1 then i0 <- i, j0 <- 0, k0 <-0.
! 1662: * if f = x^d+x^k+x^j+x^i+1 (k>j>i) then i0 <- i, j0 <- j, k0 <-k.
! 1663: * otherwise i0,j0,k0 is set to 0.
! 1664: */
! 1665:
! 1666: fd = f->b;
! 1667: w = (d>>5)+1;
! 1668: if ( f->w && (d==degup2(f)) ) {
! 1669: for ( nz = 0, i = d; i >= 0; i-- )
! 1670: if ( fd[i>>5]&(1<<(i&31)) ) nz++;
! 1671: switch ( nz ) {
! 1672: case 3:
! 1673: for ( i0 = 1; !(fd[i0>>5]&(1<<(i0&31))) ; i0++ );
! 1674: /* reset i0-th bit */
! 1675: fd[i0>>5] &= ~(1<<(i0&31));
! 1676: j0 = k0 = 0;
! 1677: break;
! 1678: case 5:
! 1679: for ( i0 = 1; !(fd[i0>>5]&(1<<(i0&31))) ; i0++ );
! 1680: /* reset i0-th bit */
! 1681: fd[i0>>5] &= ~(1<<(i0&31));
! 1682: for ( j0 = i0+1; !(fd[j0>>5]&(1<<(j0&31))) ; j0++ );
! 1683: /* reset j0-th bit */
! 1684: fd[j0>>5] &= ~(1<<(j0&31));
! 1685: for ( k0 = j0+1; !(fd[k0>>5]&(1<<(k0&31))) ; k0++ );
! 1686: /* reset k0-th bit */
! 1687: fd[k0>>5] &= ~(1<<(k0&31));
! 1688: break;
! 1689: default:
! 1690: f->w = 0; break;
! 1691: }
! 1692: } else
! 1693: f->w = 0;
! 1694:
! 1695: if ( !f->w ) {
! 1696: fd = f->b;
! 1697: f->w = w; fd[0] |= 1; fd[d>>5] |= (1<<(d&31));
! 1698: i0 = j0 = k0 = 0;
! 1699: }
! 1700: /* if j0 > 0 then f is already a pentanomial */
! 1701: if ( j0 > 0 ) goto PENTA;
! 1702:
! 1703: /* searching for an irreducible trinomial */
! 1704:
! 1705: for ( i = 1; 2*i <= d; i++ ) {
! 1706: /* skip the polynomials 'before' f */
! 1707: if ( i < i0 ) continue;
! 1708: if ( i == i0 ) { i0 = 0; continue; }
! 1709: /* set i-th bit */
! 1710: fd[i>>5] |= (1<<(i&31));
! 1711: ret = irredcheck_dddup2(f);
! 1712: if ( ret == 1 ) return 0;
! 1713: /* reset i-th bit */
! 1714: fd[i>>5] &= ~(1<<(i&31));
! 1715: }
! 1716:
! 1717: /* searching for an irreducible pentanomial */
! 1718: PENTA:
! 1719: for ( i = 1; i < d; i++ ) {
! 1720: /* skip the polynomials 'before' f */
! 1721: if ( i < i0 ) continue;
! 1722: if ( i == i0 ) i0 = 0;
! 1723: /* set i-th bit */
! 1724: fd[i>>5] |= (1<<(i&31));
! 1725: for ( j = i+1; j < d; j++ ) {
! 1726: /* skip the polynomials 'before' f */
! 1727: if ( j < j0 ) continue;
! 1728: if ( j == j0 ) j0 = 0;
! 1729: /* set j-th bit */
! 1730: fd[j>>5] |= (1<<(j&31));
! 1731: for ( k = j+1; k < d; k++ ) {
! 1732: /* skip the polynomials 'before' f */
! 1733: if ( k < k0 ) continue;
! 1734: else if ( k == k0 ) { k0 = 0; continue; }
! 1735: /* set k-th bit */
! 1736: fd[k>>5] |= (1<<(k&31));
! 1737: ret = irredcheck_dddup2(f);
! 1738: if ( ret == 1 ) return 0;
! 1739: /* reset k-th bit */
! 1740: fd[k>>5] &= ~(1<<(k&31));
! 1741: }
! 1742: /* reset j-th bit */
! 1743: fd[j>>5] &= ~(1<<(j&31));
! 1744: }
! 1745: /* reset i-th bit */
! 1746: fd[i>>5] &= ~(1<<(i&31));
! 1747: }
! 1748: /* exhausted */
! 1749: return 1;
! 1750: }
! 1751:
! 1752: /*
! 1753: * f = an irreducible trinomial or pentanomial of degree d 'after' f
! 1754: *
! 1755: * searching strategy:
! 1756: * trinomial x^d+x^i+1:
! 1757: * i is as small as possible.
! 1758: * trinomial x^d+x^i+x^j+x^k+1:
! 1759: * i is as small as possible.
! 1760: * For such i, j is as small as possible.
! 1761: * For such i and j, 'k' is as small as possible.
! 1762: *
! 1763: * return value : 0 --- exists
! 1764: * 1 --- does not exist (exhaustion)
! 1765: */
! 1766:
! 1767: int _generate_good_irreducible_polynomial(UP2 f,int d)
! 1768: {
! 1769: int ret,i,j,k,nz,i0,j0,k0;
! 1770: int w;
! 1771: unsigned int *fd;
! 1772:
! 1773: /*
! 1774: * if f = x^d+x^i+1 then i0 <- i, j0 <- 0, k0 <-0.
! 1775: * if f = x^d+x^k+x^j+x^i+1 (k>j>i) then i0 <- i, j0 <- j, k0 <-k.
! 1776: * otherwise i0,j0,k0 is set to 0.
! 1777: */
! 1778:
! 1779: fd = f->b;
! 1780: w = (d>>5)+1;
! 1781: if ( f->w && (d==degup2(f)) ) {
! 1782: for ( nz = 0, i = d; i >= 0; i-- )
! 1783: if ( fd[i>>5]&(1<<(i&31)) ) nz++;
! 1784: switch ( nz ) {
! 1785: case 3:
! 1786: for ( i0 = 1; !(fd[i0>>5]&(1<<(i0&31))) ; i0++ );
! 1787: /* reset i0-th bit */
! 1788: fd[i0>>5] &= ~(1<<(i0&31));
! 1789: j0 = k0 = 0;
! 1790: break;
! 1791: case 5:
! 1792: for ( i0 = 1; !(fd[i0>>5]&(1<<(i0&31))) ; i0++ );
! 1793: /* reset i0-th bit */
! 1794: fd[i0>>5] &= ~(1<<(i0&31));
! 1795: for ( j0 = i0+1; !(fd[j0>>5]&(1<<(j0&31))) ; j0++ );
! 1796: /* reset j0-th bit */
! 1797: fd[j0>>5] &= ~(1<<(j0&31));
! 1798: for ( k0 = j0+1; !(fd[k0>>5]&(1<<(k0&31))) ; k0++ );
! 1799: /* reset k0-th bit */
! 1800: fd[k0>>5] &= ~(1<<(k0&31));
! 1801: break;
! 1802: default:
! 1803: f->w = 0; break;
! 1804: }
! 1805: } else
! 1806: f->w = 0;
! 1807:
! 1808: if ( !f->w ) {
! 1809: fd = f->b;
! 1810: f->w = w; fd[0] |= 1; fd[d>>5] |= (1<<(d&31));
! 1811: i0 = j0 = k0 = 0;
! 1812: }
! 1813: /* if j0 > 0 then f is already a pentanomial */
! 1814: if ( j0 > 0 ) goto PENTA;
! 1815:
! 1816: /* searching for an irreducible trinomial */
! 1817:
! 1818: for ( i = 1; 2*i <= d; i++ ) {
! 1819: /* skip the polynomials 'before' f */
! 1820: if ( i < i0 ) continue;
! 1821: if ( i == i0 ) { i0 = 0; continue; }
! 1822: /* set i-th bit */
! 1823: fd[i>>5] |= (1<<(i&31));
! 1824: ret = irredcheck_dddup2(f);
! 1825: if ( ret == 1 ) return 0;
! 1826: /* reset i-th bit */
! 1827: fd[i>>5] &= ~(1<<(i&31));
! 1828: }
! 1829:
! 1830: /* searching for an irreducible pentanomial */
! 1831: PENTA:
! 1832: for ( i = 3; i < d; i++ ) {
! 1833: /* skip the polynomials 'before' f */
! 1834: if ( i < i0 ) continue;
! 1835: if ( i == i0 ) i0 = 0;
! 1836: /* set i-th bit */
! 1837: fd[i>>5] |= (1<<(i&31));
! 1838: for ( j = 2; j < i; j++ ) {
! 1839: /* skip the polynomials 'before' f */
! 1840: if ( j < j0 ) continue;
! 1841: if ( j == j0 ) j0 = 0;
! 1842: /* set j-th bit */
! 1843: fd[j>>5] |= (1<<(j&31));
! 1844: for ( k = 1; k < j; k++ ) {
! 1845: /* skip the polynomials 'before' f */
! 1846: if ( k < k0 ) continue;
! 1847: else if ( k == k0 ) { k0 = 0; continue; }
! 1848: /* set k-th bit */
! 1849: fd[k>>5] |= (1<<(k&31));
! 1850: ret = irredcheck_dddup2(f);
! 1851: if ( ret == 1 ) return 0;
! 1852: /* reset k-th bit */
! 1853: fd[k>>5] &= ~(1<<(k&31));
! 1854: }
! 1855: /* reset j-th bit */
! 1856: fd[j>>5] &= ~(1<<(j&31));
! 1857: }
! 1858: /* reset i-th bit */
! 1859: fd[i>>5] &= ~(1<<(i&31));
! 1860: }
! 1861: /* exhausted */
! 1862: return 1;
! 1863: }