Annotation of OpenXM_contrib2/asir2018/engine/Hgfs.c, Revision 1.1
1.1 ! noro 1: /* $OpenXM$ */
! 2:
! 3: #include "ca.h"
! 4: #include "inline.h"
! 5:
! 6: int debug_sfbfctr;
! 7:
! 8: void lnfsf(int n,UM p0,UM p1,struct p_pair *list,UM np0,UM np1);
! 9: void extractcoefbm(BM f,int dx,UM r);
! 10:
! 11: int comp_dum(DUM a,DUM b)
! 12: {
! 13: if ( DEG(a->f) > DEG(b->f) )
! 14: return -1;
! 15: else if ( DEG(a->f) < DEG(b->f) )
! 16: return 1;
! 17: else
! 18: return 0;
! 19: }
! 20:
! 21: void ufctrsf(P p,DCP *dcp)
! 22: {
! 23: int n,i,j,k;
! 24: DCP dc,dc0;
! 25: P lc;
! 26: UM mp;
! 27: UM *tl;
! 28: Obj obj;
! 29: struct oDUM *udc,*udc1;
! 30:
! 31: simp_ff((Obj)p,&obj); p = (P)obj;
! 32: if ( !p ) {
! 33: NEWDC(dc); COEF(dc) = 0; DEG(dc) = ONE;
! 34: NEXT(dc) = 0; *dcp = dc;
! 35: return;
! 36: }
! 37: mp = W_UMALLOC(UDEG(p));
! 38: ptosfum(p,mp);
! 39: if ( (n = DEG(mp)) < 0 ) {
! 40: *dcp = 0; return;
! 41: } else if ( n == 0 ) {
! 42: NEWDC(dc); COEF(dc) = p; DEG(dc) = ONE;
! 43: NEXT(dc) = 0; *dcp = dc;
! 44: return;
! 45: }
! 46: lc = COEF(DC(p));
! 47: if ( !_isonesf(COEF(mp)[n]) ) {
! 48: monicsfum(mp);
! 49: }
! 50:
! 51: W_CALLOC(n+1,struct oDUM,udc);
! 52: gensqfrsfum(mp,udc);
! 53:
! 54: tl = (UM *)ALLOCA((n+1)*sizeof(UM));
! 55: W_CALLOC(DEG(mp)+1,struct oDUM,udc1);
! 56:
! 57: for ( i = 0,j = 0; udc[i].f; i++ )
! 58: if ( DEG(udc[i].f) == 1 ) {
! 59: udc1[j].f = udc[i].f; udc1[j].n = udc[i].n; j++;
! 60: } else {
! 61: bzero((char *)tl,(n+1)*sizeof(UM));
! 62: czsfum(udc[i].f,tl);
! 63: for ( k = 0; tl[k]; k++, j++ ) {
! 64: udc1[j].f = tl[k]; udc1[j].n = udc[i].n;
! 65: }
! 66: }
! 67: udc = udc1;
! 68: for ( i = 0; udc[i].f; i++ );
! 69: qsort(udc,i,sizeof(struct oDUM),
! 70: (int (*)(const void *,const void *))comp_dum);
! 71:
! 72: NEWDC(dc0); COEF(dc0) = lc; DEG(dc0) = ONE; dc = dc0;
! 73: for ( n = 0; udc[n].f; n++ ) {
! 74: NEWDC(NEXT(dc)); dc = NEXT(dc);
! 75: STOQ(udc[n].n,DEG(dc)); sfumtop(VR(p),udc[n].f,&COEF(dc));
! 76: }
! 77: NEXT(dc) = 0; *dcp = dc0;
! 78: }
! 79:
! 80: void gensqfrsfum(UM p,struct oDUM *dc)
! 81: {
! 82: int n,i,j,d,mod;
! 83: UM t,s,g,f,f1,b;
! 84: GFS u,v;
! 85:
! 86: if ( (n = DEG(p)) == 1 ) {
! 87: dc[0].f = UMALLOC(DEG(p)); cpyum(p,dc[0].f); dc[0].n = 1;
! 88: return;
! 89: }
! 90: t = W_UMALLOC(n); s = W_UMALLOC(n); g = W_UMALLOC(n);
! 91: f = W_UMALLOC(n); f1 = W_UMALLOC(n); b = W_UMALLOC(n);
! 92: diffsfum(p,t); cpyum(p,s); gcdsfum(t,s,g);
! 93: if ( !DEG(g) ) {
! 94: dc[0].f = UMALLOC(DEG(p)); cpyum(p,dc[0].f); dc[0].n = 1;
! 95: return;
! 96: }
! 97: cpyum(p,b); cpyum(p,t); divsfum(t,g,f);
! 98: for ( i = 0, d = 0; DEG(f); i++ ) {
! 99: while ( 1 ) {
! 100: cpyum(b,t);
! 101: if ( divsfum(t,f,s) >= 0 )
! 102: break;
! 103: else {
! 104: cpyum(s,b); d++;
! 105: }
! 106: }
! 107: cpyum(b,t); cpyum(f,s); gcdsfum(t,s,f1);
! 108: divsfum(f,f1,s); cpyum(f1,f);
! 109: dc[i].f = UMALLOC(DEG(s)); cpyum(s,dc[i].f); dc[i].n = d;
! 110: }
! 111: mod = characteristic_sf();
! 112: if ( DEG(b) > 0 ) {
! 113: d = 1;
! 114: while ( 1 ) {
! 115: cpyum(b,t);
! 116: for ( j = DEG(t); j >= 0; j-- )
! 117: if ( COEF(t)[j] && (j % mod) )
! 118: break;
! 119: if ( j >= 0 )
! 120: break;
! 121: else {
! 122: DEG(s) = DEG(t)/mod;
! 123: for ( j = 0; j <= DEG(t); j++ ) {
! 124: iftogfs(COEF(t)[j*mod],&u);
! 125: pthrootgfs(u,&v);
! 126: COEF(s)[j] = v?FTOIF(CONT(v)):0;
! 127: }
! 128: cpyum(s,b); d *= mod;
! 129: }
! 130: }
! 131: gensqfrsfum(b,dc+i);
! 132: for ( j = i; dc[j].f; j++ )
! 133: dc[j].n *= d;
! 134: }
! 135: }
! 136:
! 137: void randsfum(int d,UM p)
! 138: {
! 139: int i;
! 140:
! 141: for ( i = 0; i < d; i++ )
! 142: COEF(p)[i] = _randomsf();
! 143: for ( i = d-1; i >= 0 && !COEF(p)[i]; i-- );
! 144: p->d = i;
! 145: }
! 146:
! 147: void pwrmodsfum(UM p,int e,UM f,UM pr)
! 148: {
! 149: UM wt,ws,q;
! 150:
! 151: if ( e == 0 ) {
! 152: DEG(pr) = 0; COEF(pr)[0] = _onesf();
! 153: } else if ( DEG(p) < 0 )
! 154: DEG(pr) = -1;
! 155: else if ( e == 1 ) {
! 156: q = W_UMALLOC(DEG(p)); cpyum(p,pr);
! 157: DEG(pr) = divsfum(pr,f,q);
! 158: } else if ( DEG(p) == 0 ) {
! 159: DEG(pr) = 0; COEF(pr)[0] = _pwrsf(COEF(p)[0],e);
! 160: } else {
! 161: wt = W_UMALLOC(2*DEG(f)); ws = W_UMALLOC(2*DEG(f));
! 162: q = W_UMALLOC(2*DEG(f));
! 163: pwrmodsfum(p,e/2,f,wt);
! 164: if ( !(e%2) ) {
! 165: mulsfum(wt,wt,pr); DEG(pr) = divsfum(pr,f,q);
! 166: } else {
! 167: mulsfum(wt,wt,ws);
! 168: DEG(ws) = divsfum(ws,f,q);
! 169: mulsfum(ws,p,pr);
! 170: DEG(pr) = divsfum(pr,f,q);
! 171: }
! 172: }
! 173: }
! 174:
! 175: void spwrsfum(UM m,UM f,Z e,UM r)
! 176: {
! 177: UM t,s,q;
! 178: Z e1,rem,two;
! 179:
! 180: if ( !e ) {
! 181: DEG(r) = 0; COEF(r)[0] = _onesf();
! 182: } else if ( UNIZ(e) )
! 183: cpyum(f,r);
! 184: else {
! 185: STOQ(2,two);
! 186: divqrz(e,two,&e1,&rem);
! 187: t = W_UMALLOC(2*DEG(m)); spwrsfum(m,f,e1,t);
! 188: s = W_UMALLOC(2*DEG(m)); q = W_UMALLOC(2*DEG(m));
! 189: mulsfum(t,t,s); DEG(s) = divsfum(s,m,q);
! 190: if ( rem ) {
! 191: mulsfum(s,f,t); DEG(t) = divsfum(t,m,q); cpyum(t,r);
! 192: } else
! 193: cpyum(s,r);
! 194: }
! 195: }
! 196:
! 197: void tracemodsfum(UM m,UM f,int e,UM r)
! 198: {
! 199: UM t,s,q,u;
! 200: int i;
! 201:
! 202: q = W_UMALLOC(2*DEG(m)+DEG(f)); /* XXX */
! 203: t = W_UMALLOC(2*DEG(m));
! 204: s = W_UMALLOC(2*DEG(m));
! 205: u = W_UMALLOC(2*DEG(m));
! 206: DEG(f) = divsfum(f,m,q);
! 207: cpyum(f,s);
! 208: cpyum(f,t);
! 209: for ( i = 1; i < e; i++ ) {
! 210: mulsfum(t,t,u);
! 211: DEG(u) = divsfum(u,m,q); cpyum(u,t);
! 212: addsfum(t,s,u); cpyum(u,s);
! 213: }
! 214: cpyum(s,r);
! 215: }
! 216:
! 217: void make_qmatsf(UM p,UM *tab,int ***mp)
! 218: {
! 219: int n,i,j;
! 220: int *c;
! 221: UM q,r;
! 222: int **mat;
! 223: int one;
! 224:
! 225: n = DEG(p);
! 226: *mp = mat = almat(n,n);
! 227: for ( j = 0; j < n; j++ ) {
! 228: r = W_UMALLOC(DEG(tab[j])); q = W_UMALLOC(DEG(tab[j]));
! 229: cpyum(tab[j],r); DEG(r) = divsfum(r,p,q);
! 230: for ( i = 0, c = COEF(r); i <= DEG(r); i++ )
! 231: mat[i][j] = c[i];
! 232: }
! 233: one = _onesf();
! 234: for ( i = 0; i < n; i++ )
! 235: mat[i][i] = _subsf(mat[i][i],one);
! 236: }
! 237:
! 238: void nullsf(int **mat,int n,int *ind)
! 239: {
! 240: int i,j,l,s,h,inv;
! 241: int *t,*u;
! 242:
! 243: bzero((char *)ind,n*sizeof(int));
! 244: ind[0] = 0;
! 245: for ( i = j = 0; j < n; i++, j++ ) {
! 246: for ( ; j < n; j++ ) {
! 247: for ( l = i; l < n; l++ )
! 248: if ( mat[l][j] )
! 249: break;
! 250: if ( l < n ) {
! 251: t = mat[i]; mat[i] = mat[l]; mat[l] = t; break;
! 252: } else
! 253: ind[j] = 1;
! 254: }
! 255: if ( j == n )
! 256: break;
! 257: inv = _invsf(mat[i][j]);
! 258: for ( s = j, t = mat[i]; s < n; s++ )
! 259: t[s] = _mulsf(t[s],inv);
! 260: for ( l = 0; l < n; l++ ) {
! 261: if ( l == i )
! 262: continue;
! 263: u = mat[l]; h = _chsgnsf(u[j]);
! 264: for ( s = j; s < n; s++ )
! 265: u[s] = _addsf(_mulsf(h,t[s]),u[s]);
! 266: }
! 267: }
! 268: }
! 269:
! 270: void null_to_solsf(int **mat,int *ind,int n,UM *r)
! 271: {
! 272: int i,j,k,l;
! 273: int *c;
! 274: UM w;
! 275:
! 276: for ( i = 0, l = 0; i < n; i++ ) {
! 277: if ( !ind[i] )
! 278: continue;
! 279: w = UMALLOC(n);
! 280: for ( j = k = 0, c = COEF(w); j < n; j++ )
! 281: if ( ind[j] )
! 282: c[j] = 0;
! 283: else
! 284: c[j] = mat[k++][i];
! 285: c[i] = _chsgnsf(_onesf());
! 286: for ( j = n; j >= 0; j-- )
! 287: if ( c[j] )
! 288: break;
! 289: DEG(w) = j;
! 290: r[l++] = w;
! 291: }
! 292: }
! 293: /*
! 294: make_qmatsf(p,tab,mp)
! 295: nullsf(mat,n,ind)
! 296: null_to_solsf(ind,n,r)
! 297: */
! 298:
! 299: void czsfum(UM f,UM *r)
! 300: {
! 301: int i,j;
! 302: int d,n,ord;
! 303: UM s,t,u,v,w,g,x,m,q;
! 304: UM *base;
! 305:
! 306: n = DEG(f); base = (UM *)ALLOCA(n*sizeof(UM));
! 307: bzero((char *)base,n*sizeof(UM));
! 308:
! 309: w = W_UMALLOC(2*n); q = W_UMALLOC(2*n); m = W_UMALLOC(2*n);
! 310:
! 311: base[0] = W_UMALLOC(0); DEG(base[0]) = 0; COEF(base[0])[0] = _onesf();
! 312:
! 313: t = W_UMALLOC(1); DEG(t) = 1; COEF(t)[0] = 0; COEF(t)[1] = _onesf();
! 314:
! 315: ord = field_order_sf();
! 316: pwrmodsfum(t,ord,f,w);
! 317: base[1] = W_UMALLOC(DEG(w));
! 318: cpyum(w,base[1]);
! 319:
! 320: for ( i = 2; i < n; i++ ) {
! 321: mulsfum(base[i-1],base[1],m);
! 322: DEG(m) = divsfum(m,f,q);
! 323: base[i] = W_UMALLOC(DEG(m)); cpyum(m,base[i]);
! 324: }
! 325:
! 326: v = W_UMALLOC(n); cpyum(f,v);
! 327: DEG(w) = 1; COEF(w)[0] = 0; COEF(w)[1] = _onesf();
! 328: x = W_UMALLOC(1); DEG(x) = 1; COEF(x)[0] = 0; COEF(x)[1] = _onesf();
! 329: t = W_UMALLOC(n); s = W_UMALLOC(n); u = W_UMALLOC(n); g = W_UMALLOC(n);
! 330:
! 331: for ( j = 0, d = 1; 2*d <= DEG(v); d++ ) {
! 332: for ( DEG(t) = -1, i = 0; i <= DEG(w); i++ )
! 333: if ( COEF(w)[i] ) {
! 334: mulssfum(base[i],COEF(w)[i],s);
! 335: addsfum(s,t,u); cpyum(u,t);
! 336: }
! 337: cpyum(t,w); cpyum(v,s); subsfum(w,x,t);
! 338: gcdsfum(s,t,g);
! 339: if ( DEG(g) >= 1 ) {
! 340: berlekampsf(g,d,base,r+j); j += DEG(g)/d;
! 341: divsfum(v,g,q); cpyum(q,v);
! 342: DEG(w) = divsfum(w,v,q);
! 343: for ( i = 0; i < DEG(v); i++ )
! 344: DEG(base[i]) = divsfum(base[i],v,q);
! 345: }
! 346: }
! 347: if ( DEG(v) ) {
! 348: r[j] = UMALLOC(DEG(v)); cpyum(v,r[j]); j++;
! 349: }
! 350: r[j] = 0;
! 351: }
! 352:
! 353: int berlekampsf(UM p,int df,UM *tab,UM *r)
! 354: {
! 355: int n,i,j,k,nf,d,nr;
! 356: int **mat;
! 357: int *ind;
! 358: UM mp,w,q,gcd,w1,w2;
! 359: UM *u;
! 360: int *root;
! 361:
! 362: n = DEG(p);
! 363: ind = ALLOCA(n*sizeof(int));
! 364: make_qmatsf(p,tab,&mat);
! 365: nullsf(mat,n,ind);
! 366: for ( i = 0, d = 0; i < n; i++ )
! 367: if ( ind[i] )
! 368: d++;
! 369: if ( d == 1 ) {
! 370: r[0] = UMALLOC(n); cpyum(p,r[0]); return 1;
! 371: }
! 372: u = ALLOCA(d*sizeof(UM *));
! 373: r[0] = UMALLOC(n); cpyum(p,r[0]);
! 374: null_to_solsf(mat,ind,n,u);
! 375: root = ALLOCA(d*sizeof(int));
! 376: w = W_UMALLOC(n); mp = W_UMALLOC(d);
! 377: w1 = W_UMALLOC(n); w2 = W_UMALLOC(n);
! 378: for ( i = 1, nf = 1; i < d; i++ ) {
! 379: minipolysf(u[i],p,mp);
! 380: nr = find_rootsf(mp,root);
! 381: for ( j = 0; j < nf; j++ ) {
! 382: if ( DEG(r[j]) == df )
! 383: continue;
! 384: for ( k = 0; k < nr; k++ ) {
! 385: cpyum(u[i],w1); cpyum(r[j],w2);
! 386: COEF(w1)[0] = _chsgnsf(root[k]);
! 387: gcdsfum(w1,w2,w);
! 388: if ( DEG(w) > 0 && DEG(w) < DEG(r[j]) ) {
! 389: gcd = UMALLOC(DEG(w));
! 390: q = UMALLOC(DEG(r[j])-DEG(w));
! 391: cpyum(w,gcd); divsfum(r[j],w,q);
! 392: r[j] = q; r[nf++] = gcd;
! 393: }
! 394: if ( nf == d )
! 395: return d;
! 396: }
! 397: }
! 398: }
! 399: /* NOT REACHED */
! 400: error("berlekampsf : cannot happen");
! 401: return 0;
! 402: }
! 403:
! 404: void minipolysf(UM f,UM p,UM mp)
! 405: {
! 406: struct p_pair *list,*l,*l1,*lprev;
! 407: int n,d;
! 408: UM u,p0,p1,np0,np1,q,w;
! 409:
! 410: list = (struct p_pair *)MALLOC(sizeof(struct p_pair));
! 411: list->p0 = u = W_UMALLOC(0); DEG(u) = 0; COEF(u)[0] = _onesf();
! 412: list->p1 = W_UMALLOC(0); cpyum(list->p0,list->p1);
! 413: list->next = 0;
! 414: n = DEG(p); w = UMALLOC(2*n);
! 415: p0 = UMALLOC(2*n); cpyum(list->p0,p0);
! 416: p1 = UMALLOC(2*n); cpyum(list->p1,p1);
! 417: q = W_UMALLOC(2*n);
! 418: while ( 1 ) {
! 419: COEF(p0)[DEG(p0)] = 0; DEG(p0)++; COEF(p0)[DEG(p0)] = _onesf();
! 420: mulsfum(f,p1,w); DEG(w) = divsfum(w,p,q); cpyum(w,p1);
! 421: np0 = UMALLOC(n); np1 = UMALLOC(n);
! 422: lnfsf(n,p0,p1,list,np0,np1);
! 423: if ( DEG(np1) < 0 ) {
! 424: cpyum(np0,mp); return;
! 425: } else {
! 426: l1 = (struct p_pair *)MALLOC(sizeof(struct p_pair));
! 427: l1->p0 = np0; l1->p1 = np1;
! 428: for ( l = list, lprev = 0, d = DEG(np1);
! 429: l && (DEG(l->p1) > d); lprev = l, l = l->next );
! 430: if ( lprev ) {
! 431: lprev->next = l1; l1->next = l;
! 432: } else {
! 433: l1->next = list; list = l1;
! 434: }
! 435: }
! 436: }
! 437: }
! 438:
! 439: void lnfsf(int n,UM p0,UM p1,struct p_pair *list,UM np0,UM np1)
! 440: {
! 441: int h,d1;
! 442: UM t0,t1,s0,s1;
! 443: struct p_pair *l;
! 444:
! 445: cpyum(p0,np0); cpyum(p1,np1);
! 446: t0 = W_UMALLOC(n); t1 = W_UMALLOC(n);
! 447: s0 = W_UMALLOC(n); s1 = W_UMALLOC(n);
! 448: for ( l = list; l; l = l->next ) {
! 449: d1 = DEG(np1);
! 450: if ( d1 == DEG(l->p1) ) {
! 451: h = _divsf(COEF(np1)[d1],_chsgnsf(COEF(l->p1)[d1]));
! 452: mulssfum(l->p0,h,t0); addsfum(np0,t0,s0); cpyum(s0,np0);
! 453: mulssfum(l->p1,h,t1); addsfum(np1,t1,s1); cpyum(s1,np1);
! 454: }
! 455: }
! 456: }
! 457:
! 458: int find_rootsf(UM p,int *root)
! 459: {
! 460: UM *r;
! 461: int i,n;
! 462:
! 463: n = DEG(p);
! 464: r = ALLOCA((DEG(p))*sizeof(UM));
! 465: canzassf(p,1,r);
! 466: for ( i = 0; i < n; i++ )
! 467: root[i] = _chsgnsf(COEF(r[i])[0]);
! 468: return n;
! 469: }
! 470:
! 471: void canzassf(UM f,int d,UM *r)
! 472: {
! 473: UM t,s,u,w,g,o;
! 474: Z n1,n2,n3,n4,n5,z;
! 475: UM *b;
! 476: int n,q,ed;
! 477:
! 478: if ( DEG(f) == d ) {
! 479: r[0] = UMALLOC(d); cpyum(f,r[0]);
! 480: return;
! 481: } else {
! 482: n = DEG(f); b = (UM *)ALLOCA(n*sizeof(UM));
! 483: bzero((char *)b,n*sizeof(UM));
! 484:
! 485: t = W_UMALLOC(2*d);
! 486: s = W_UMALLOC(DEG(f)); u = W_UMALLOC(DEG(f));
! 487: w = W_UMALLOC(DEG(f)); g = W_UMALLOC(DEG(f));
! 488: o = W_UMALLOC(0); DEG(o) = 0; COEF(o)[0] = _onesf();
! 489: q = field_order_sf();
! 490: if ( q % 2 ) {
! 491: STOQ(q,n1); STOQ(d,z); pwrz(n1,z,&n2); subz(n2,ONE,&n3);
! 492: STOQ(2,n4); divz(n3,n4,&n5);
! 493: } else
! 494: ed = d*extdeg_sf();
! 495: while ( 1 ) {
! 496: randsfum(2*d,t);
! 497: if ( q % 2 ) {
! 498: spwrsfum(f,t,n5,s); subsfum(s,o,u);
! 499: } else
! 500: tracemodsfum(f,t,ed,u);
! 501: cpyum(f,w);
! 502: gcdsfum(w,u,g);
! 503: if ( (DEG(g) >= 1) && (DEG(g) < DEG(f)) ) {
! 504: canzassf(g,d,r);
! 505: cpyum(f,w); divsfum(w,g,s);
! 506: canzassf(s,d,r+DEG(g)/d);
! 507: return;
! 508: }
! 509: }
! 510: }
! 511: }
! 512:
! 513: /* Hensel related functions */
! 514:
! 515: int sfberle(V,V,P,int,GFS *,DCP *);
! 516: void sfgcdgen(P,ML,ML *);
! 517: void sfhenmain2(BM,UM,UM,int,BM *);
! 518: void ptosfbm(int,P,BM);
! 519: void sfhensel(int count,P f,V x,V y,int degbound,GFS *evp,P *sfp,ML *listp);
! 520:
! 521: /* f = f(x,y) */
! 522:
! 523: void sfhensel(int count,P f,V x,V y,int degbound,GFS *evp,P *sfp,ML *listp)
! 524: {
! 525: int i;
! 526: int fn;
! 527: ML rlist;
! 528: BM fl;
! 529: VL vl,nvl;
! 530: int dx,dy,bound;
! 531: GFS ev;
! 532: P f1,t,c,sf;
! 533: DCP dc,dct,dc0;
! 534: UM q,fm,hm;
! 535: UM *gm;
! 536: struct oEGT tmp0,tmp1,eg_hensel,eg_hensel_t;
! 537:
! 538: clctv(CO,f,&vl);
! 539: if ( vl->v != x ) {
! 540: reordvar(vl,x,&nvl); reorderp(nvl,vl,f,&f1);
! 541: vl = nvl; f = f1;
! 542: }
! 543: if ( vl->next )
! 544: y = vl->next->v;
! 545: dx = getdeg(x,f);
! 546: dy = getdeg(y,f);
! 547: if ( dx == 1 ) {
! 548: *listp = rlist = MLALLOC(1); rlist->n = 1; rlist->c[0] = 0;
! 549: return;
! 550: }
! 551: fn = sfberle(x,y,f,count,&ev,&dc);
! 552: if ( fn <= 1 ) {
! 553: /* fn == 0 => short of evaluation points */
! 554: *listp = rlist = MLALLOC(1); rlist->n = fn; rlist->c[0] = 0;
! 555: return;
! 556: }
! 557: if ( degbound >= 0 ) {
! 558: /*
! 559: * reconstruct dc so that
! 560: * dc[1],... : factors satisfy degree bound
! 561: * dc[0] : product of others
! 562: */
! 563: c = dc->c; dc = NEXT(dc);
! 564: dc0 = 0;
! 565: fn = 0;
! 566: while ( dc ) {
! 567: if ( getdeg(x,COEF(dc)) <= degbound ) {
! 568: dct = NEXT(dc); NEXT(dc) = dc0; dc0 = dc; dc = dct;
! 569: fn++;
! 570: } else {
! 571: mulp(vl,COEF(dc),c,&t); c = t;
! 572: dc = NEXT(dc);
! 573: }
! 574: }
! 575: if ( OID(c) == O_P ) {
! 576: NEWDC(dc); COEF(dc) = c; DEG(dc) = ONE; NEXT(dc) = dc0;
! 577: fn++;
! 578: } else {
! 579: mulp(vl,dc0->c,c,&t); dc0->c = t; dc = dc0;
! 580: }
! 581: } else {
! 582: /* pass the the leading coeff. to the first element */
! 583: c = dc->c; dc = NEXT(dc);
! 584: mulp(vl,dc->c,c,&t); dc->c = t;
! 585: }
! 586:
! 587: /* convert mod y-a factors into UM */
! 588: gm = (UM *)ALLOCA(fn*sizeof(UM));
! 589: for ( i = 0; i < fn; i++, dc = NEXT(dc) ) {
! 590: gm[i] = W_UMALLOC(UDEG(dc->c));
! 591: ptosfum(dc->c,gm[i]);
! 592: }
! 593:
! 594: /* set bound */
! 595: /* g | f, lc_y(g) = lc_y(f) => deg_y(g) <= deg_y(f) */
! 596: /* so, bound = dy is sufficient, but we use slightly large value */
! 597: bound = dy+2;
! 598:
! 599: /* f(x,y) -> f(x,y+ev) */
! 600: fl = BMALLOC(dx,bound);
! 601: ptosfbm(bound,f,fl);
! 602: if ( ev ) shiftsfbm(fl,FTOIF(CONT(ev)));
! 603:
! 604: /* sf = f(x+ev) */
! 605: sfbmtop(fl,x,y,&sf);
! 606:
! 607: /* fm = fl mod y */
! 608: fm = W_UMALLOC(dx);
! 609: cpyum(COEF(fl)[0],fm);
! 610: hm = W_UMALLOC(dx);
! 611:
! 612: q = W_UMALLOC(dx);
! 613: rlist = MLALLOC(fn); rlist->n = fn; rlist->bound = bound;
! 614: if ( debug_sfbfctr )
! 615: fprintf(asir_out,"%d candidates\n",fn);
! 616: init_eg(&eg_hensel);
! 617: for ( i = 0; i < fn-1; i++ ) {
! 618: if ( debug_sfbfctr )
! 619: fprintf(asir_out,"deg(fm) = %d, deg(gm[%d]) = %d\n",
! 620: DEG(fm),i,DEG(gm[i]));
! 621: init_eg(&eg_hensel_t);
! 622: get_eg(&tmp0);
! 623: /* fl = gm[i]*hm mod y */
! 624: divsfum(fm,gm[i],hm);
! 625: /* fl is replaced by the cofactor of gk mod y^bound */
! 626: /* rlist->c[i] = gk */
! 627: sfhenmain2(fl,gm[i],hm,bound,(BM *)&rlist->c[i]);
! 628: cpyum(hm,fm);
! 629: get_eg(&tmp1); add_eg(&eg_hensel_t,&tmp0,&tmp1);
! 630: add_eg(&eg_hensel,&tmp0,&tmp1);
! 631: if ( debug_sfbfctr) {
! 632: print_eg("Hensel",&eg_hensel_t);
! 633: fprintf(asir_out,"\n");
! 634: }
! 635: }
! 636: if ( debug_sfbfctr) {
! 637: print_eg("Hensel total",&eg_hensel);
! 638: fprintf(asir_out,"\n");
! 639: }
! 640: /* finally, fl must be the lift of gm[fn-1] */
! 641: rlist->c[i] = fl;
! 642:
! 643: #if 0
! 644: /* y -> y-a */
! 645: mev = _chsgnsf(FTOIF(CONT(ev)));
! 646: for ( i = 0; i < fn; i++ )
! 647: shiftsfbm((BM)(rlist->c[i]),mev);
! 648: #endif
! 649: *evp = ev;
! 650: *sfp = sf;
! 651: *listp = rlist;
! 652: }
! 653:
! 654: /* main variable of f = x */
! 655:
! 656: int sfberle(V x,V y,P f,int count,GFS *ev,DCP *dcp)
! 657: {
! 658: UM wf,wf1,wf2,wfs,gcd;
! 659: int fn,n;
! 660: GFS m,fm;
! 661: DCP dc,dct,dc0;
! 662: VL vl;
! 663: P lc,lc0,f0;
! 664: Obj obj;
! 665: int j,q,index,i;
! 666:
! 667: NEWVL(vl); vl->v = x;
! 668: NEWVL(NEXT(vl)); NEXT(vl)->v = y;
! 669: NEXT(NEXT(vl)) =0;
! 670: simp_ff((Obj)f,&obj); f = (P)obj;
! 671: n = QTOS(DEG(DC(f)));
! 672: wf = W_UMALLOC(n); wf1 = W_UMALLOC(n); wf2 = W_UMALLOC(n);
! 673: wfs = W_UMALLOC(n); gcd = W_UMALLOC(n);
! 674: q = field_order_sf();
! 675: lc = DC(f)->c;
! 676: for ( j = 0, fn = n + 1, index = 0;
! 677: index < q && j < count && fn > 1; index++ ) {
! 678: indextogfs(index,&m);
! 679: substp(vl,lc,y,(P)m,&lc0);
! 680: if ( lc0 ) {
! 681: substp(vl,f,y,(P)m,&f0);
! 682: ptosfum(f0,wf); cpyum(wf,wf1);
! 683: diffsfum(wf1,wf2); gcdsfum(wf1,wf2,gcd);
! 684: if ( DEG(gcd) == 0 ) {
! 685: ufctrsf(f0,&dc);
! 686: for ( dct = NEXT(dc), i = 0; dct; dct = NEXT(dct), i++ );
! 687: if ( i < fn ) {
! 688: dc0 = dc; fn = i; fm = m;
! 689: }
! 690: j++;
! 691: }
! 692: }
! 693: }
! 694: if ( index == q )
! 695: return 0;
! 696: else if ( fn == 1 )
! 697: return 1;
! 698: else {
! 699: *dcp = dc0;
! 700: *ev = fm;
! 701: return fn;
! 702: }
! 703: }
! 704:
! 705: void sfgcdgen(P f,ML blist,ML *clistp)
! 706: {
! 707: int i;
! 708: int n,d,np;
! 709: UM wf,wm,wx,wy,wu,wv,wa,wb,wg,q,tum;
! 710: UM *in,*out;
! 711: ML clist;
! 712:
! 713: n = UDEG(f); np = blist->n;
! 714: d = 2*n;
! 715: q = W_UMALLOC(d); wf = W_UMALLOC(d);
! 716: wm = W_UMALLOC(d); wx = W_UMALLOC(d);
! 717: wy = W_UMALLOC(d); wu = W_UMALLOC(d);
! 718: wv = W_UMALLOC(d); wg = W_UMALLOC(d);
! 719: wa = W_UMALLOC(d); wb = W_UMALLOC(d);
! 720: ptosfum(f,wf); DEG(wg) = 0; COEF(wg)[0] = _onesf();
! 721: *clistp = clist = MLALLOC(np); clist->n = np;
! 722: for ( i = 0, in = (UM *)blist->c, out = (UM *)clist->c; i < np; i++ ) {
! 723: divsfum(wf,in[i],q); tum = wf; wf = q; q = tum;
! 724: cpyum(wf,wx); cpyum(in[i],wy);
! 725: eucsfum(wx,wy,wa,wb); mulsfum(wa,wg,wm);
! 726: DEG(wm) = divsfum(wm,in[i],q); out[i] = UMALLOC(DEG(wm));
! 727: cpyum(wm,out[i]); mulsfum(q,wf,wu);
! 728: mulsfum(wg,wb,wv); addsfum(wu,wv,wg);
! 729: }
! 730: }
! 731:
! 732: /* f = g0*h0 mod y -> f = gk*hk mod y^(dy+1), f is replaced by hk */
! 733:
! 734: void sfhenmain2(BM f,UM g0,UM h0,int dy,BM *gp)
! 735: {
! 736: int i,k;
! 737: int dx;
! 738: UM wt,wa,wb,q,w1,w2,wh1,wg1,ws;
! 739: UM wc,wd,we,wz;
! 740: BM wb0,wb1;
! 741: int dg,dh;
! 742: BM fk,gk,hk;
! 743:
! 744: if ( DEG(f) < dy )
! 745: error("sfhenmain2 : invalid input");
! 746:
! 747: dx = degbm(f);
! 748: dg = DEG(g0);
! 749: dh = DEG(h0);
! 750:
! 751: W_BMALLOC(dx,dy,wb0); W_BMALLOC(dx,dy,wb1);
! 752: wt = W_UMALLOC(dx); ws = W_UMALLOC(dx); q = W_UMALLOC(2*dx);
! 753: wg1 = W_UMALLOC(2*dx); wh1 = W_UMALLOC(2*dx);
! 754:
! 755: /* fk = gk*hk mod y^k */
! 756: W_BMALLOC(dx,dy,fk);
! 757: cpyum(COEF(f)[0],COEF(fk)[0]);
! 758: gk = BMALLOC(dg,dy);
! 759: cpyum(g0,COEF(gk)[0]);
! 760: W_BMALLOC(dh,dy,hk);
! 761: cpyum(h0,COEF(hk)[0]);
! 762:
! 763: wc = W_UMALLOC(2*dx); wd = W_UMALLOC(2*dx);
! 764: we = W_UMALLOC(2*dx); wz = W_UMALLOC(2*dx);
! 765:
! 766: /* compute wa,wb s.t. wa*g0+wb*h0 = 1 mod y */
! 767: w1 = W_UMALLOC(dg); cpyum(g0,w1);
! 768: w2 = W_UMALLOC(dh); cpyum(h0,w2);
! 769: wa = W_UMALLOC(2*dx); wb = W_UMALLOC(2*dx); /* XXX */
! 770: eucsfum(w1,w2,wa,wb);
! 771:
! 772: if ( debug_sfbfctr)
! 773: fprintf(stderr,"dy=%d\n",dy);
! 774: for ( k = 1; k <= dy; k++ ) {
! 775: if ( debug_sfbfctr)
! 776: fprintf(stderr,".");
! 777:
! 778: /* at this point, f = gk*hk mod y^k */
! 779:
! 780: /* clear wt */
! 781: clearum(wt,dx);
! 782:
! 783: /* wt = (f-gk*hk)/y^k */
! 784: subsfum(COEF(f)[k],COEF(fk)[k],wt);
! 785:
! 786: /* compute wf1,wg1 s.t. wh1*g0+wg1*h0 = wt */
! 787: mulsfum(wa,wt,wh1); DEG(wh1) = divsfum(wh1,h0,q);
! 788: mulsfum(wh1,g0,wc); subsfum(wt,wc,wd); DEG(wd) = divsfum(wd,h0,wg1);
! 789:
! 790: /* check */
! 791: #if 0
! 792: if ( DEG(wd) >= 0 || DEG(wg1) > ng )
! 793: error("henmain2 : cannot happen(adj)");
! 794:
! 795: mulsfum(wg1,h0,wc); mulsfum(wh1,g0,wd); addsfum(wc,wd,we);
! 796: subsfum(we,wt,wz);
! 797: if ( DEG(wz) >= 0 )
! 798: error("henmain2 : cannot happen");
! 799: #endif
! 800:
! 801: /* fk += ((wg1*hk+wh1*gk)*y^k+wg1*wh1*y^(2*k) mod y^(dy+1) */
! 802: /* wb0 = wh1*y^k */
! 803: clearbm(dx,wb0);
! 804: cpyum(wh1,COEF(wb0)[k]);
! 805:
! 806: /* wb1 = gk*wb0 mod y^(dy+1) */
! 807: clearbm(dx,wb1);
! 808: mulsfbm(gk,wb0,wb1);
! 809: /* fk += wb1 */
! 810: addtosfbm(wb1,fk);
! 811:
! 812: /* wb0 = wg1*y^k */
! 813: clearbm(dx,wb0);
! 814: cpyum(wg1,COEF(wb0)[k]);
! 815:
! 816: /* wb1 = hk*wb0 mod y^(dy+1) */
! 817: clearbm(dx,wb1);
! 818: mulsfbm(hk,wb0,wb1);
! 819: /* fk += wb1 */
! 820: addtosfbm(wb1,fk);
! 821:
! 822: /* fk += wg1*wh1*y^(2*k) mod y^(dy+1) */
! 823: if ( 2*k <= dy ) {
! 824: mulsfum(wg1,wh1,wt); addsfum(COEF(fk)[2*k],wt,ws);
! 825: cpyum(ws,COEF(fk)[2*k]);
! 826: }
! 827:
! 828: /* gk += wg1*y^k, hk += wh1*y^k */
! 829: cpyum(wg1,COEF(gk)[k]);
! 830: cpyum(wh1,COEF(hk)[k]);
! 831: }
! 832: if ( debug_sfbfctr)
! 833: fprintf(stderr,"\n");
! 834: *gp = gk;
! 835: DEG(f) = dy;
! 836: for ( i = 0; i <= dy; i++ )
! 837: cpyum(COEF(hk)[i],COEF(f)[i]);
! 838: }
! 839:
! 840: /* a0*g+b0*h = 1 mod y -> a*g+b*h = 1 mod y^(dy+1) */
! 841:
! 842: void sfexgcd_by_hensel(BM g,BM h,int dy,BM *ap,BM *bp)
! 843: {
! 844: int i,k;
! 845: int dx;
! 846: UM wt,wa,wb,q,w1,w2,ws;
! 847: UM wc,wd,we,wz,wa1,wb1;
! 848: BM wz0,wz1;
! 849: int dg,dh;
! 850: BM a,b,c;
! 851:
! 852: dg = degbm(g);
! 853: dh = degbm(h);
! 854: dx = dg+dh;
! 855:
! 856: a = BMALLOC(dh,dy);
! 857: b = BMALLOC(dg,dy);
! 858: /* c holds a*g+b*h-1 */
! 859: c = BMALLOC(dg+dh,dy);
! 860:
! 861: W_BMALLOC(dx,dy,wz0); W_BMALLOC(dx,dy,wz1);
! 862:
! 863: wt = W_UMALLOC(dx); ws = W_UMALLOC(dx); q = W_UMALLOC(2*dx);
! 864: wa1 = W_UMALLOC(2*dx); wb1 = W_UMALLOC(2*dx);
! 865: wc = W_UMALLOC(2*dx); wd = W_UMALLOC(2*dx);
! 866: we = W_UMALLOC(2*dx); wz = W_UMALLOC(2*dx);
! 867:
! 868: /* compute wa,wb s.t. wa*g0+wb*h0 = 1 mod y */
! 869: w1 = W_UMALLOC(dg); cpyum(COEF(g)[0],w1);
! 870: w2 = W_UMALLOC(dh); cpyum(COEF(h)[0],w2);
! 871: wa = W_UMALLOC(2*dx); wb = W_UMALLOC(2*dx); /* XXX */
! 872: eucsfum(w1,w2,wa,wb);
! 873: cpyum(wa,COEF(a)[0]); cpyum(wb,COEF(b)[0]);
! 874:
! 875: /* initialize c to a*g+b*h-1 */
! 876: mulsfbm(a,g,c); mulsfbm(b,h,wz0); addtosfbm(wz0,c);
! 877: COEF(COEF(c)[0])[0] = 0;
! 878:
! 879: if ( debug_sfbfctr)
! 880: fprintf(stderr,"dy=%d\n",dy);
! 881: for ( k = 1; k <= dy; k++ ) {
! 882: if ( debug_sfbfctr)
! 883: fprintf(stderr,".");
! 884:
! 885: /* at this point, a*g+b*h = 1 mod y^k, c = a*g+b*h-1 */
! 886:
! 887: /* wt = -((a*g+b*h-1)/y^k) */
! 888: cpyum(COEF(c)[k],wt);
! 889: for ( i = DEG(wt); i >= 0; i-- )
! 890: COEF(wt)[i] = _chsgnsf(COEF(wt)[i]);
! 891:
! 892: /* compute wa1,wb1 s.t. wa1*g0+wb1*h0 = wt */
! 893: mulsfum(wa,wt,wa1); DEG(wa1) = divsfum(wa1,COEF(h)[0],q);
! 894: mulsfum(wa1,COEF(g)[0],wc); subsfum(wt,wc,wd);
! 895: DEG(wd) = divsfum(wd,COEF(h)[0],wb1);
! 896:
! 897: /* c += ((wa1*g+wb1*h)*y^k mod y^(dy+1) */
! 898: /* wz0 = wa1*y^k */
! 899: clearbm(dx,wz0);
! 900: cpyum(wa1,COEF(wz0)[k]);
! 901:
! 902: /* wz1 = wz0*g mod y^(dy+1) */
! 903: clearbm(dx,wz1);
! 904: mulsfbm(g,wz0,wz1);
! 905: /* c += wz1 */
! 906: addtosfbm(wz1,c);
! 907:
! 908: /* wz0 = wb1*y^k */
! 909: clearbm(dx,wz0);
! 910: cpyum(wb1,COEF(wz0)[k]);
! 911:
! 912: /* wz1 = wz0*h mod y^(dy+1) */
! 913: clearbm(dx,wz1);
! 914: mulsfbm(h,wz0,wz1);
! 915: /* c += wz1 */
! 916: addtosfbm(wz1,c);
! 917:
! 918: /* a += wa1*y^k, b += wb1*y^k */
! 919: cpyum(wa1,COEF(a)[k]);
! 920: cpyum(wb1,COEF(b)[k]);
! 921: }
! 922: if ( debug_sfbfctr)
! 923: fprintf(stderr,"\n");
! 924: DEG(a) = dy;
! 925: DEG(b) = dy;
! 926: *ap = a;
! 927: *bp = b;
! 928: }
! 929:
! 930: /* fl->c[i] = coef_y(f,i) */
! 931:
! 932: void ptosfbm(int dy,P f,BM fl)
! 933: {
! 934: DCP dc;
! 935: int d,i,dx;
! 936: UM t;
! 937:
! 938: dx = QTOS(DEG(DC(f)));
! 939: if ( DEG(fl) < dy )
! 940: error("ptosfbm : invalid input");
! 941: DEG(fl) = dy;
! 942: clearbm(dx,fl);
! 943: t = UMALLOC(dy);
! 944: for ( dc = DC(f); dc; dc = NEXT(dc) ) {
! 945: d = QTOS(DEG(dc));
! 946: ptosfum(COEF(dc),t);
! 947: for ( i = 0; i <= DEG(t); i++ )
! 948: COEF(COEF(fl)[i])[d] = COEF(t)[i];
! 949: }
! 950: for ( i = 0; i <= dy; i++ )
! 951: degum(COEF(fl)[i],dx);
! 952: }
! 953:
! 954: /* x : main variable */
! 955:
! 956: void sfbmtop(BM f,V x,V y,P *fp)
! 957: {
! 958: UM *c;
! 959: int i,j,d,a,dy;
! 960: GFS b;
! 961: DCP dc0,dc,dct;
! 962:
! 963: dy = DEG(f);
! 964: c = COEF(f);
! 965: d = degbm(f);
! 966:
! 967: dc0 = 0;
! 968: for ( i = 0; i <= d; i++ ) {
! 969: dc = 0;
! 970: for ( j = 0; j <= dy; j++ ) {
! 971: if ( DEG(c[j]) >= i && (a = COEF(c[j])[i]) ) {
! 972: NEWDC(dct);
! 973: STOQ(j,DEG(dct));
! 974: iftogfs(a,&b);
! 975: COEF(dct) = (P)b;
! 976: NEXT(dct) = dc;
! 977: dc = dct;
! 978: }
! 979: }
! 980: if ( dc ) {
! 981: NEWDC(dct);
! 982: STOQ(i,DEG(dct));
! 983: MKP(y,dc,COEF(dct));
! 984: NEXT(dct) = dc0;
! 985: dc0 = dct;
! 986: }
! 987: }
! 988: if ( dc0 )
! 989: MKP(x,dc0,*fp);
! 990: else
! 991: *fp = 0;
! 992: }
! 993:
! 994: void sfsqfr(P f,DCP *dcp)
! 995: {
! 996: Obj obj;
! 997: DCP dc;
! 998: VL vl;
! 999:
! 1000: simp_ff((Obj)f,&obj); f = (P)obj;
! 1001: clctv(CO,f,&vl);
! 1002: if ( !vl ) {
! 1003: /* f is a const */
! 1004: NEWDC(dc); DEG(dc) = ONE; COEF(dc) = f; NEXT(dc) = 0; *dcp = dc;
! 1005: } else if ( !NEXT(vl) )
! 1006: sfusqfr(f,dcp);
! 1007: else
! 1008: sqfrsf(vl,f,dcp);
! 1009: }
! 1010:
! 1011: void sfusqfr(P f,DCP *dcp)
! 1012: {
! 1013: DCP dc,dct;
! 1014: struct oDUM *udc;
! 1015: V x;
! 1016: P lc;
! 1017: int n,i;
! 1018: UM mf;
! 1019:
! 1020: x = VR(f);
! 1021: n = getdeg(x,f);
! 1022: mf = W_UMALLOC(n);
! 1023: ptosfum(f,mf);
! 1024: lc = COEF(DC(f));
! 1025: if ( !_isonesf(COEF(mf)[n]) ) {
! 1026: monicsfum(mf);
! 1027: }
! 1028: W_CALLOC(n+1,struct oDUM,udc);
! 1029: gensqfrsfum(mf,udc);
! 1030: for ( i = 0, dc = 0; udc[i].f; i++ ) {
! 1031: NEWDC(dct); STOQ(udc[i].n,DEG(dct));
! 1032: sfumtop(x,udc[i].f,&COEF(dct));
! 1033: NEXT(dct) = dc; dc = dct;
! 1034: }
! 1035: NEWDC(dct); DEG(dct) = ONE; COEF(dct) = (P)lc; NEXT(dct) = dc;
! 1036: *dcp = dct;
! 1037: }
! 1038:
! 1039: #if 0
! 1040: void sfbsqfrmain(P f,V x,V y,DCP *dcp)
! 1041: {
! 1042: /* XXX*/
! 1043: }
! 1044:
! 1045: /* f is bivariate */
! 1046:
! 1047: void sfbsqfr(P f,V x,V y,DCP *dcp)
! 1048: {
! 1049: P t,rf,cx,cy;
! 1050: VL vl,rvl;
! 1051: DCP dcx,dcy,dct,dc;
! 1052: struct oVL vl0,vl1;
! 1053:
! 1054: /* cy(y) = cont(f,x), f /= cy */
! 1055: cont_pp_sfp(vl,f,&cy,&t); f = t;
! 1056: /* rvl = [y,x] */
! 1057: reordvar(vl,y,&rvl); reorderp(rvl,vl,f,&rf);
! 1058: /* cx(x) = cont(rf,y), Rf /= cy */
! 1059: cont_pp_sfp(rvl,rf,&cx,&t); rf = t;
! 1060: reorderp(vl,rvl,rf,&f);
! 1061:
! 1062: /* f -> cx*cy*f */
! 1063: sfsqfr(cx,&dcx); dcx = NEXT(dcx);
! 1064: sfsqfr(cy,&dcy); dcy = NEXT(dcy);
! 1065: if ( dcx ) {
! 1066: for ( dct = dcx; NEXT(dct); dct = NEXT(dct) );
! 1067: NEXT(dct) = dcy;
! 1068: } else
! 1069: dcx = dcy;
! 1070: if ( OID(f) == O_N )
! 1071: *dcp = dcx;
! 1072: else {
! 1073: /* f must be bivariate */
! 1074: sfbsqfrmain(f,x,y,&dc);
! 1075: if ( dcx ) {
! 1076: for ( dct = dcx; NEXT(dct); dct = NEXT(dct) );
! 1077: NEXT(dct) = dc;
! 1078: } else
! 1079: dcx = dc;
! 1080: *dcp = dcx;
! 1081: }
! 1082: }
! 1083: #endif
! 1084:
! 1085: void sfdtest(P,ML,V,V,DCP *);
! 1086:
! 1087: /* if degbound >= 0 find factor s.t. deg_x(factor) <= degbound */
! 1088:
! 1089: void sfbfctr(P f,V x,V y,int degbound,DCP *dcp)
! 1090: {
! 1091: ML list;
! 1092: P sf;
! 1093: GFS ev;
! 1094: DCP dc,dct;
! 1095: BM fl;
! 1096: int dx,dy;
! 1097:
! 1098: /* sf(x) = f(x+ev) = list->c[0]*list->c[1]*... */
! 1099: sfhensel(5,f,x,y,degbound,&ev,&sf,&list);
! 1100: if ( list->n == 0 )
! 1101: error("sfbfctr : short of evaluation points");
! 1102: else if ( list->n == 1 ) {
! 1103: /* f is irreducible */
! 1104: NEWDC(dc); DEG(dc) = ONE; COEF(dc) = f; NEXT(dc) = 0;
! 1105: *dcp = dc;
! 1106: return;
! 1107: }
! 1108: sfdtest(sf,list,x,y,&dc);
! 1109: if ( ev ) {
! 1110: dx = getdeg(x,sf);
! 1111: dy = getdeg(y,sf);
! 1112: W_BMALLOC(dx,dy,fl);
! 1113: for ( dct = dc; dct; dct = NEXT(dct) ) {
! 1114: ptosfbm(dy,COEF(dct),fl);
! 1115: shiftsfbm(fl,_chsgnsf(FTOIF(CONT(ev))));
! 1116: sfbmtop(fl,x,y,&COEF(dct));
! 1117: }
! 1118: }
! 1119: *dcp = dc;
! 1120: }
! 1121:
! 1122: /* returns shifted f, shifted factors and the eval pt */
! 1123:
! 1124: void sfbfctr_shift(P f,V x,V y,int degbound,GFS *evp,P *sfp,DCP *dcp)
! 1125: {
! 1126: ML list;
! 1127: P sf;
! 1128: GFS ev;
! 1129: DCP dc,dct;
! 1130: int dx,dy;
! 1131:
! 1132: /* sf(x) = f(x+ev) = list->c[0]*list->c[1]*... */
! 1133: sfhensel(5,f,x,y,degbound,&ev,&sf,&list);
! 1134: if ( list->n == 0 )
! 1135: error("sfbfctr_shift : short of evaluation points");
! 1136: else if ( list->n == 1 ) {
! 1137: /* f is irreducible */
! 1138: NEWDC(dc); DEG(dc) = ONE; COEF(dc) = f; NEXT(dc) = 0;
! 1139: *evp = 0;
! 1140: *sfp = f;
! 1141: *dcp = dc;
! 1142: } else {
! 1143: sfdtest(sf,list,x,y,dcp);
! 1144: *evp = ev;
! 1145: *sfp = sf;
! 1146: }
! 1147: }
! 1148:
! 1149: /* f = f(x,y) = list->c[0]*list->c[1]*... mod y^(list->bound+1) */
! 1150:
! 1151: void sfdtest(P f,ML list,V x,V y,DCP *dcp)
! 1152: {
! 1153: int np,dx,dy;
! 1154: int i,j,k,bound;
! 1155: int *win;
! 1156: P g,lcg,factor,cofactor,lcyx;
! 1157: P csum;
! 1158: DCP dcf,dcf0,dc;
! 1159: BM *c;
! 1160: BM lcy;
! 1161: UM lcg0,lcy0,w;
! 1162: UM *d1c;
! 1163: ML wlist;
! 1164: struct oVL vl1,vl0;
! 1165: VL vl;
! 1166: int z,dt,dtok;
! 1167:
! 1168: /* vl = [x,y] */
! 1169: vl0.v = x; vl0.next = &vl1; vl1.v = y; vl1.next = 0; vl = &vl0;
! 1170:
! 1171: /* setup various structures and arrays */
! 1172: dx = getdeg(x,f);
! 1173: dy = getdeg(y,f);
! 1174: np = list->n;
! 1175: win = W_ALLOC(np+1);
! 1176: wlist = W_MLALLOC(np);
! 1177: wlist->n = list->n;
! 1178: bound = wlist->bound = list->bound;
! 1179: c = (BM *)COEF(wlist);
! 1180: bcopy((char *)COEF(list),(char *)c,(int)(sizeof(BM)*np));
! 1181:
! 1182: lcg0 = W_UMALLOC(2*dy);
! 1183:
! 1184: /* initialize g by f */
! 1185: g = f;
! 1186:
! 1187: /* initialize lcg */
! 1188: mulp(vl,g,COEF(DC(g)),&lcg);
! 1189:
! 1190: /* initialize lcg0 */
! 1191: const_term(lcg,lcg0);
! 1192:
! 1193: /* initialize csum = lcg(1) */
! 1194: sfcsump(vl,lcg,&csum);
! 1195:
! 1196: /* initialize lcy by LC(f) */
! 1197: W_BMALLOC(0,dy,lcy);
! 1198: NEWDC(dc); COEF(dc) = COEF(DC(g)); DEG(dc) = 0;
! 1199: NEWP(lcyx); VR(lcyx) = x; DC(lcyx) = dc;
! 1200: ptosfbm(dy,lcyx,lcy);
! 1201:
! 1202: /* initialize lcy0 by LC(f) */
! 1203: lcy0 = W_UMALLOC(bound);
! 1204: ptosfum(COEF(DC(g)),lcy0);
! 1205:
! 1206: /* ((d-1 coefs)*lcy0 */
! 1207: d1c = (UM *)W_ALLOC(np*sizeof(UM));
! 1208: w = W_UMALLOC(2*bound);
! 1209: for ( i = 1; i < np; i++ ) {
! 1210: extractcoefbm(c[i],degbm(c[i])-1,w);
! 1211: d1c[i] = W_UMALLOC(2*bound);
! 1212: mulsfum(w,lcy0,d1c[i]);
! 1213: /* d1c[i] = d1c[i] mod y^(bound+1) */
! 1214: if ( DEG(d1c[i]) > bound ) {
! 1215: for ( j = DEG(d1c[i]); j > bound; j-- )
! 1216: COEF(d1c[i])[j] = 0;
! 1217: degum(d1c[i],bound);
! 1218: }
! 1219: }
! 1220:
! 1221: if ( debug_sfbfctr)
! 1222: fprintf(stderr,"np = %d\n",np);
! 1223: dtok = 0;
! 1224: for ( g = f, k = 1, dcf = dcf0 = 0, win[0] = 1, --np, z = 0; ; z++ ) {
! 1225: if ( debug_sfbfctr && !(z % 1000) ) fprintf(stderr,".");
! 1226: dt = sfdegtest(dy,bound,d1c,k,win);
! 1227: if ( dt )
! 1228: dtok++;
! 1229: if ( dt && sfdtestmain(vl,lcg,lcg0,lcy,csum,wlist,
! 1230: k,win,&factor,&cofactor) ) {
! 1231: NEXTDC(dcf0,dcf); DEG(dcf) = ONE; COEF(dcf) = factor;
! 1232: g = cofactor;
! 1233:
! 1234: /* update lcg */
! 1235: mulp(vl,g,COEF(DC(g)),&lcg);
! 1236:
! 1237: /* update lcg0 */
! 1238: const_term(lcg,lcg0);
! 1239:
! 1240: /* update csum */
! 1241: sfcsump(vl,lcg,&csum);
! 1242:
! 1243: /* update dy */
! 1244: dy = getdeg(y,g);
! 1245:
! 1246: /* update lcy */
! 1247: clearbm(0,lcy);
! 1248: COEF(dc) = COEF(DC(g));
! 1249: ptosfbm(dy,lcyx,lcy);
! 1250:
! 1251: for ( i = 0; i < k - 1; i++ )
! 1252: for ( j = win[i] + 1; j < win[i + 1]; j++ )
! 1253: c[j-i-1] = c[j];
! 1254: for ( j = win[k-1] + 1; j <= np; j++ )
! 1255: c[j-k] = c[j];
! 1256: if ( ( np -= k ) < k )
! 1257: break;
! 1258: if ( np - win[0] + 1 < k )
! 1259: if ( ++k > np )
! 1260: break;
! 1261: else
! 1262: for ( i = 0; i < k; i++ )
! 1263: win[i] = i + 1;
! 1264: else
! 1265: for ( i = 1; i < k; i++ )
! 1266: win[i] = win[0] + i;
! 1267:
! 1268:
! 1269: /* update lcy0 */
! 1270: ptosfum(COEF(DC(g)),lcy0);
! 1271:
! 1272: /* update d-1 coeffs */
! 1273: for ( i = 1; i <= np; i++ ) {
! 1274: extractcoefbm(c[i],degbm(c[i])-1,w);
! 1275: mulsfum(w,lcy0,d1c[i]);
! 1276: /* d1c[i] = d1c[1] mod y^(bound+1) */
! 1277: if ( DEG(d1c[i]) > bound ) {
! 1278: for ( j = DEG(d1c[i]); j > bound; j-- )
! 1279: COEF(d1c[i])[j] = 0;
! 1280: degum(d1c[i],bound);
! 1281: }
! 1282: }
! 1283: } else if ( !ncombi(1,np,k,win) ) {
! 1284: if ( k == np )
! 1285: break;
! 1286: else
! 1287: for ( i = 0, ++k; i < k; i++ )
! 1288: win[i] = i + 1;
! 1289: }
! 1290: }
! 1291: if ( debug_sfbfctr )
! 1292: fprintf(stderr,"total %d, omitted by degtest %d\n",z,z-dtok);
! 1293: NEXTDC(dcf0,dcf); COEF(dcf) = g;
! 1294: DEG(dcf) = ONE; NEXT(dcf) = 0; *dcp = dcf0;
! 1295: }
! 1296:
! 1297: void extractcoefbm(BM f,int dx,UM r)
! 1298: {
! 1299: int j;
! 1300: UM fj;
! 1301:
! 1302: for ( j = DEG(f); j >= 0; j-- ) {
! 1303: fj = COEF(f)[j];
! 1304: if ( fj && DEG(fj) >= dx ) {
! 1305: COEF(r)[j] = COEF(fj)[dx];
! 1306: } else
! 1307: COEF(r)[j] = 0;
! 1308: }
! 1309: degum(r,DEG(f));
! 1310: }
! 1311:
! 1312: /* deg_y(prod mod y^(bound+1)) <= dy ? */
! 1313:
! 1314: int sfdegtest(int dy,int bound,UM *d1c,int k,int *in)
! 1315: {
! 1316: int i,j;
! 1317: UM w,w1,wt;
! 1318: BM t;
! 1319:
! 1320: w = W_UMALLOC(bound);
! 1321: w1 = W_UMALLOC(bound);
! 1322: clearum(w,bound);
! 1323: for ( i = 0; i < k; i++ ) {
! 1324: addsfum(w,d1c[in[i]],w1); wt = w; w = w1; w1 = wt;
! 1325: }
! 1326: return DEG(w) <= dy ? 1 : 0;
! 1327: }
! 1328:
! 1329: /* lcy = LC(g), lcg = lcy*g, lcg0 = const part of lcg */
! 1330:
! 1331: int sfdtestmain(VL vl,P lcg,UM lcg0,BM lcy,P csum,ML list,
! 1332: int k,int *in,P *fp,P *cofp)
! 1333: {
! 1334: P fmul,csumg,q,cont;
! 1335: V x,y;
! 1336:
! 1337: x = vl->v;
! 1338: y = vl->next->v;
! 1339: if (!sfctest(lcg0,lcy,list,k,in))
! 1340: return 0;
! 1341: mulsfbmarray(UDEG(lcg),lcy,list,k,in,x,y,&fmul);
! 1342: if ( csum ) {
! 1343: sfcsump(vl,fmul,&csumg);
! 1344: if ( csumg ) {
! 1345: if ( !divtp(vl,csum,csumg,&q) )
! 1346: return 0;
! 1347: }
! 1348: }
! 1349: if ( divtp_by_sfbm(vl,lcg,fmul,&q) ) {
! 1350: cont_pp_sfp(vl,fmul,&cont,fp);
! 1351: cont_pp_sfp(vl,q,&cont,cofp);
! 1352: return 1;
! 1353: } else
! 1354: return 0;
! 1355: }
! 1356:
! 1357: void const_term(P f,UM c)
! 1358: {
! 1359: DCP dc;
! 1360:
! 1361: for ( dc = DC(f); dc && DEG(dc); dc = NEXT(dc) );
! 1362: if ( dc )
! 1363: ptosfum(COEF(dc),c);
! 1364: else
! 1365: DEG(c) = -1;
! 1366: }
! 1367:
! 1368: void const_term_sfbm(BM f,UM c)
! 1369: {
! 1370: int i,dy;
! 1371:
! 1372: dy = DEG(f);
! 1373: for ( i = 0; i <= dy; i++ )
! 1374: if ( DEG(COEF(f)[i]) >= 0 )
! 1375: COEF(c)[i] = COEF(COEF(f)[i])[0];
! 1376: else
! 1377: COEF(c)[i] = 0;
! 1378: degum(c,dy);
! 1379: }
! 1380:
! 1381: /* lcy*(product of const part) | lcg0 ? */
! 1382:
! 1383: int sfctest(UM lcg0,BM lcy,ML list,int k,int *in)
! 1384: {
! 1385: int dy,i,dr;
! 1386: UM t,s,u,w;
! 1387: BM *l;
! 1388:
! 1389: dy = list->bound;
! 1390: t = W_UMALLOC(2*dy);
! 1391: s = W_UMALLOC(2*dy);
! 1392: u = W_UMALLOC(2*dy);
! 1393: const_term_sfbm(lcy,t);
! 1394: if ( DEG(t) < 0 )
! 1395: return 1;
! 1396:
! 1397: l = (BM *)list->c;
! 1398: for ( i = 0; i < k; i++ ) {
! 1399: const_term_sfbm(l[in[i]],s);
! 1400: mulsfum(t,s,u);
! 1401: if ( DEG(u) > dy )
! 1402: degum(u,dy);
! 1403: w = t; t = u; u = w;
! 1404: }
! 1405: cpyum(lcg0,s);
! 1406: dr = divsfum(s,t,u);
! 1407: if ( dr >= 0 )
! 1408: return 0;
! 1409: else
! 1410: return 1;
! 1411: }
! 1412:
! 1413: /* main var of f is x */
! 1414:
! 1415: void mulsfbmarray(int dx,BM lcy,ML list,int k,int *in,V x,V y,P *g)
! 1416: {
! 1417: int dy,i;
! 1418: BM wb0,wb1,t;
! 1419: BM *l;
! 1420:
! 1421: dy = list->bound;
! 1422: W_BMALLOC(dx,dy,wb0); W_BMALLOC(dx,dy,wb1);
! 1423: l = (BM *)list->c;
! 1424: clearbm(dx,wb0);
! 1425: mulsfbm(lcy,l[in[0]],wb0);
! 1426: for ( i = 1; i < k; i++ ) {
! 1427: clearbm(dx,wb1);
! 1428: mulsfbm(l[in[i]],wb0,wb1);
! 1429: t = wb0; wb0 = wb1; wb1 = t;
! 1430: }
! 1431: sfbmtop(wb0,x,y,g);
! 1432: }
! 1433:
! 1434: void sfcsump(VL vl,P f,P *s)
! 1435: {
! 1436: P t,u;
! 1437: DCP dc;
! 1438:
! 1439: for ( dc = DC(f), t = 0; dc; dc = NEXT(dc) ) {
! 1440: addp(vl,COEF(dc),t,&u); t = u;
! 1441: }
! 1442: *s = t;
! 1443: }
! 1444:
! 1445: /* *fp = primitive part of f w.r.t. x */
! 1446:
! 1447: void cont_pp_sfp(VL vl,P f,P *cp,P *fp)
! 1448: {
! 1449: V x,y;
! 1450: int d;
! 1451: UM t,s,gcd;
! 1452: DCP dc;
! 1453: GFS g;
! 1454:
! 1455: x = vl->v;
! 1456: y = vl->next->v;
! 1457: d = getdeg(y,f);
! 1458: if ( d == 0 ) {
! 1459: itogfs(1,&g);
! 1460: *cp = (P)g;
! 1461: *fp = f; /* XXX */
! 1462: } else {
! 1463: t = W_UMALLOC(2*d);
! 1464: s = W_UMALLOC(2*d);
! 1465: gcd = W_UMALLOC(2*d);
! 1466: dc = DC(f);
! 1467: ptosfum(COEF(dc),gcd);
! 1468: for ( dc = NEXT(dc); dc; dc = NEXT(dc) ) {
! 1469: ptosfum(COEF(dc),t);
! 1470: gcdsfum(gcd,t,s);
! 1471: cpyum(s,gcd);
! 1472: }
! 1473: sfumtop(y,gcd,cp);
! 1474: divsp(vl,f,*cp,fp);
! 1475: }
! 1476: }
! 1477:
! 1478: int divtp_by_sfbm(VL vl,P f,P g,P *qp)
! 1479: {
! 1480: V x,y;
! 1481: int fx,fy,gx,gy;
! 1482: BM fl,gl,ql;
! 1483: UM *cf,*cg,*cq;
! 1484: UM hg,q,t,s;
! 1485: int i,j,dr;
! 1486:
! 1487: x = vl->v; y = vl->next->v;
! 1488: fx = getdeg(x,f); fy = getdeg(y,f);
! 1489: gx = getdeg(x,g); gy = getdeg(y,g);
! 1490:
! 1491: if ( fx < gx || fy < gy )
! 1492: return 0;
! 1493: W_BMALLOC(fx,fy,fl); ptosfbm(fy,f,fl); cf = COEF(fl);
! 1494: W_BMALLOC(gx,gy,gl); ptosfbm(gy,g,gl); cg = COEF(gl);
! 1495: W_BMALLOC(fx-gx,fy-gy,ql); cq = COEF(ql);
! 1496:
! 1497: hg = cg[gy];
! 1498: q = W_UMALLOC(fx); t = W_UMALLOC(fx); s = W_UMALLOC(fx);
! 1499:
! 1500: for ( i = fy; i >= gy; i-- ) {
! 1501: if ( DEG(cf[i]) < 0 )
! 1502: continue;
! 1503: dr = divsfum(cf[i],hg,q);
! 1504: if ( dr >= 0 )
! 1505: return 0;
! 1506: if ( DEG(q) > fx-gx )
! 1507: return 0;
! 1508: cpyum(q,cq[i-gy]);
! 1509: for ( j = 0; j <= gy; j++ ) {
! 1510: mulsfum(cg[j],q,t);
! 1511: subsfum(cf[j+i-gy],t,s);
! 1512: cpyum(s,cf[j+i-gy]);
! 1513: }
! 1514: }
! 1515: for ( j = gy-1; j >= 0 && DEG(cf[j]) < 0; j-- );
! 1516: if ( j >= 0 )
! 1517: return 0;
! 1518: sfbmtop(ql,x,y,qp);
! 1519: return 1;
! 1520: }
! 1521:
! 1522: /* XXX generate an irreducible poly of degree n */
! 1523:
! 1524: extern int current_gfs_q1;
! 1525: extern int *current_gfs_ntoi;
! 1526:
! 1527: void generate_defpoly_sfum(int n,UM *dp)
! 1528: {
! 1529: UM r,dr,t,g;
! 1530: UM *f;
! 1531: int *c,*w;
! 1532: int max,i,j;
! 1533:
! 1534: *dp = r = UMALLOC(n);
! 1535: DEG(r) = n;
! 1536: c = COEF(r);
! 1537: c[n] = _onesf();
! 1538: max = current_gfs_q1;
! 1539: w = (int *)ALLOCA(n*sizeof(int));
! 1540: bzero(w,n*sizeof(int));
! 1541:
! 1542: dr = W_UMALLOC(n); t = W_UMALLOC(n); g = W_UMALLOC(n);
! 1543: f = (UM *)ALLOCA((n+1)*sizeof(UM));
! 1544: while ( 1 ) {
! 1545: for ( i = 0; i < n && w[i] == max; i++ );
! 1546: if ( i == n ) {
! 1547: /* XXX cannot happen */
! 1548: error("generate_defpoly_sfum : cannot happen");
! 1549: }
! 1550: for ( j = 0; j < i; j++ )
! 1551: w[j] = 0;
! 1552: w[i]++;
! 1553: if ( !current_gfs_ntoi )
! 1554: for ( i = 0; i < n; i++ )
! 1555: c[i] = w[i]?FTOIF(w[i]):0;
! 1556: else
! 1557: for ( i = 0; i < n; i++ )
! 1558: c[i] = w[i]?FTOIF(w[i]-1):0;
! 1559: if ( !c[0] )
! 1560: continue;
! 1561: diffsfum(r,dr); cpyum(r,t); gcdsfum(t,dr,g);
! 1562: if ( DEG(g) > 0 )
! 1563: continue;
! 1564:
! 1565: czsfum(r,f);
! 1566: for ( i = 0; f[i]; i++ );
! 1567: if ( i == 1 )
! 1568: return;
! 1569: }
! 1570: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to Hgfs.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib2 / asir2018 / engine