Annotation of OpenXM_contrib2/asir2018/builtin/gf.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 "parse.h"
! 52:
! 53: struct resf_dlist {
! 54: int ib,id;
! 55: };
! 56:
! 57: int resf_degtest(int,int *,int,struct resf_dlist *);
! 58: void uhensel(P,NODE,int,int,NODE *);
! 59: void uhensel_incremental(P,NODE,int,int,int,NODE *);
! 60: void resf_hensel(int,P,int,P *,ML *);
! 61: void resf_dtest(P,ML,int,int *,int *,DCP *);
! 62: void resf_dtest_special(P,ML,int,int *,int *,DCP *);
! 63: void dtest_special(P,ML,P *);
! 64: void hensel_special(int,int,P,P *,ML *);
! 65:
! 66: void nullspace(UM **,UM,int,int,int *);
! 67: void nullspace_lm(LM **,int,int *);
! 68: void nullspace_gf2n(GF2N **,int,int *);
! 69: void nullspace_gfpn(GFPN **,int,int *);
! 70: void nullspace_gfs(GFS **,int,int *);
! 71: void nullspace_gfsn(GFSN **,int,int *);
! 72: void null_to_sol(int **,int *,int,int,UM *);
! 73:
! 74: void showgfmat(UM **,int);
! 75: void pwr_mod(P,P,V,P,int,Z,P *);
! 76: void rem_mod(P,P,V,P,int,P *);
! 77:
! 78: void Pnullspace(),Pgcda_mod(),Pftest(),Presfmain(),Ppwr_mod(),Puhensel();
! 79: void Puhensel_incremental();
! 80: void Psfuhensel();
! 81:
! 82: void Pnullspace_ff();
! 83:
! 84: void Psolve_linear_equation_gf2n();
! 85: void Plinear_form_to_vect(),Pvect_to_linear_form();
! 86:
! 87: void solve_linear_equation_gf2n(GF2N **,int,int,int *);
! 88: void linear_form_to_array(P,VL,int,Num *);
! 89: void array_to_linear_form(Num *,VL,int,P *);
! 90: void sfuhensel(P,NODE,GFS,int,NODE *);
! 91:
! 92: extern int current_ff;
! 93:
! 94: struct ftab gf_tab[] = {
! 95: {"solve_linear_equation_gf2n",Psolve_linear_equation_gf2n,1},
! 96: {"nullspace",Pnullspace,3},
! 97: {"nullspace_ff",Pnullspace_ff,1},
! 98: /* {"gcda_mod",Pgcda_mod,5}, */
! 99: {"ftest",Pftest,2},
! 100: {"resfmain",Presfmain,4},
! 101: {"pwr_mod",Ppwr_mod,6},
! 102: {"uhensel",Puhensel,4},
! 103: {"uhensel_incremental",Puhensel_incremental,5},
! 104: {"sfuhensel",Psfuhensel,4},
! 105: {0,0,0},
! 106: };
! 107:
! 108: int resf_degtest(int k,int *in,int nb,struct resf_dlist *dlist)
! 109: {
! 110: int i,d0,d;
! 111: int dl_i;
! 112: struct resf_dlist *t;
! 113:
! 114: if ( k < nb )
! 115: return 0;
! 116: if ( nb == 1 )
! 117: return 1;
! 118: d0 = 0; d = 0; dl_i = 0;
! 119: for ( i = 0; i < k; i++ ) {
! 120: t = &dlist[in[i]];
! 121: if ( t->ib > dl_i + 1 )
! 122: return 0;
! 123: else if ( t->ib == dl_i )
! 124: d += t->id;
! 125: else if ( !d || (dl_i && d0 != d) )
! 126: return 0;
! 127: else {
! 128: d0 = d; dl_i++; d = t->id;
! 129: }
! 130: }
! 131: if ( dl_i != nb - 1 || d0 != d )
! 132: return 0;
! 133: else
! 134: return 1;
! 135: }
! 136:
! 137: void Puhensel(NODE arg,LIST *rp)
! 138: {
! 139: P f;
! 140: NODE mfl,r;
! 141: int mod,bound;
! 142:
! 143: f = (P)ARG0(arg);
! 144: mfl = BDY((LIST)ARG1(arg));
! 145: mod = QTOS((Q)ARG2(arg));
! 146: bound = QTOS((Q)ARG3(arg));
! 147: uhensel(f,mfl,mod,bound,&r);
! 148: MKLIST(*rp,r);
! 149: }
! 150:
! 151: void Puhensel_incremental(NODE arg,LIST *rp)
! 152: {
! 153: P f;
! 154: NODE mfl,r;
! 155: int mod,bound,start;
! 156:
! 157: f = (P)ARG0(arg);
! 158: mfl = BDY((LIST)ARG1(arg));
! 159: mod = QTOS((Q)ARG2(arg));
! 160: start = QTOS((Q)ARG3(arg));
! 161: bound = QTOS((Q)ARG4(arg));
! 162: uhensel_incremental(f,mfl,mod,start,bound,&r);
! 163: MKLIST(*rp,r);
! 164: }
! 165:
! 166: void uhensel(P f,NODE mfl,int mod,int bound,NODE *rp)
! 167: {
! 168: ML blist,clist,rlist;
! 169: LUM fl;
! 170: int nf,i;
! 171: P s;
! 172: V v;
! 173: NODE t,top;
! 174:
! 175: nf = length(mfl);
! 176: blist = MLALLOC(nf); blist->n = nf; blist->mod = mod;
! 177: for ( i = 0, t = mfl; i < nf; i++, t = NEXT(t) ) {
! 178: blist->c[i] = (pointer)UMALLOC(UDEG((P)BDY(t)));
! 179: ptoum(mod,(P)BDY(t),blist->c[i]);
! 180: }
! 181: gcdgen(f,blist,&clist);
! 182: blist->bound = clist->bound = bound;
! 183: W_LUMALLOC((int)UDEG(f),bound,fl);
! 184: ptolum(mod,bound,f,fl);
! 185: henmain(fl,blist,clist,&rlist);
! 186: v = VR(f);
! 187: for ( i = nf-1, top = 0; i >= 0; i-- ) {
! 188: lumtop(v,mod,bound,rlist->c[i],&s);
! 189: MKNODE(t,s,top); top = t;
! 190: }
! 191: *rp = top;
! 192: }
! 193:
! 194: void uhensel_incremental(P f,NODE mfl,int mod,int start,int bound,NODE *rp)
! 195: {
! 196: ML blist,clist,rlist;
! 197: LUM fl;
! 198: LUM *lblist;
! 199: int nf,i,j,k;
! 200: int **p;
! 201: P s;
! 202: V v;
! 203: NODE t,top;
! 204:
! 205: nf = length(mfl);
! 206: blist = MLALLOC(nf); blist->n = nf; blist->mod = mod;
! 207: lblist = (LUM *)MALLOC(nf*sizeof(LUM));
! 208: for ( i = 0, t = mfl; i < nf; i++, t = NEXT(t) ) {
! 209: blist->c[i] = (pointer)UMALLOC(UDEG((P)BDY(t)));
! 210: ptoum(mod,(P)BDY(t),blist->c[i]);
! 211: W_LUMALLOC((int)UDEG((P)BDY(t)),bound,lblist[i]);
! 212: ptolum(mod,start,(P)BDY(t),lblist[i]);
! 213: p = lblist[i]->c;
! 214: for ( j = DEG(lblist[i]); j >= 0; j-- )
! 215: for ( k = start; k < bound; k++ )
! 216: p[j][k] = 0;
! 217: }
! 218: gcdgen(f,blist,&clist);
! 219: clist->bound = bound;
! 220: W_LUMALLOC((int)UDEG(f),bound,fl);
! 221: ptolum(mod,bound,f,fl);
! 222: henmain_incremental(fl,lblist,clist,nf,mod,start,bound);
! 223: v = VR(f);
! 224: for ( i = nf-1, top = 0; i >= 0; i-- ) {
! 225: lumtop_unsigned(v,mod,bound,lblist[i],&s);
! 226: MKNODE(t,s,top); top = t;
! 227: }
! 228: *rp = top;
! 229: }
! 230:
! 231: void Psfuhensel(NODE arg,LIST *rp)
! 232: {
! 233: P f;
! 234: int bound;
! 235: NODE r,mfl;
! 236: GFS ev;
! 237:
! 238: f = (P)ARG0(arg);
! 239: mfl = BDY((LIST)ARG1(arg));
! 240: ev = (GFS)ARG2(arg);
! 241: bound = QTOS((Q)ARG3(arg));
! 242: sfuhensel(f,mfl,ev,bound,&r);
! 243: MKLIST(*rp,r);
! 244: }
! 245:
! 246: void sfuhensel(P f,NODE mfl,GFS ev,int bound,NODE *rp)
! 247: {
! 248: BM fl;
! 249: BM *r;
! 250: VL vl,nvl;
! 251: int i,fn,dx,dy,d;
! 252: NODE t,top;
! 253: UM fm,hm,q;
! 254: UM *gm;
! 255: V x,y;
! 256: P g,s,u;
! 257:
! 258: clctv(CO,f,&vl);
! 259: if ( !vl || !vl->next || vl->next->next )
! 260: error("sfuhensel : f must be a bivariate poly");
! 261:
! 262: for ( i = 0, t = mfl; t; i++, t = NEXT(t) );
! 263: fn = i;
! 264:
! 265: gm = (UM *)MALLOC(fn*sizeof(UM));
! 266:
! 267: /* XXX : more severe check is necessary */
! 268: x = VR((P)BDY(mfl));
! 269: y = vl->v == x ? vl->next->v : vl->v;
! 270:
! 271: for ( i = 0, t = mfl, d = 0; i < fn; i++, t = NEXT(t) ) {
! 272: gm[i] = (pointer)UMALLOC(getdeg(x,(P)BDY(t)));
! 273: ptosfum((P)BDY(t),gm[i]);
! 274: d += DEG(gm[i]);
! 275: }
! 276:
! 277: /* reorder f if necessary */
! 278: if ( vl->v != x ) {
! 279: reordvar(vl,x,&nvl); reorderp(nvl,vl,f,&g);
! 280: vl = nvl; f = g;
! 281: }
! 282: dx = getdeg(x,f);
! 283: if ( dx != d )
! 284: error("sfuhensel : product of factors has incompatible degree");
! 285:
! 286: dy = getdeg(y,f);
! 287: dy = MAX(dy,bound);
! 288: fl = BMALLOC(dx,dy);
! 289: ptosfbm(dy,f,fl);
! 290: if ( ev ) shiftsfbm(fl,FTOIF(CONT(ev)));
! 291:
! 292: /* fm = fl mod y */
! 293: fm = W_UMALLOC(dx);
! 294: cpyum(COEF(fl)[0],fm);
! 295: hm = W_UMALLOC(dx);
! 296:
! 297: q = W_UMALLOC(dx);
! 298: r = (BM *)MLALLOC(fn*sizeof(BM));
! 299: for ( i = 0; i < fn-1; i++ ) {
! 300: /* fl = gm[i]*hm mod y */
! 301: divsfum(fm,gm[i],hm);
! 302: /* fl is replaced by the cofactor of gk mod y^bound */
! 303: /* r[i] = gk */
! 304: sfhenmain2(fl,gm[i],hm,bound,r+i);
! 305: cpyum(hm,fm);
! 306: }
! 307: /* finally, fl must be the lift of gm[fn-1] */
! 308: r[i] = fl;
! 309:
! 310: for ( i = fn-1, top = 0; i >= 0; i-- ) {
! 311: sfbmtop(r[i],x,y,&s);
! 312: reorderp(CO,vl,s,&u);
! 313: MKNODE(t,u,top); top = t;
! 314: }
! 315: *rp = top;
! 316: }
! 317:
! 318: void Presfmain(NODE arg,LIST *rp)
! 319: {
! 320: P f;
! 321: int mod,n,nb,i,j,k;
! 322: int *nf,*md;
! 323: P *mf;
! 324: NODE mfl,mdl,t,s,u;
! 325: ML list;
! 326: DCP dc;
! 327: int sflag;
! 328:
! 329: f = (P)ARG0(arg);
! 330: mod = QTOS((Q)ARG1(arg));
! 331: mfl = BDY((LIST)ARG2(arg));
! 332: mdl = BDY((LIST)ARG3(arg));
! 333: for ( n = nb = 0, t = mfl; t; nb++, t = NEXT(t) )
! 334: for ( s = BDY((LIST)BDY(t)); s; n++, s = NEXT(s) );
! 335: if ( n == nb ) {
! 336: /* f must be irreducible! */
! 337: NEWDC(dc);
! 338: dc->c = f; dc->d = ONE;
! 339: } else {
! 340: nf = W_ALLOC(nb); md = W_ALLOC(nb); mf = (P *)ALLOCA(n*sizeof(P));
! 341: for ( k = i = 0, t = mfl, u = mdl, sflag = 1; t;
! 342: i++, t = NEXT(t), u = NEXT(u) ) {
! 343: if ( sflag && length(BDY((LIST)BDY(t))) != 2 )
! 344: sflag = 0;
! 345: for ( j = 0, s = BDY((LIST)BDY(t)); s; j++, s = NEXT(s) )
! 346: mf[k++] = (P)BDY(s);
! 347: nf[i] = j; md[i] = QTOS((Q)BDY(u));
! 348: }
! 349: resf_hensel(mod,f,n,mf,&list);
! 350: if ( sflag )
! 351: resf_dtest_special(f,list,nb,nf,md,&dc);
! 352: else
! 353: resf_dtest(f,list,nb,nf,md,&dc);
! 354: }
! 355: dcptolist(dc,rp);
! 356: }
! 357:
! 358: void resf_hensel(int mod,P f,int nf,P *mfl,ML *listp)
! 359: {
! 360: register int i,j;
! 361: int q,n,bound,inv,lc;
! 362: int *p;
! 363: int **pp;
! 364: ML blist,clist,bqlist,cqlist,rlist;
! 365: UM *b;
! 366: LUM fl,tl;
! 367: LUM *l;
! 368: UM w;
! 369:
! 370: w = W_UMALLOC(UDEG(f));
! 371: blist = MLALLOC(nf); blist->n = nf; blist->mod = mod;
! 372:
! 373: /* c[0] must have lc(f) */
! 374: blist->c[0] = (pointer)UMALLOC(UDEG(mfl[0]));
! 375: ptoum(mod,mfl[0],w);
! 376: inv = invm(w->c[UDEG(mfl[0])],mod);
! 377: lc = remqi((Q)LC(f),mod);
! 378: lc = dmar(inv,lc,0,mod);
! 379: if ( lc == 1 )
! 380: copyum(w,blist->c[0]);
! 381: else
! 382: mulsum(mod,w,lc,blist->c[0]);
! 383:
! 384: /* c[i] (i=1,...,nf-1) must be monic */
! 385: for ( i = 1; i < nf; i++ ) {
! 386: blist->c[i] = (pointer)UMALLOC(UDEG(mfl[i]));
! 387: ptoum(mod,mfl[i],w);
! 388: inv = invm(w->c[UDEG(mfl[i])],mod);
! 389: if ( inv == 1 )
! 390: copyum(w,blist->c[i]);
! 391: else
! 392: mulsum(mod,w,inv,blist->c[i]);
! 393: }
! 394:
! 395: gcdgen(f,blist,&clist); henprep(f,blist,clist,&bqlist,&cqlist);
! 396: n = bqlist->n; q = bqlist->mod;
! 397: bqlist->bound = cqlist->bound = bound = mignotte(q,f);
! 398: if ( bound == 1 ) {
! 399: *listp = rlist = MLALLOC(n);
! 400: rlist->n = n; rlist->mod = q; rlist->bound = bound;
! 401: for ( i = 0, b = (UM *)bqlist->c, l = (LUM *)rlist->c; i < n; i++ ) {
! 402: tl = LUMALLOC(DEG(b[i]),1); l[i] = tl; p = COEF(b[i]);
! 403: for ( j = 0, pp = COEF(tl); j <= DEG(tl); j++ )
! 404: pp[j][0] = p[j];
! 405: }
! 406: } else {
! 407: W_LUMALLOC((int)UDEG(f),bound,fl);
! 408: ptolum(q,bound,f,fl); henmain(fl,bqlist,cqlist,listp);
! 409: }
! 410: }
! 411:
! 412: void resf_dtest(P f,ML list,int nb,int *nfl,int *mdl,DCP *dcp)
! 413: {
! 414: int n,np,bound,q;
! 415: int i,j,k;
! 416: int *win;
! 417: P g,factor,cofactor;
! 418: Q csum,csumt;
! 419: DCP dcf,dcf0;
! 420: LUM *c;
! 421: ML wlist;
! 422: struct resf_dlist *dlist;
! 423: int z;
! 424:
! 425: n = UDEG(f); np = list->n; bound = list->bound; q = list->mod;
! 426: win = W_ALLOC(np+1);
! 427: ucsump(f,&csum); mulq(csum,(Q)COEF(DC(f)),&csumt); csum = csumt;
! 428: wlist = W_MLALLOC(np); wlist->n = list->n;
! 429: wlist->mod = list->mod; wlist->bound = list->bound;
! 430: c = (LUM *)COEF(wlist);
! 431: bcopy((char *)COEF(list),(char *)c,(int)(sizeof(LUM)*np));
! 432: dlist = (struct resf_dlist *)ALLOCA(np*sizeof(struct resf_dlist));
! 433: for ( i = 0, j = 0; j < nb; j++ )
! 434: for ( k = 0; k < nfl[j]; k++, i++ ) {
! 435: dlist[i].ib = j; dlist[i].id = DEG(c[i])/mdl[j];
! 436: }
! 437: k = nb;
! 438: for ( i = 0; i < nb; i++ )
! 439: win[i] = i+1;
! 440: for ( g = f, dcf = dcf0 = 0, --np, z = 0; ; ) {
! 441: #if 0
! 442: if ( !(z++ % 10000) )
! 443: fputc('.',stderr);
! 444: #endif
! 445: if ( resf_degtest(k,win,nb,dlist) &&
! 446: dtestmain(g,csum,wlist,k,win,&factor,&cofactor) ) {
! 447: NEXTDC(dcf0,dcf); DEG(dcf) = ONE; COEF(dcf) = factor;
! 448: g = cofactor;
! 449: ucsump(g,&csum); mulq(csum,(Q)COEF(DC(g)),&csumt); csum = csumt;
! 450: for ( i = 0; i < k - 1; i++ )
! 451: for ( j = win[i] + 1; j < win[i + 1]; j++ ) {
! 452: c[j-i-1] = c[j];
! 453: dlist[j-i-1] = dlist[j];
! 454: }
! 455: for ( j = win[k-1] + 1; j <= np; j++ ) {
! 456: c[j-k] = c[j];
! 457: dlist[j-k] = dlist[j];
! 458: }
! 459: if ( ( np -= k ) < k )
! 460: break;
! 461: if ( np - win[0] + 1 < k )
! 462: if ( ++k > np )
! 463: break;
! 464: else
! 465: for ( i = 0; i < k; i++ )
! 466: win[i] = i + 1;
! 467: else
! 468: for ( i = 1; i < k; i++ )
! 469: win[i] = win[0] + i;
! 470: } else if ( !ncombi(1,np,k,win) )
! 471: if ( k == np )
! 472: break;
! 473: else
! 474: for ( i = 0, ++k; i < k; i++ )
! 475: win[i] = i + 1;
! 476: }
! 477: NEXTDC(dcf0,dcf); COEF(dcf) = g;
! 478: DEG(dcf) = ONE; NEXT(dcf) = 0;*dcp = dcf0;
! 479: }
! 480:
! 481: void resf_dtest_special(P f,ML list,int nb,int *nfl,int *mdl,DCP *dcp)
! 482: {
! 483: int n,np,bound,q;
! 484: int i,j;
! 485: int *win;
! 486: P g,factor,cofactor;
! 487: Q csum,csumt;
! 488: DCP dcf,dcf0;
! 489: LUM *c;
! 490: ML wlist;
! 491: int max;
! 492:
! 493: n = UDEG(f); np = list->n; bound = list->bound; q = list->mod;
! 494: win = W_ALLOC(np+1);
! 495: ucsump(f,&csum); mulq(csum,(Q)COEF(DC(f)),&csumt); csum = csumt;
! 496: wlist = W_MLALLOC(np); wlist->n = list->n;
! 497: wlist->mod = list->mod; wlist->bound = list->bound;
! 498: c = (LUM *)COEF(wlist); bcopy((char *)COEF(list),(char *)c,(int)(sizeof(LUM)*np));
! 499: max = 1<<nb;
! 500: for ( g = f, dcf = dcf0 = 0, i = 0; i < max; i++ ) {
! 501: #if 0
! 502: if ( !(i % 1000) )
! 503: fprintf(stderr,"i=%d\n",i);
! 504: #endif
! 505: for ( j = 0; j < nb; j++ )
! 506: win[j] = (i&(1<<j)) ? 2*j+1 : 2*j;
! 507: if ( dtestmain(g,csum,wlist,nb,win,&factor,&cofactor) ) {
! 508: #if 0
! 509: fprintf(stderr,"success : i=%d\n",i);
! 510: #endif
! 511: NEXTDC(dcf0,dcf); DEG(dcf) = ONE; COEF(dcf) = factor;
! 512: NEXTDC(dcf0,dcf); DEG(dcf) = ONE; COEF(dcf) = cofactor;
! 513: NEXT(dcf) = 0;*dcp = dcf0;
! 514: return;
! 515: }
! 516: }
! 517: NEXTDC(dcf0,dcf); COEF(dcf) = g;
! 518: DEG(dcf) = ONE; NEXT(dcf) = 0;*dcp = dcf0;
! 519: }
! 520:
! 521: #if 0
! 522: void Pftest(arg,rp)
! 523: NODE arg;
! 524: P *rp;
! 525: {
! 526: ML list;
! 527: DCP dc;
! 528: P p;
! 529: P *mfl;
! 530:
! 531: p = (P)ARG0(arg); mfl = (P *)(((VECT)ARG1(arg))->body);
! 532: hensel_special(4,1,p,mfl,&list);
! 533: dtest_special(p,list,rp);
! 534: }
! 535:
! 536: void dtest_special(f,list,pr)
! 537: P f;
! 538: ML list;
! 539: P *pr;
! 540: {
! 541: int n,np,bound,q;
! 542: int i,j,k;
! 543: int *win;
! 544: P g,factor,cofactor;
! 545: Q csum,csumt;
! 546: DCP dc,dcf,dcf0;
! 547: LUM *c;
! 548: ML wlist;
! 549:
! 550: n = UDEG(f); np = list->n; bound = list->bound; q = list->mod;
! 551: win = W_ALLOC(np+1);
! 552: ucsump(f,&csum); mulq(csum,(Q)COEF(DC(f)),&csumt); csum = csumt;
! 553: wlist = W_MLALLOC(np); wlist->n = list->n;
! 554: wlist->mod = list->mod; wlist->bound = list->bound;
! 555: c = (LUM *)COEF(wlist); bcopy((char *)COEF(list),(char *)c,(int)(sizeof(LUM)*np));
! 556: for ( g = f, i = 130000; i < (1<<20); i++ ) {
! 557: #if 0
! 558: if ( !(i % 1000) )
! 559: fprintf(stderr,"i=%d\n",i);
! 560: #endif
! 561: for ( j = 0; j < 20; j++ )
! 562: win[j] = (i&(1<<j)) ? 2*j+1 : 2*j;
! 563: if ( dtestmain(g,csum,wlist,20,win,&factor,&cofactor) ) {
! 564: #if 0
! 565: fprintf(stderr,"success : i=%d\n",i);
! 566: #endif
! 567: *pr = factor; return;
! 568: }
! 569: }
! 570: }
! 571:
! 572: void hensel_special(index,count,f,mfl,listp)
! 573: int index,count;
! 574: P f;
! 575: P *mfl;
! 576: ML *listp;
! 577: {
! 578: register int i,j;
! 579: int q,n,t,d,r,u,br,tmp,bound;
! 580: int *c,*p,*m,*w;
! 581: int **pp;
! 582: DCP dc;
! 583: ML blist,clist,bqlist,cqlist,rlist;
! 584: UM *b;
! 585: LUM fl,tl;
! 586: LUM *l;
! 587:
! 588: blist = MLALLOC(40); blist->n = 40; blist->mod = 11;
! 589: for ( i = 0; i < 40; i++ ) {
! 590: blist->c[i] = (pointer)UMALLOC(6);
! 591: ptoum(11,mfl[i],blist->c[i]);
! 592: }
! 593: gcdgen(f,blist,&clist); henprep(f,blist,clist,&bqlist,&cqlist);
! 594: n = bqlist->n; q = bqlist->mod;
! 595: bqlist->bound = cqlist->bound = bound = mignotte(q,f);
! 596: if ( bound == 1 ) {
! 597: *listp = rlist = MLALLOC(n);
! 598: rlist->n = n; rlist->mod = q; rlist->bound = bound;
! 599: for ( i = 0, b = (UM *)bqlist->c, l = (LUM *)rlist->c; i < n; i++ ) {
! 600: tl = LUMALLOC(DEG(b[i]),1); l[i] = tl; p = COEF(b[i]);
! 601: for ( j = 0, pp = COEF(tl); j <= DEG(tl); j++ )
! 602: pp[j][0] = p[j];
! 603: }
! 604: } else {
! 605: W_LUMALLOC(UDEG(f),bound,fl);
! 606: ptolum(q,bound,f,fl); henmain(fl,bqlist,cqlist,listp);
! 607: }
! 608: }
! 609: #endif
! 610:
! 611: #if 0
! 612: void Pftest(arg,rp)
! 613: NODE arg;
! 614: P *rp;
! 615: {
! 616: ML list;
! 617: DCP dc;
! 618: P p;
! 619: P *mfl;
! 620:
! 621: p = (P)ARG0(arg); mfl = (P *)(((VECT)ARG1(arg))->body);
! 622: hensel_special(2,1,p,mfl,&list);
! 623: dtest_special(p,list,rp);
! 624: }
! 625:
! 626: void dtest_special(f,list,pr)
! 627: P f;
! 628: ML list;
! 629: P *pr;
! 630: {
! 631: int n,np,bound,q;
! 632: int i,j,k,t,b0;
! 633: int *win;
! 634: P g,factor,cofactor;
! 635: Q csum,csumt;
! 636: DCP dc,dcf,dcf0;
! 637: LUM *c;
! 638: ML wlist;
! 639: static int nbits[16] = {0,1,1,2,1,2,2,3,1,2,2,3,2,3,3,4};
! 640:
! 641: n = UDEG(f); np = list->n; bound = list->bound; q = list->mod;
! 642: win = W_ALLOC(np+1);
! 643: ucsump(f,&csum); mulq(csum,(Q)COEF(DC(f)),&csumt); csum = csumt;
! 644: wlist = W_MLALLOC(np); wlist->n = list->n;
! 645: wlist->mod = list->mod; wlist->bound = list->bound;
! 646: c = (LUM *)COEF(wlist); bcopy((char *)COEF(list),(char *)c,(int)(sizeof(LUM)*np));
! 647: for ( g = f, i = 0; i < (1<<23); i++ ) {
! 648: #if 0
! 649: if ( !(i % 1000) )
! 650: fprintf(stderr,"i=%d\n",i);
! 651: #endif
! 652: t = i>>20; b0 = nbits[t];
! 653: if ( !b0 )
! 654: continue;
! 655: for ( j = 1; j < 6; j++ ) {
! 656: t = (i>>(20-4*j))&0xf;
! 657: if ( nbits[t] != b0 )
! 658: break;
! 659: }
! 660: if ( j != 6 )
! 661: continue;
! 662: for ( j = k = 0; j < 24; j++ )
! 663: if ( i & (1<<(23-j)) )
! 664: win[k++] = j;
! 665: if ( dtestmain(g,csum,wlist,k,win,&factor,&cofactor) ) {
! 666: #if 0
! 667: fprintf(stderr,"success : i=%d\n",i);
! 668: #endif
! 669: *pr = factor; return;
! 670: }
! 671: }
! 672: *pr = f;
! 673: }
! 674:
! 675: void hensel_special(index,count,f,mfl,listp)
! 676: int index,count;
! 677: P f;
! 678: P *mfl;
! 679: ML *listp;
! 680: {
! 681: register int i,j;
! 682: int q,n,t,d,r,u,br,tmp,bound;
! 683: int *c,*p,*m,*w;
! 684: int **pp;
! 685: DCP dc;
! 686: ML blist,clist,bqlist,cqlist,rlist;
! 687: UM *b;
! 688: LUM fl,tl;
! 689: LUM *l;
! 690:
! 691: blist = MLALLOC(24); blist->n = 24; blist->mod = 5;
! 692: for ( i = 0; i < 24; i++ ) {
! 693: blist->c[i] = (pointer)UMALLOC(7);
! 694: ptoum(5,mfl[i],blist->c[i]);
! 695: }
! 696: gcdgen(f,blist,&clist); henprep(f,blist,clist,&bqlist,&cqlist);
! 697: n = bqlist->n; q = bqlist->mod;
! 698: bqlist->bound = cqlist->bound = bound = mignotte(q,f);
! 699: if ( bound == 1 ) {
! 700: *listp = rlist = MLALLOC(n);
! 701: rlist->n = n; rlist->mod = q; rlist->bound = bound;
! 702: for ( i = 0, b = (UM *)bqlist->c, l = (LUM *)rlist->c; i < n; i++ ) {
! 703: tl = LUMALLOC(DEG(b[i]),1); l[i] = tl; p = COEF(b[i]);
! 704: for ( j = 0, pp = COEF(tl); j <= DEG(tl); j++ )
! 705: pp[j][0] = p[j];
! 706: }
! 707: } else {
! 708: W_LUMALLOC(UDEG(f),bound,fl);
! 709: ptolum(q,bound,f,fl); henmain(fl,bqlist,cqlist,listp);
! 710: }
! 711: }
! 712: #endif
! 713:
! 714: void Pftest(NODE arg,P *rp)
! 715: {
! 716: ML list;
! 717: P p;
! 718: P *mfl;
! 719:
! 720: p = (P)ARG0(arg); mfl = (P *)(((VECT)ARG1(arg))->body);
! 721: hensel_special(5,1,p,mfl,&list);
! 722: dtest_special(p,list,rp);
! 723: }
! 724:
! 725: int nbits(int a)
! 726: {
! 727: int i,s;
! 728:
! 729: for ( i = 0, s = 0; a && (i < 20); i++, a >>= 1 )
! 730: if ( a & 1 ) s++;
! 731: return s;
! 732: }
! 733:
! 734: void dtest_special(P f,ML list,P *pr)
! 735: {
! 736: int n,np,bound,q;
! 737: int i,j,k,b0;
! 738: int *win;
! 739: P g,factor,cofactor;
! 740: Q csum,csumt;
! 741: LUM *c;
! 742: ML wlist;
! 743:
! 744: n = UDEG(f); np = list->n; bound = list->bound; q = list->mod;
! 745: win = W_ALLOC(np+1);
! 746: ucsump(f,&csum); mulq(csum,(Q)COEF(DC(f)),&csumt); csum = csumt;
! 747: wlist = W_MLALLOC(np); wlist->n = list->n;
! 748: wlist->mod = list->mod; wlist->bound = list->bound;
! 749: c = (LUM *)COEF(wlist); bcopy((char *)COEF(list),(char *)c,(int)(sizeof(LUM)*np));
! 750: for ( g = f, i = 0; i < (1<<19); i++ ) {
! 751: #if 0
! 752: if ( !(i % 10000) )
! 753: fprintf(stderr,"i=%d\n",i);
! 754: #endif
! 755: b0 = nbits(i>>10);
! 756: if ( !b0 || (nbits(i&0x3ff) != b0) )
! 757: continue;
! 758: for ( j = k = 0; j < 20; j++ )
! 759: if ( i & (1<<(19-j)) )
! 760: win[k++] = j;
! 761: if ( dtestmain(g,csum,wlist,k,win,&factor,&cofactor) ) {
! 762: #if 0
! 763: fprintf(stderr,"success : i=%d\n",i);
! 764: #endif
! 765: *pr = factor; return;
! 766: }
! 767: }
! 768: *pr = f;
! 769: }
! 770:
! 771: void hensel_special(int index,int count,P f,P *mfl,ML *listp)
! 772: {
! 773: register int i,j;
! 774: int q,n,bound;
! 775: int *p;
! 776: int **pp;
! 777: ML blist,clist,bqlist,cqlist,rlist;
! 778: UM *b;
! 779: LUM fl,tl;
! 780: LUM *l;
! 781:
! 782: blist = MLALLOC(20); blist->n = 20; blist->mod = 11;
! 783: for ( i = 0; i < 20; i++ ) {
! 784: blist->c[i] = (pointer)UMALLOC(10);
! 785: ptoum(11,mfl[i],blist->c[i]);
! 786: }
! 787: gcdgen(f,blist,&clist); henprep(f,blist,clist,&bqlist,&cqlist);
! 788: n = bqlist->n; q = bqlist->mod;
! 789: bqlist->bound = cqlist->bound = bound = mignotte(q,f);
! 790: if ( bound == 1 ) {
! 791: *listp = rlist = MLALLOC(n);
! 792: rlist->n = n; rlist->mod = q; rlist->bound = bound;
! 793: for ( i = 0, b = (UM *)bqlist->c, l = (LUM *)rlist->c; i < n; i++ ) {
! 794: tl = LUMALLOC(DEG(b[i]),1); l[i] = tl; p = COEF(b[i]);
! 795: for ( j = 0, pp = COEF(tl); j <= DEG(tl); j++ )
! 796: pp[j][0] = p[j];
! 797: }
! 798: } else {
! 799: W_LUMALLOC((int)UDEG(f),bound,fl);
! 800: ptolum(q,bound,f,fl); henmain(fl,bqlist,cqlist,listp);
! 801: }
! 802: }
! 803:
! 804: void Pnullspace(NODE arg,LIST *rp)
! 805: {
! 806: int i,j,n,mod;
! 807: MAT mat,r;
! 808: VECT u;
! 809: V v;
! 810: P p,z;
! 811: Z q;
! 812: UM **w;
! 813: UM mp;
! 814: P *t;
! 815: UM *s;
! 816: int *ind;
! 817: NODE n0,n1;
! 818:
! 819: mat = (MAT)ARG0(arg);
! 820: p = (P)ARG1(arg);
! 821: v = VR(p);
! 822: mod = QTOS((Q)ARG2(arg));
! 823: n = mat->row;
! 824: w = (UM **)almat_pointer(n,n);
! 825: for ( i = 0; i < n; i++ )
! 826: for ( j = 0, t = (P *)mat->body[i], s = w[i]; j < n; j++ ) {
! 827: ptomp(mod,t[j],&z);
! 828: s[j] = W_UMALLOC((z&&!NUM(z))?UDEG(z):0);
! 829: mptoum(z,s[j]);
! 830: }
! 831: mp = W_UMALLOC(UDEG(p)); ptoum(mod,p,mp);
! 832: ind = (int *)ALLOCA(n*sizeof(int));
! 833: nullspace(w,mp,mod,n,ind);
! 834: MKMAT(r,n,n);
! 835: for ( i = 0; i < n; i++ )
! 836: for ( j = 0, t = (P *)r->body[i], s = w[i]; j < n; j++ )
! 837: umtop(v,s[j],&t[j]);
! 838: MKVECT(u,n);
! 839: for ( i = 0; i < n; i++ ) {
! 840: STOQ(ind[i],q); u->body[i] = (pointer)q;
! 841: }
! 842: MKNODE(n1,u,0); MKNODE(n0,r,n1); MKLIST(*rp,n0);
! 843: }
! 844:
! 845: void nullspace(UM **mat,UM p,int mod,int n,int *ind)
! 846: {
! 847: int i,j,l,s,d;
! 848: UM q,w,w1,h,inv;
! 849: UM *t,*u;
! 850:
! 851: d = DEG(p); inv = W_UMALLOC(d); q = W_UMALLOC(2*d);
! 852: w = W_UMALLOC(2*d); w1 = W_UMALLOC(2*d); h = W_UMALLOC(d);
! 853: bzero(ind,n*sizeof(int));
! 854: ind[0] = 0;
! 855: for ( i = j = 0; j < n; i++, j++ ) {
! 856: for ( ; j < n; j++ ) {
! 857: for ( l = i; l < n; l++ )
! 858: if ( DEG(mat[l][j])>=0 )
! 859: break;
! 860: if ( l < n ) {
! 861: t = mat[i]; mat[i] = mat[l]; mat[l] = t; break;
! 862: } else
! 863: ind[j] = 1;
! 864: }
! 865: if ( j == n )
! 866: break;
! 867: invum(mod,p,mat[i][j],inv);
! 868: for ( s = j, t = mat[i]; s < n; s++ ) {
! 869: mulum(mod,t[s],inv,w);
! 870: DEG(w) = divum(mod,w,p,q);
! 871: cpyum(w,t[s]);
! 872: }
! 873: for ( l = 0; l < n; l++ ) {
! 874: if ( l == i )
! 875: continue;
! 876: u = mat[l]; DEG(w) = -1; subum(mod,w,u[j],h);
! 877: for ( s = j; s < n; s++ ) {
! 878: mulum(mod,h,t[s],w); addum(mod,w,u[s],w1);
! 879: DEG(w1) = divum(mod,w1,p,q); cpyum(w1,u[s]);
! 880: }
! 881: }
! 882: }
! 883: }
! 884:
! 885: void Pnullspace_ff(NODE arg,LIST *rp)
! 886: {
! 887: int i,j,n;
! 888: MAT mat,r;
! 889: VECT u;
! 890: Z q;
! 891: Obj **w;
! 892: Obj *t;
! 893: Obj *s;
! 894: int *ind;
! 895: NODE n0,n1;
! 896:
! 897: mat = (MAT)ARG0(arg);
! 898: n = mat->row;
! 899: w = (Obj **)almat_pointer(n,n);
! 900: for ( i = 0; i < n; i++ )
! 901: for ( j = 0, t = (Obj *)mat->body[i], s = w[i]; j < n; j++ )
! 902: s[j] = t[j];
! 903: ind = (int *)ALLOCA(n*sizeof(int));
! 904: switch ( current_ff ) {
! 905: case FF_GFP:
! 906: nullspace_lm((LM **)w,n,ind); break;
! 907: case FF_GF2N:
! 908: nullspace_gf2n((GF2N **)w,n,ind); break;
! 909: case FF_GFPN:
! 910: nullspace_gfpn((GFPN **)w,n,ind); break;
! 911: case FF_GFS:
! 912: nullspace_gfs((GFS **)w,n,ind); break;
! 913: case FF_GFSN:
! 914: nullspace_gfsn((GFSN **)w,n,ind); break;
! 915: default:
! 916: error("nullspace_ff : current_ff is not set");
! 917: }
! 918: MKMAT(r,n,n);
! 919: for ( i = 0; i < n; i++ )
! 920: for ( j = 0, t = (Obj *)r->body[i], s = w[i]; j < n; j++ )
! 921: t[j] = s[j];
! 922: MKVECT(u,n);
! 923: for ( i = 0; i < n; i++ ) {
! 924: STOQ(ind[i],q); u->body[i] = (pointer)q;
! 925: }
! 926: MKNODE(n1,u,0); MKNODE(n0,r,n1); MKLIST(*rp,n0);
! 927: }
! 928:
! 929: void nullspace_lm(LM **mat,int n,int *ind)
! 930: {
! 931: int i,j,l,s;
! 932: Z q,z,mod;
! 933: LM w,w1,h,inv;
! 934: LM *t,*u;
! 935:
! 936: getmod_lm(&mod);
! 937:
! 938: bzero(ind,n*sizeof(int));
! 939: ind[0] = 0;
! 940: for ( i = j = 0; j < n; i++, j++ ) {
! 941: for ( ; j < n; j++ ) {
! 942: for ( l = i; l < n; l++ ) {
! 943: simplm(mat[l][j],&w); mat[l][j] = w;
! 944: if ( mat[l][j] )
! 945: break;
! 946: }
! 947: if ( l < n ) {
! 948: t = mat[i]; mat[i] = mat[l]; mat[l] = t; break;
! 949: } else
! 950: ind[j] = 1;
! 951: }
! 952: if ( j == n )
! 953: break;
! 954: lmtolf(mat[i][j],&q); invz(q,mod,&z); MKLM(BDY(z),inv);
! 955: for ( s = j, t = mat[i]; s < n; s++ ) {
! 956: mullm(t[s],inv,&w); t[s] = w;
! 957: }
! 958: for ( l = 0; l < n; l++ ) {
! 959: if ( l == i )
! 960: continue;
! 961: u = mat[l]; chsgnlm(u[j],&h);
! 962: for ( s = j; s < n; s++ ) {
! 963: mullm(h,t[s],&w); addlm(w,u[s],&w1); u[s] = w1;
! 964: }
! 965: }
! 966: }
! 967: }
! 968:
! 969: void nullspace_gf2n(GF2N **mat,int n,int *ind)
! 970: {
! 971: int i,j,l,s;
! 972: GF2N w,w1,h,inv;
! 973: GF2N *t,*u;
! 974: extern int gf2n_lazy;
! 975:
! 976: bzero(ind,n*sizeof(int));
! 977: ind[0] = 0;
! 978: for ( i = j = 0; j < n; i++, j++ ) {
! 979: for ( ; j < n; j++ ) {
! 980: for ( l = i; l < n; l++ ) {
! 981: simpgf2n(mat[l][j],&w); mat[l][j] = w;
! 982: if ( mat[l][j] )
! 983: break;
! 984: }
! 985: if ( l < n ) {
! 986: t = mat[i]; mat[i] = mat[l]; mat[l] = t; break;
! 987: } else
! 988: ind[j] = 1;
! 989: }
! 990: if ( j == n )
! 991: break;
! 992: invgf2n(mat[i][j],&inv);
! 993: for ( s = j, t = mat[i]; s < n; s++ ) {
! 994: mulgf2n(t[s],inv,&w); t[s] = w;
! 995: }
! 996: for ( l = 0; l < n; l++ ) {
! 997: if ( l == i )
! 998: continue;
! 999: u = mat[l]; h = u[j];
! 1000: for ( s = j; s < n; s++ ) {
! 1001: mulgf2n(h,t[s],&w); addgf2n(w,u[s],&w1); u[s] = w1;
! 1002: }
! 1003: }
! 1004: }
! 1005: }
! 1006:
! 1007: void nullspace_gfpn(GFPN **mat,int n,int *ind)
! 1008: {
! 1009: int i,j,l,s;
! 1010: GFPN w,w1,h,inv;
! 1011: GFPN *t,*u;
! 1012:
! 1013: bzero(ind,n*sizeof(int));
! 1014: ind[0] = 0;
! 1015: for ( i = j = 0; j < n; i++, j++ ) {
! 1016: for ( ; j < n; j++ ) {
! 1017: for ( l = i; l < n; l++ ) {
! 1018: simpgfpn(mat[l][j],&w); mat[l][j] = w;
! 1019: if ( mat[l][j] )
! 1020: break;
! 1021: }
! 1022: if ( l < n ) {
! 1023: t = mat[i]; mat[i] = mat[l]; mat[l] = t; break;
! 1024: } else
! 1025: ind[j] = 1;
! 1026: }
! 1027: if ( j == n )
! 1028: break;
! 1029: divgfpn((GFPN)ONE,(GFPN)mat[i][j],&inv);
! 1030: for ( s = j, t = mat[i]; s < n; s++ ) {
! 1031: mulgfpn(t[s],inv,&w); t[s] = w;
! 1032: }
! 1033: for ( l = 0; l < n; l++ ) {
! 1034: if ( l == i )
! 1035: continue;
! 1036: u = mat[l]; chsgngfpn(u[j],&h);
! 1037: for ( s = j; s < n; s++ ) {
! 1038: mulgfpn(h,t[s],&w); addgfpn(w,u[s],&w1); u[s] = w1;
! 1039: }
! 1040: }
! 1041: }
! 1042: }
! 1043:
! 1044: void nullspace_gfs(GFS **mat,int n,int *ind)
! 1045: {
! 1046: int i,j,l,s;
! 1047: GFS w,w1,h,inv;
! 1048: GFS *t,*u;
! 1049: GFS one;
! 1050:
! 1051: bzero(ind,n*sizeof(int));
! 1052: ind[0] = 0;
! 1053: mqtogfs(ONEM,&one);
! 1054:
! 1055: for ( i = j = 0; j < n; i++, j++ ) {
! 1056: for ( ; j < n; j++ ) {
! 1057: for ( l = i; l < n; l++ )
! 1058: if ( mat[l][j] )
! 1059: break;
! 1060: if ( l < n ) {
! 1061: t = mat[i]; mat[i] = mat[l]; mat[l] = t; break;
! 1062: } else
! 1063: ind[j] = 1;
! 1064: }
! 1065: if ( j == n )
! 1066: break;
! 1067: divgfs(one,mat[i][j],&inv);
! 1068: for ( s = j, t = mat[i]; s < n; s++ ) {
! 1069: mulgfs(t[s],inv,&w); t[s] = w;
! 1070: }
! 1071: for ( l = 0; l < n; l++ ) {
! 1072: if ( l == i )
! 1073: continue;
! 1074: u = mat[l];
! 1075: chsgngfs(u[j],&h);
! 1076: for ( s = j; s < n; s++ ) {
! 1077: mulgfs(h,t[s],&w); addgfs(w,u[s],&w1); u[s] = w1;
! 1078: }
! 1079: }
! 1080: }
! 1081: }
! 1082:
! 1083: void nullspace_gfsn(GFSN **mat,int n,int *ind)
! 1084: {
! 1085: int i,j,l,s;
! 1086: GFSN w,w1,h,inv;
! 1087: GFSN *t,*u;
! 1088:
! 1089: bzero(ind,n*sizeof(int));
! 1090: ind[0] = 0;
! 1091:
! 1092: for ( i = j = 0; j < n; i++, j++ ) {
! 1093: for ( ; j < n; j++ ) {
! 1094: for ( l = i; l < n; l++ )
! 1095: if ( mat[l][j] )
! 1096: break;
! 1097: if ( l < n ) {
! 1098: t = mat[i]; mat[i] = mat[l]; mat[l] = t; break;
! 1099: } else
! 1100: ind[j] = 1;
! 1101: }
! 1102: if ( j == n )
! 1103: break;
! 1104: invgfsn(mat[i][j],&inv);
! 1105: for ( s = j, t = mat[i]; s < n; s++ ) {
! 1106: mulgfsn(t[s],inv,&w); t[s] = w;
! 1107: }
! 1108: for ( l = 0; l < n; l++ ) {
! 1109: if ( l == i )
! 1110: continue;
! 1111: u = mat[l];
! 1112: chsgngfsn(u[j],&h);
! 1113: for ( s = j; s < n; s++ ) {
! 1114: mulgfsn(h,t[s],&w); addgfsn(w,u[s],&w1); u[s] = w1;
! 1115: }
! 1116: }
! 1117: }
! 1118: }
! 1119:
! 1120: /* p = a(0)vl[0]+a(1)vl[1]+...+a(m-1)vl[m-1]+a(m) -> array = [a(0) a(1) ... a(m)] */
! 1121:
! 1122: void linear_form_to_array(P p,VL vl,int m,Num *array)
! 1123: {
! 1124: int i;
! 1125: DCP dc;
! 1126:
! 1127: bzero((char *)array,(m+1)*sizeof(Num *));
! 1128: for ( i = 0; p && vl; vl = NEXT(vl), i++ ) {
! 1129: if ( ID(p) == O_N )
! 1130: break;
! 1131: else if ( VR(p) == vl->v ) {
! 1132: dc = DC(p);
! 1133: array[i] = (Num)COEF(dc);
! 1134: dc = NEXT(dc);
! 1135: p = dc ? COEF(dc) : 0;
! 1136: }
! 1137: }
! 1138: array[m] = (Num)p;
! 1139: }
! 1140:
! 1141: void array_to_linear_form(Num *array,VL vl,int m,P *r)
! 1142: {
! 1143: P t;
! 1144: DCP dc0,dc1;
! 1145:
! 1146: if ( !m )
! 1147: *r = (P)array[0];
! 1148: else {
! 1149: array_to_linear_form(array+1,NEXT(vl),m-1,&t);
! 1150: if ( !array[0] )
! 1151: *r = t;
! 1152: else {
! 1153: NEWDC(dc0); DEG(dc0) = ONE; COEF(dc0) = (P)array[0];
! 1154: if ( !t )
! 1155: NEXT(dc0) = 0;
! 1156: else {
! 1157: NEWDC(dc1); DEG(dc1) = 0; COEF(dc1) = t;
! 1158: NEXT(dc1) = 0;
! 1159: NEXT(dc0) = dc1;
! 1160: }
! 1161: MKP(vl->v,dc0,*r);
! 1162: }
! 1163: }
! 1164: }
! 1165:
! 1166: void Psolve_linear_equation_gf2n(NODE arg,LIST *rp)
! 1167: {
! 1168: NODE eqs,tn;
! 1169: VL vars,tvl;
! 1170: int i,j,n,m,dim,codim;
! 1171: GF2N **w;
! 1172: int *ind;
! 1173: NODE n0,n1;
! 1174:
! 1175: get_vars(ARG0(arg),&vars);
! 1176: eqs = BDY((LIST)ARG0(arg));
! 1177: for ( n = 0, tn = eqs; tn; tn = NEXT(tn), n++);
! 1178: for ( m = 0, tvl = vars; tvl; tvl = NEXT(tvl), m++);
! 1179: w = (GF2N **)almat_pointer(n,m+1);
! 1180: for ( i = 0, tn = eqs; i < n; i++, tn = NEXT(tn) )
! 1181: linear_form_to_array(BDY(tn),vars,m,(Num *)w[i]);
! 1182: ind = (int *)ALLOCA(m*sizeof(int));
! 1183: solve_linear_equation_gf2n(w,n,m,ind);
! 1184: for ( j = 0, dim = 0; j < m; j++ )
! 1185: if ( ind[j] )
! 1186: dim++;
! 1187: codim = m-dim;
! 1188: for ( i = codim; i < n; i++ )
! 1189: if ( w[i][m] ) {
! 1190: MKLIST(*rp,0); return;
! 1191: }
! 1192: for ( i = 0, n0 = 0; i < codim; i++ ) {
! 1193: NEXTNODE(n0,n1);
! 1194: array_to_linear_form((Num *)w[i],vars,m,(P *)&BDY(n1));
! 1195: }
! 1196: if ( n0 )
! 1197: NEXT(n1) = 0;
! 1198: MKLIST(*rp,n0);
! 1199: }
! 1200:
! 1201: void solve_linear_equation_gf2n(GF2N **mat,int n,int m,int *ind)
! 1202: {
! 1203: int i,j,l,s;
! 1204: GF2N w,w1,h,inv;
! 1205: GF2N *t,*u;
! 1206: extern int gf2n_lazy;
! 1207:
! 1208: bzero(ind,m*sizeof(int));
! 1209: ind[0] = 0;
! 1210: for ( i = j = 0; j < m; i++, j++ ) {
! 1211: for ( ; j < m; j++ ) {
! 1212: for ( l = i; l < n; l++ ) {
! 1213: simpgf2n(mat[l][j],&w); mat[l][j] = w;
! 1214: if ( mat[l][j] )
! 1215: break;
! 1216: }
! 1217: if ( l < n ) {
! 1218: t = mat[i]; mat[i] = mat[l]; mat[l] = t; break;
! 1219: } else
! 1220: ind[j] = 1;
! 1221: }
! 1222: if ( j == m )
! 1223: break;
! 1224: invgf2n(mat[i][j],&inv);
! 1225: for ( s = j, t = mat[i]; s <= m; s++ ) {
! 1226: mulgf2n(t[s],inv,&w); t[s] = w;
! 1227: }
! 1228: for ( l = 0; l < n; l++ ) {
! 1229: if ( l == i )
! 1230: continue;
! 1231: u = mat[l]; h = u[j];
! 1232: for ( s = j; s <= m; s++ ) {
! 1233: mulgf2n(h,t[s],&w); addgf2n(w,u[s],&w1); u[s] = w1;
! 1234: }
! 1235: }
! 1236: }
! 1237: }
! 1238:
! 1239: /*
! 1240: void null_to_sol(mat,ind,mod,n,r)
! 1241: int **mat;
! 1242: int *ind;
! 1243: int mod,n;
! 1244: UM *r;
! 1245: {
! 1246: int i,j,k,l;
! 1247: int *c;
! 1248: UM w;
! 1249:
! 1250: for ( i = 0, l = 0; i < n; i++ ) {
! 1251: if ( !ind[i] )
! 1252: continue;
! 1253: w = UMALLOC(n);
! 1254: for ( j = k = 0, c = COEF(w); j < n; j++ )
! 1255: if ( ind[j] )
! 1256: c[j] = 0;
! 1257: else
! 1258: c[j] = mat[k++][i];
! 1259: c[i] = mod-1;
! 1260: for ( j = n; j >= 0; j-- )
! 1261: if ( c[j] )
! 1262: break;
! 1263: DEG(w) = j;
! 1264: r[l++] = w;
! 1265: }
! 1266: }
! 1267: */
! 1268:
! 1269: void showgfmat(UM **mat,int n)
! 1270: {
! 1271: int i,j,k;
! 1272: int *c;
! 1273: UM p;
! 1274:
! 1275: for ( i = 0; i < n; i++ ) {
! 1276: for ( j = 0; j < n; j++ ) {
! 1277: p = mat[i][j];
! 1278: if ( DEG(p) < 0 )
! 1279: fprintf(asir_out,"0");
! 1280: else
! 1281: for ( p = mat[i][j], k = DEG(p), c = COEF(p); k >= 0; k-- ) {
! 1282: if ( c[k] )
! 1283: fprintf(asir_out,"+%d",c[k]);
! 1284: if ( k > 1 )
! 1285: fprintf(asir_out,"a^%d",k);
! 1286: else if ( k == 1 )
! 1287: fprintf(asir_out,"a",k);
! 1288: }
! 1289: fprintf(asir_out," ");
! 1290: }
! 1291: fprintf(asir_out,"\n");
! 1292: }
! 1293: }
! 1294:
! 1295: #if 0
! 1296: void Pgcda_mod(arg,rp)
! 1297: NODE arg;
! 1298: P *rp;
! 1299: {
! 1300: p1 = (P)ARG0(arg);
! 1301: p2 = (P)ARG1(arg);
! 1302: v = VR((P)ARG2(arg));
! 1303: d = (P)ARG3(arg);
! 1304: m = QTOS((Q)ARG4(arg));
! 1305: reordvar(CO,v,&vl);
! 1306: reorderp(vl,CO,p1,&t); ptomp(m,t,&m1);
! 1307: reorderp(vl,CO,p2,&t); ptomp(m,t,&m2);
! 1308: if ( NUM(m1) || NUM(m2) || VR(m1) != v || VR(m2) != v ) {
! 1309: *rp = ONE; return;
! 1310: }
! 1311: if ( deg(v,m1) >= deg(v,m2) ) {
! 1312: t = m1; m1 = m2; m2 = t;
! 1313: }
! 1314: while ( 1 ) {
! 1315: inva_mod(COEF(DC(m2)),d,m,&inv);
! 1316: NEWDC(dc); NEXT(dc) = 0; MKP(v,dc,h);
! 1317: d0 = deg(v,m1)-deg(v,m2); STOQ(d0,DEG(dc));
! 1318: mulgq(m,d,inv,COEF(DC(m1)),&COEF(dc));
! 1319: mulgp(vl,m,d,m2,h,&t); subgp(vl,m,d,m1,t,&s);
! 1320: }
! 1321: }
! 1322: #endif
! 1323:
! 1324: void Ppwr_mod(NODE arg,P *rp)
! 1325: {
! 1326: P p,a,d,r;
! 1327: int m;
! 1328: Z n;
! 1329:
! 1330: m = QTOS((Q)ARG4(arg)); n = (Z)ARG5(arg);
! 1331: ptomp(m,(P)ARG0(arg),&p); ptomp(m,(P)ARG1(arg),&a);
! 1332: ptomp(m,(P)ARG3(arg),&d);
! 1333: pwr_mod(p,a,VR((P)ARG2(arg)),d,m,n,&r);
! 1334: mptop(r,rp);
! 1335: }
! 1336:
! 1337: void pwr_mod(P p,P a,V v,P d,int m,Z n,P *rp)
! 1338: {
! 1339: P t,s,r;
! 1340: Z n1,two,b;
! 1341:
! 1342: if ( !n )
! 1343: *rp = (P)ONEM;
! 1344: else if ( UNIZ(n) )
! 1345: *rp = p;
! 1346: else {
! 1347: STOQ(2,two);
! 1348: divqrz(n,two,&n1,&b);
! 1349: pwr_mod(p,a,v,d,m,n1,&t);
! 1350: mulmp(CO,m,t,t,&s); rem_mod(s,a,v,d,m,&r);
! 1351: if ( b ) {
! 1352: mulmp(CO,m,r,p,&t); rem_mod(t,a,v,d,m,rp);
! 1353: } else
! 1354: *rp = r;
! 1355: }
! 1356: }
! 1357:
! 1358: void rem_mod(P p,P a,V v,P d,int m,P *rp)
! 1359: {
! 1360: P tmp,r1,r2;
! 1361:
! 1362: divsrmp(CO,m,p,d,&tmp,&r1);
! 1363: divsrmp(CO,m,r1,a,&tmp,&r2);
! 1364: divsrmp(CO,m,r2,d,&tmp,rp);
! 1365: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to gf.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib2 / asir2018 / builtin