Annotation of OpenXM_contrib2/asir2018/engine/E.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:
! 52: void henmv(VL vl,VN vn,P f,P g0,P h0,P a0,P b0,P lg,P lh,P lg0,P lh0,Z q,int k,P *gp,P *hp)
! 53: {
! 54: P g1,h1,a1,b1;
! 55: Z q2,r2,two;
! 56:
! 57: STOQ(2,two);
! 58: divqrz(q,two,&q2,&r2);
! 59: adjc(vl,g0,a0,lg0,q,&g1,&a1); adjc(vl,h0,b0,lh0,q,&h1,&b1);
! 60: henmvmain(vl,vn,f,g1,h1,b1,a1,lg,lh,q,q2,k,gp,hp);
! 61: }
! 62:
! 63: void henmvmain(VL vl,VN vn,P f,P fi0,P fi1,P gi0,P gi1,P l0,P l1,Z mod,Z mod2,int k,P *fr0,P *fr1)
! 64: {
! 65: V v;
! 66: int n,i,j;
! 67: int *md;
! 68: P x,m,m1,c,q,r,a,s,u,ff,f0,f1;
! 69: P w0,w1,cf,cfi,t,q1,dvr;
! 70: P *c0,*c1;
! 71: P *f0h,*f1h;
! 72:
! 73: v = VR(f); n = deg(v,f); MKV(v,x);
! 74: c0 = (P *)ALLOCA((n+1)*sizeof(P));
! 75: c1 = (P *)ALLOCA((n+1)*sizeof(P));
! 76: invz((Z)LC(fi1),mod,(Z *)&c); mulp(vl,fi1,c,&t); cmp(mod,t,&dvr);
! 77: cm2p(mod,mod2,gi0,&c0[0]); cm2p(mod,mod2,gi1,&c1[0]);
! 78: for ( i = 1; i <= n; i++ ) {
! 79: mulp(vl,x,c1[i-1],&m); divsrp(vl,m,dvr,&q,&r); mulp(vl,q,c,&q1);
! 80: cm2p(mod,mod2,r,&c1[i]);
! 81: mulp(vl,x,c0[i-1],&m); mulp(vl,q1,fi0,&m1); addp(vl,m,m1,&a);
! 82: cm2p(mod,mod2,a,&c0[i]);
! 83: }
! 84: affine(vl,f,vn,&t); cm2p(mod,mod2,t,&ff);
! 85: for ( i = 0; vn[i].v; i++ );
! 86: md = ( int *) ALLOCA((i+1)*sizeof(int));
! 87: for ( i = 0; vn[i].v; i++ )
! 88: md[i] = getdeg(vn[i].v,ff);
! 89: cm2p(mod,mod2,fi0,&f0); affine(vl,l0,vn,&t);
! 90: if ( NUM(f0) )
! 91: cm2p(mod,mod2,t,&f0);
! 92: else
! 93: cm2p(mod,mod2,t,&COEF(DC(f0)));
! 94: cm2p(mod,mod2,fi1,&f1); affine(vl,l1,vn,&t);
! 95: if ( NUM(f1) )
! 96: cm2p(mod,mod2,t,&f1);
! 97: else
! 98: cm2p(mod,mod2,t,&COEF(DC(f1)));
! 99: W_CALLOC(k,P,f0h); W_CALLOC(k,P,f1h);
! 100: for ( i = 0; i <= k; i++ ) {
! 101: exthpc(vl,v,f0,i,&f0h[i]); exthpc(vl,v,f1,i,&f1h[i]);
! 102: }
! 103: for ( j = 1; j <= k; j++ ) {
! 104: for ( i = 0; vn[i].v; i++ )
! 105: if ( getdeg(vn[i].v,f0)+getdeg(vn[i].v,f1) > md[i] )
! 106: goto END;
! 107: for ( i = 0, s = 0; i <= j; i++ ) {
! 108: mulp(vl,f0h[i],f1h[j-i],&t); addp(vl,s,t,&u); s = u;
! 109: }
! 110: cm2p(mod,mod2,s,&t);
! 111: exthpc(vl,v,ff,j,&u); subm2p(vl,mod,mod2,u,t,&cf);
! 112: for ( i = 0, w0 = 0, w1 = 0; i <= n; i++ ) {
! 113: if ( !cf )
! 114: cfi = 0;
! 115: else if ( VR(cf) == v )
! 116: coefp(cf,i,&cfi);
! 117: else if ( i == 0 )
! 118: cfi = cf;
! 119: else
! 120: cfi = 0;
! 121: if ( cfi ) {
! 122: mulp(vl,cfi,c0[i],&m); addp(vl,w0,m,&a); w0 = a;
! 123: mulp(vl,cfi,c1[i],&m); addp(vl,w1,m,&a); w1 = a;
! 124: }
! 125: }
! 126: cm2p(mod,mod2,w0,&t); addm2p(vl,mod,mod2,f0,t,&a);
! 127: addm2p(vl,mod,mod2,f0h[j],t,&s); f0h[j] = s; f0 = a;
! 128: cm2p(mod,mod2,w1,&u); addm2p(vl,mod,mod2,f1,u,&a);
! 129: addm2p(vl,mod,mod2,f1h[j],u,&s); f1h[j] = s; f1 = a;
! 130: if ( !t ) {
! 131: restore(vl,f0,vn,&s); cm2p(mod,mod2,s,&t);
! 132: if ( divtpz(vl,f,t,&s) ) {
! 133: *fr0 = t; *fr1 = s;
! 134: return;
! 135: }
! 136: }
! 137: if ( !u ) {
! 138: restore(vl,f1,vn,&s); cm2p(mod,mod2,s,&t);
! 139: if ( divtpz(vl,f,t,&s) ) {
! 140: *fr0 = s; *fr1 = t;
! 141: return;
! 142: }
! 143: }
! 144: }
! 145: END:
! 146: restore(vl,f0,vn,&t); cm2p(mod,mod2,t,fr0);
! 147: restore(vl,f1,vn,&t); cm2p(mod,mod2,t,fr1);
! 148: }
! 149:
! 150: /*
! 151: input : f, fi0, fi1, p, k; f = fi0 * fi1 mod p; ( p, k : integer )
! 152: output : gr1 * fr0 + gr0 * fr1 = 1 mod qr; ( qr = p^(2^(k-1)) )
! 153: */
! 154:
! 155: void henzq(P f,P i0,UM fi0,P i1,UM fi1,int p,int k,P *fr0p,P *fr1p,P *gr0p,P *gr1p,Z *qrp)
! 156: {
! 157: Z q,qq,q2,two,r2;
! 158: int n,i;
! 159: UM wg0,wg1,wf0,wf1;
! 160: P f0,f1,g0,g1,m,m1,s,c,d,q1,r1,rm,rm1,a,a1,a2;
! 161:
! 162: n = UDEG(f);
! 163: wf0 = W_UMALLOC(n); wf1 = W_UMALLOC(n);
! 164: wg0 = W_UMALLOC(n); wg1 = W_UMALLOC(n);
! 165: cpyum(fi0,wf0); cpyum(fi1,wf1);
! 166: eucum(p,wf0,wf1,wg1,wg0);
! 167: umtop(VR(f),fi0,&f0); umtop(VR(f),fi1,&f1);
! 168: umtop(VR(f),wg0,&g0); umtop(VR(f),wg1,&g1);
! 169:
! 170: STOQ(2,two); STOQ(p,q); divqrz(q,two,&q2,&r2);
! 171: for ( i = 1; i < k; i++ ) {
! 172: /* c = ((f - f0*f1)/q) mod q;
! 173: q1 = (c*g1) / f1;
! 174: r1 = (c*g1) mod f1;
! 175: f1 += (r1 mod q)*q;
! 176: f0 += ((c*g0 + q1*f0) mod q)*q;
! 177:
! 178: d = ((1 - (f1*g0 + f0*g1))/q) mod q;
! 179: q1 = (d*g0) / f1;
! 180: r1 = (d*g0) mod f1;
! 181: g1 += (r1 mod q)*q;
! 182: g0 += ((c*g0 + q1*f0) mod q)*q;
! 183: q = q^2;
! 184: */
! 185:
! 186: /* c = ((f - f0*f1)/q) mod q */
! 187: mulp(CO,f0,f1,&m); subp(CO,f,m,&s);
! 188: divcp(s,(Q)q,&m); cm2p(q,q2,m,&c);
! 189:
! 190: /* q1 = (c*g1) / f1; r1 = (c*g1) mod f1; */
! 191: mulp(CO,c,g1,&m); cm2p(q,q2,m,&s);
! 192: udivpwm(q,s,f1,&q1,&r1);
! 193:
! 194: /* f1 = f1 + (r1 mod q)*q; */
! 195: cm2p(q,q2,r1,&rm);
! 196: mulpq(rm,(P)q,&rm1); addp(CO,f1,rm1,&a);
! 197: f1 = a;
! 198:
! 199: /* a1 = (c*g0 + q1*f0) mod q; */
! 200: mulp(CO,c,g0,&m); mulp(CO,q1,f0,&m1); addp(CO,m,m1,&a);
! 201: cm2p(q,q2,a,&a1);
! 202:
! 203: /* f0 = f0 + a1*q; */
! 204: mulpq(a1,(P)q,&a2);
! 205: addp(CO,f0,a2,&a);
! 206: f0 = a;
! 207:
! 208: /* d = ((1 - (f1*g0 + f0*g1))/q) mod q; */
! 209: mulp(CO,f1,g0,&m); mulp(CO,f0,g1,&m1); addp(CO,m,m1,&a);
! 210: subp(CO,(P)ONE,a,&s);
! 211: divcp(s,(Q)q,&m); cm2p(q,q2,m,&d);
! 212:
! 213: /* q1 = (d*g0) / f1; r1 = (d*g0) mod f1; */
! 214: mulp(CO,d,g1,&m); udivpwm(q,m,f1,&q1,&r1);
! 215:
! 216: /* g1 = g1 + (r1 mod q )*q; */
! 217: cm2p(q,q2,r1,&rm);
! 218: mulpq(rm,(P)q,&rm1); addp(CO,g1,rm1,&a);
! 219: g1 = a;
! 220:
! 221: /* g0 = g0 + (c*g0 + q1*f0 mod q)*q; */
! 222: mulp(CO,d,g0,&m); mulp(CO,q1,f0,&m1); addp(CO,m,m1,&a);
! 223: cm2p(q,q2,a,&a1); mulpq(a1,(P)q,&a2);
! 224: addp(CO,g0,a2,&a);
! 225: g0 = a;
! 226:
! 227: /* q = q^2; */
! 228: mulz(q,q,&qq); q = qq; divqrz(q,two,&q2,&r2);
! 229: }
! 230: *fr0p = f0; *fr1p = f1; *gr0p = g0; *gr1p = g1; *qrp = q;
! 231: }
! 232:
! 233: void henzq1(P g,P h,Q bound,P *gcp,P *hcp,Z *qp)
! 234: {
! 235: V v;
! 236: Z w,q,q1;
! 237: Q f;
! 238: Q u,t,a,b,s,ts;
! 239: P c,c1;
! 240: P tg,th,tr;
! 241: UM wg,wh,ws,wt,wm;
! 242: int n,m,modres,mod,index,i;
! 243: P gc0,hc0;
! 244: P z,zz,zzz;
! 245:
! 246:
! 247: v = VR(g); n=deg(v,g); m=deg(v,h);
! 248: norm(g,&a); norm(h,&b);
! 249: STOQ(m,w); pwrq(a,(Q)w,&t);
! 250: STOQ(n,w); pwrq(b,(Q)w,&s);
! 251: mulq(t,s,&ts);
! 252:
! 253: factorialz(n+m,&w); mulq(ts,(Q)w,&s);
! 254: addq(s,s,&f);
! 255:
! 256: wg = W_UMALLOC(m+n); wh = W_UMALLOC(m+n);
! 257: wt = W_UMALLOC(m+n); ws = W_UMALLOC(m+n);
! 258: wm = W_UMALLOC(m+n);
! 259:
! 260: for ( q = ONE, t = 0, c = 0, index = 0; ; ) {
! 261: mod = get_lprime(index++);
! 262: if ( !remqi((Q)LC(g),mod) || !remqi((Q)LC(h),mod) )
! 263: continue;
! 264: ptomp(mod,g,&tg); ptomp(mod,h,&th);
! 265: srchump(mod,tg,th,&tr);
! 266: if ( !tr )
! 267: continue;
! 268: else
! 269: modres = CONT((MQ)tr);
! 270:
! 271: mptoum(tg,wg); mptoum(th,wh);
! 272: eucum(mod,wg,wh,ws,wt); /* wg * ws + wh * wt = 1 */
! 273: DEG(wm) = 0; COEF(wm)[0] = modres;
! 274: mulum(mod,ws,wm,wt);
! 275: for ( i = DEG(wt); i >= 0; i-- )
! 276: if ( ( COEF(wt)[i] * 2 ) > mod )
! 277: COEF(wt)[i] -= mod;
! 278: chnrem(mod,v,c,q,wt,&c1,&q1);
! 279: if ( !ucmpp(c,c1) ) {
! 280: mulp(CO,c,g,&z); divsrp(CO,z,h,&zz,&zzz);
! 281: if ( NUM(zzz) ) {
! 282: q = q1; c = c1;
! 283: break;
! 284: }
! 285: }
! 286: q = q1; c = c1;
! 287:
! 288: if ( cmpq(f,(Q)q) < 0 )
! 289: break;
! 290: }
! 291: ptozp(c,1,&s,&gc0);
! 292: /* gc0*g = zz*h + zzz -> c*g + (-zz)*h = zzz */
! 293: mulp(CO,gc0,g,&z);
! 294: divsrp(CO,z,h,&zz,&zzz);
! 295: ptozp(zz,1,&s,(P *)&t);
! 296: if ( INT(s) )
! 297: chsgnp(zz,&hc0);
! 298: else {
! 299: MPZTOZ(mpq_denref(BDY(s)),q); mulq((Q)q,(Q)zzz,(Q *)&q1); zzz = (P)q1;
! 300: mulp(CO,gc0,(P)q,&c); gc0 = c; mulp(CO,zz,(P)q,&c); chsgnp(c,&hc0);
! 301: }
! 302: if ( !INT((Q)zzz) ) {
! 303: MPZTOZ(mpq_denref(BDY((Q)zzz)),q); MPZTOZ(mpq_numref(BDY((Q)zzz)),q1); zzz = (P)q1;
! 304: mulp(CO,gc0,(P)q,&c); gc0 = c; mulp(CO,hc0,(P)q,&c); hc0 = c;
! 305: }
! 306: for ( index = 0; ; ) {
! 307: mod = get_lprime(index++);
! 308: if ( !remqi((Q)zzz,mod) ||
! 309: !remqi((Q)LC(g),mod) ||
! 310: !remqi((Q)LC(h),mod) )
! 311: continue;
! 312: for ( STOQ(mod,q); cmpq((Q)q,bound) < 0; ) {
! 313: mulz(q,q,&q1); q = q1;
! 314: }
! 315: *qp = q;
! 316: invz((Z)zzz,q,&q1);
! 317: mulp(CO,gc0,(P)q1,hcp); mulp(CO,hc0,(P)q1,gcp);
! 318: return;
! 319: }
! 320: }
! 321:
! 322: void addm2p(VL vl,Z mod,Z mod2,P n1,P n2,P *nr)
! 323: {
! 324: P t;
! 325:
! 326: addp(vl,n1,n2,&t);
! 327: if ( !t )
! 328: *nr = 0;
! 329: else
! 330: cm2p(mod,mod2,t,nr);
! 331: }
! 332:
! 333: void subm2p(VL vl,Z mod,Z mod2,P n1,P n2,P *nr)
! 334: {
! 335: P t;
! 336:
! 337: subp(vl,n1,n2,&t);
! 338: if ( !t )
! 339: *nr = 0;
! 340: else
! 341: cm2p(mod,mod2,t,nr);
! 342: }
! 343:
! 344: void mulm2p(VL vl,Z mod,Z mod2,P n1,P n2,P *nr)
! 345: {
! 346: P t;
! 347:
! 348: mulp(vl,n1,n2,&t);
! 349: if ( !t )
! 350: *nr = 0;
! 351: else
! 352: cm2p(mod,mod2,t,nr);
! 353: }
! 354:
! 355: void cmp(Z mod,P p,P *pr)
! 356: {
! 357: P t;
! 358: DCP dc,dcr,dcr0;
! 359:
! 360: if ( !p )
! 361: *pr = 0;
! 362: else if ( NUM(p) )
! 363: remz((Z)p,mod,(Z *)pr);
! 364: else {
! 365: for ( dc = DC(p), dcr0 = 0; dc; dc = NEXT(dc) ) {
! 366: cmp(mod,COEF(dc),&t);
! 367: if ( t ) {
! 368: NEXTDC(dcr0,dcr);
! 369: DEG(dcr) = DEG(dc);
! 370: COEF(dcr) = t;
! 371: }
! 372: }
! 373: if ( !dcr0 )
! 374: *pr = 0;
! 375: else {
! 376: NEXT(dcr) = 0;
! 377: MKP(VR(p),dcr0,*pr);
! 378: }
! 379: }
! 380: }
! 381:
! 382: void cm2p(Z mod,Z m,P p,P *pr)
! 383: {
! 384: P t;
! 385: DCP dc,dcr,dcr0;
! 386:
! 387: if ( !p )
! 388: *pr = 0;
! 389: else if ( NUM(p) )
! 390: rem2q((Z)p,mod,m,(Z *)pr);
! 391: else {
! 392: for ( dc = DC(p), dcr0 = 0; dc; dc = NEXT(dc) ) {
! 393: cm2p(mod,m,COEF(dc),&t);
! 394: if ( t ) {
! 395: NEXTDC(dcr0,dcr);
! 396: DEG(dcr) = DEG(dc);
! 397: COEF(dcr) = t;
! 398: }
! 399: }
! 400: if ( !dcr0 )
! 401: *pr = 0;
! 402: else {
! 403: NEXT(dcr) = 0;
! 404: MKP(VR(p),dcr0,*pr);
! 405: }
! 406: }
! 407: }
! 408:
! 409: void addm2q(Z mod,Z mod2,Z n1,Z n2,Z *nr)
! 410: {
! 411: Z t;
! 412:
! 413: addz(n1,n2,&t);
! 414: if ( !t )
! 415: *nr = 0;
! 416: else
! 417: rem2q(t,mod,mod2,nr);
! 418: }
! 419:
! 420: void subm2q(Z mod,Z mod2,Z n1,Z n2,Z *nr)
! 421: {
! 422: Z t;
! 423:
! 424: subz(n1,n2,&t);
! 425: if ( !t )
! 426: *nr = 0;
! 427: else
! 428: rem2q(t,mod,mod2,nr);
! 429: }
! 430:
! 431: void mulm2q(Z mod,Z mod2,Z n1,Z n2,Z *nr)
! 432: {
! 433: Z t;
! 434:
! 435: mulz(n1,n2,&t);
! 436: if ( !t )
! 437: *nr = 0;
! 438: else
! 439: rem2q(t,mod,mod2,nr);
! 440: }
! 441:
! 442: void rem2q(Z n,Z m,Z m2,Z *nr)
! 443: {
! 444: Z r;
! 445: int sgn;
! 446:
! 447: remz(n,m,&r);
! 448: if ( !r )
! 449: *nr = 0;
! 450: else {
! 451: sgn = cmpz(r,m2);
! 452: if ( sgn > 0 )
! 453: subz(r,m,nr);
! 454: else
! 455: *nr = r;
! 456: }
! 457: }
! 458:
! 459: /*
! 460: extract d-homogeneous part with respect to vl - {v}
! 461: */
! 462:
! 463: void exthpc_generic(VL vl,P p,int d,V v,P *pr)
! 464: {
! 465: P w,x,t,t1,a,xd;
! 466: V v0;
! 467: DCP dc;
! 468:
! 469: if ( d < 0 || !p )
! 470: *pr = 0;
! 471: else if ( NUM(p) )
! 472: if ( d == 0 )
! 473: *pr = p;
! 474: else
! 475: *pr = 0;
! 476: else if ( v == VR(p) )
! 477: exthpc(vl,v,p,d,pr);
! 478: else {
! 479: v0 = VR(p);
! 480: for ( MKV(v0,x), dc = DC(p), w = 0; dc; dc = NEXT(dc) ) {
! 481: exthpc_generic(vl,COEF(dc),d-QTOS(DEG(dc)),v,&t);
! 482: pwrp(vl,x,DEG(dc),&xd);
! 483: mulp(vl,xd,t,&t1); addp(vl,w,t1,&a); w = a;
! 484: }
! 485: *pr = w;
! 486: }
! 487: }
! 488:
! 489: void exthp(VL vl,P p,int d,P *pr)
! 490: {
! 491: P t,t1,a,w,x,xd;
! 492: DCP dc;
! 493:
! 494: if ( d < 0 )
! 495: *pr = 0;
! 496: else if ( NUM(p) )
! 497: if ( d == 0 )
! 498: *pr = p;
! 499: else
! 500: *pr = 0;
! 501: else {
! 502: for ( MKV(VR(p),x), dc = DC(p), w = 0; dc; dc = NEXT(dc) ) {
! 503: exthp(vl,COEF(dc),d - QTOS(DEG(dc)),&t);
! 504: pwrp(vl,x,DEG(dc),&xd);
! 505: mulp(vl,xd,t,&t1); addp(vl,w,t1,&a); w = a;
! 506: }
! 507: *pr = w;
! 508: }
! 509: }
! 510:
! 511: void exthpc(VL vl,V v,P p,int d,P *pr)
! 512: {
! 513: P t,t1,a,w,x,xd;
! 514: DCP dc;
! 515:
! 516: if ( v != VR(p) )
! 517: exthp(vl,p,d,pr);
! 518: else if ( d < 0 )
! 519: *pr = 0;
! 520: else {
! 521: for ( MKV(v,x), dc = DC(p), w = 0; dc; dc = NEXT(dc) ) {
! 522: exthp(vl,COEF(dc),d,&t);
! 523: pwrp(vl,x,DEG(dc),&xd);
! 524: mulp(vl,xd,t,&t1); addp(vl,w,t1,&a); w = a;
! 525: }
! 526: *pr = w;
! 527: }
! 528: }
! 529:
! 530: void cbound(VL vl,P p,Q *b)
! 531: {
! 532: Q n,m;
! 533: Z t,e;
! 534: int k;
! 535:
! 536: cmax(p,&n);
! 537: addq(n,n,&m);
! 538:
! 539: k = geldb(vl,p);
! 540: STOQ(3,t); STOQ(k,e);
! 541:
! 542: pwrq((Q)t,(Q)e,&n);
! 543: mulq(m,n,b);
! 544: }
! 545:
! 546: int geldb(VL vl,P p)
! 547: {
! 548: int m;
! 549:
! 550: for ( m = 0; vl; vl = NEXT(vl) )
! 551: m += getdeg(vl->v,p);
! 552: return ( m );
! 553: }
! 554:
! 555: int getdeg(V v,P p)
! 556: {
! 557: int m,t;
! 558: DCP dc;
! 559:
! 560: if ( !p || NUM(p) )
! 561: return ( 0 );
! 562: else if ( v == VR(p) )
! 563: return ( deg(v,p) );
! 564: else {
! 565: for ( dc = DC(p), m = 0; dc; dc = NEXT(dc) ) {
! 566: t = getdeg(v,COEF(dc));
! 567: m = MAX(m,t);
! 568: }
! 569: return ( m );
! 570: }
! 571: }
! 572:
! 573: /* YYY */
! 574:
! 575: void cmax(P p,Q *b)
! 576: {
! 577: DCP dc;
! 578: Q m,m1;
! 579:
! 580: if ( NUM(p) )
! 581: absq((Q)p,b);
! 582: else {
! 583: for ( dc = DC(p), m = 0; dc; dc = NEXT(dc) ) {
! 584: cmax(COEF(dc),&m1);
! 585: if ( cmpq(m1,m) > 0 )
! 586: m = m1;
! 587: }
! 588: *b = m;
! 589: }
! 590: }
! 591:
! 592: int nextbin(VN vn,int n)
! 593: {
! 594: int tmp,i,carry;
! 595:
! 596: if ( n == 0 )
! 597: return ( 1 );
! 598:
! 599: for ( i = n - 1, carry = 1; i >= 0; i-- ) {
! 600: tmp = vn[i].n + carry;
! 601: vn[i].n = tmp % 2;
! 602: carry = tmp / 2;
! 603: }
! 604: return ( carry );
! 605: }
! 606:
! 607: void mulsgn(VN vn,VN vnt,int n,VN vn1)
! 608: {
! 609: int i;
! 610:
! 611: for ( i = 0; vn[i].v; i++ )
! 612: vn1[i].n = vn[i].n;
! 613: for ( i = 0; i < n; i++ )
! 614: if ( vnt[i].n )
! 615: vn1[(long)vnt[i].v].n *= -1;
! 616: }
! 617:
! 618: void next(VN vn)
! 619: {
! 620: int i,m,n,tmp,carry;
! 621:
! 622: for ( m = 0, i = 0; vn[i].v; i++ )
! 623: m = MAX(m,ABS(vn[i].n));
! 624: if ( m == 0 ) {
! 625: vn[--i].n = 1;
! 626: return;
! 627: }
! 628: for ( carry = 1, m++, i--; ( i >= 0 ) && carry; i-- ) {
! 629: tmp = vn[i].n + carry;
! 630: vn[i].n = tmp % m;
! 631: carry = tmp / m;
! 632: }
! 633: if ( ( i == -1 ) && carry ) {
! 634: for ( i = 0; vn[i].v; i++ )
! 635: vn[i].n = 0;
! 636: vn[--i].n = m;
! 637: } else {
! 638: for ( n = 0, i = 0; vn[i].v; i++ )
! 639: n = MAX(n,ABS(vn[i].n));
! 640: if ( n < m - 1 )
! 641: vn[--i].n = m - 1;
! 642: }
! 643: }
! 644:
! 645: void clctv(VL vl,P p,VL *nvlp)
! 646: {
! 647: int i,n;
! 648: VL tvl;
! 649: VN tvn;
! 650:
! 651: if ( !p || NUM(p) ) {
! 652: *nvlp = 0;
! 653: return;
! 654: }
! 655:
! 656: for ( n = 0, tvl = vl; tvl; tvl = NEXT(tvl), n++ );
! 657: tvn = (VN) ALLOCA((n+1)*sizeof(struct oVN));
! 658: for ( i = 0, tvl = vl; tvl; tvl = NEXT(tvl), i++ ) {
! 659: tvn[i].v = tvl->v;
! 660: tvn[i].n = 0;
! 661: }
! 662:
! 663: markv(tvn,n,p);
! 664: vntovl(tvn,n,nvlp);
! 665: }
! 666:
! 667: void markv(VN vn,int n,P p)
! 668: {
! 669: V v;
! 670: DCP dc;
! 671: int i;
! 672:
! 673: if ( NUM(p) )
! 674: return;
! 675: v = VR(p);
! 676: for ( i = 0, v = VR(p); i < n; i++ )
! 677: if ( v == vn[i].v ) {
! 678: vn[i].n = 1;
! 679: break;
! 680: }
! 681: for ( dc = DC(p); dc; dc = NEXT(dc) )
! 682: markv(vn,n,COEF(dc));
! 683: }
! 684:
! 685: void vntovl(VN vn,int n,VL *vlp)
! 686: {
! 687: int i;
! 688: VL tvl,tvl0;
! 689:
! 690: for ( i = 0, tvl0 = 0; ; ) {
! 691: while ( ( i < n ) && ( vn[i].n == 0 ) ) i++;
! 692: if ( i == n )
! 693: break;
! 694: else {
! 695: if ( !tvl0 ) {
! 696: NEWVL(tvl0);
! 697: tvl = tvl0;
! 698: } else {
! 699: NEWVL(NEXT(tvl));
! 700: tvl = NEXT(tvl);
! 701: }
! 702: tvl->v = vn[i++].v;
! 703: }
! 704: }
! 705: if ( tvl0 )
! 706: NEXT(tvl) = 0;
! 707: *vlp = tvl0;
! 708: }
! 709:
! 710: int dbound(V v,P f)
! 711: {
! 712: int m,t;
! 713: DCP dc;
! 714:
! 715: if ( !f )
! 716: return ( -1 );
! 717: else if ( v != VR(f) )
! 718: return homdeg(f);
! 719: else {
! 720: for ( dc = DC(f), m = 0; dc; dc = NEXT(dc) ) {
! 721: t = homdeg(COEF(dc));
! 722: m = MAX(m,t);
! 723: }
! 724: return ( m );
! 725: }
! 726: }
! 727:
! 728: int homdeg(P f)
! 729: {
! 730: int m,t;
! 731: DCP dc;
! 732:
! 733: if ( !f )
! 734: return ( -1 );
! 735: else if ( NUM(f) )
! 736: return( 0 );
! 737: else {
! 738: for ( dc = DC(f), m = 0; dc; dc = NEXT(dc) ) {
! 739: t = QTOS(DEG(dc))+homdeg(COEF(dc));
! 740: m = MAX(m,t);
! 741: }
! 742: return ( m );
! 743: }
! 744: }
! 745:
! 746: int minhomdeg(P f)
! 747: {
! 748: int m,t;
! 749: DCP dc;
! 750:
! 751: if ( !f )
! 752: return ( -1 );
! 753: else if ( NUM(f) )
! 754: return( 0 );
! 755: else {
! 756: for ( dc = DC(f), m = homdeg(f); dc; dc = NEXT(dc) ) {
! 757: t = QTOS(DEG(dc))+minhomdeg(COEF(dc));
! 758: m = MIN(m,t);
! 759: }
! 760: return ( m );
! 761: }
! 762: }
! 763:
! 764: void adjc(VL vl,P f,P a,P lc0,Z q,P *fr,P *ar)
! 765: {
! 766: Z m,t;
! 767: P m1;
! 768:
! 769: invz((Z)LC(f),q,&t);
! 770: mulz((Z)lc0,t,&m);
! 771: mulpq(f,(P)m,&m1); cmp(q,m1,fr);
! 772: invz(m,q,&t);
! 773: mulpq(a,(P)t,&m1);
! 774: cmp(q,m1,ar);
! 775: }
! 776:
! 777: void affinemain(VL vl,P p,V v0,int n,P *pl,P *pr)
! 778: {
! 779: P x,t,m,c,s,a;
! 780: DCP dc;
! 781: Z d;
! 782:
! 783: if ( !p )
! 784: *pr = 0;
! 785: else if ( NUM(p) )
! 786: *pr = p;
! 787: else if ( VR(p) != v0 ) {
! 788: MKV(VR(p),x);
! 789: for ( c = 0, dc = DC(p); dc; dc = NEXT(dc) ) {
! 790: affinemain(vl,COEF(dc),v0,n,pl,&t);
! 791: if ( DEG(dc) ) {
! 792: pwrp(vl,x,DEG(dc),&s); mulp(vl,s,t,&m);
! 793: addp(vl,m,c,&a); c = a;
! 794: } else {
! 795: addp(vl,t,c,&a); c = a;
! 796: }
! 797: }
! 798: *pr = c;
! 799: } else {
! 800: dc = DC(p);
! 801: c = COEF(dc);
! 802: for ( d = DEG(dc), dc = NEXT(dc);
! 803: dc; d = DEG(dc), dc = NEXT(dc) ) {
! 804: mulp(vl,pl[QTOS(d)-QTOS(DEG(dc))],c,&m);
! 805: addp(vl,m,COEF(dc),&c);
! 806: }
! 807: if ( d ) {
! 808: mulp(vl,pl[QTOS(d)],c,&m); c = m;
! 809: }
! 810: *pr = c;
! 811: }
! 812: }
! 813:
! 814: void restore(VL vl,P f,VN vn,P *fr)
! 815: {
! 816: int i;
! 817: P vv,g,g1,t;
! 818: Z s;
! 819:
! 820: g = f;
! 821: for ( i = 0; vn[i].v; i++ ) {
! 822: MKV(vn[i].v,t);
! 823: if ( vn[i].n ) {
! 824: STOQ(-vn[i].n,s);
! 825: addp(vl,t,(P)s,&vv);
! 826: } else
! 827: vv = t;
! 828:
! 829: substp(vl,g,vn[i].v,vv,&g1); g = g1;
! 830: }
! 831: *fr = g;
! 832: }
! 833:
! 834: void mergev(VL vl,VL vl1,VL vl2,VL *nvlp)
! 835: {
! 836: int i,n;
! 837: VL tvl;
! 838: VN vn;
! 839:
! 840: if ( !vl1 ) {
! 841: *nvlp = vl2; return;
! 842: } else if ( !vl2 ) {
! 843: *nvlp = vl1; return;
! 844: }
! 845: for ( i = 0, tvl = vl; tvl; i++, tvl = NEXT(tvl) );
! 846: n = i;
! 847: W_CALLOC(n,struct oVN,vn);
! 848: for ( i = 0, tvl = vl; tvl; tvl = NEXT(tvl), i++ )
! 849: vn[i].v = tvl->v;
! 850: for ( i = 0, tvl = vl1; (i < n) && tvl; i++, tvl = NEXT(tvl)) {
! 851: while ( ( i < n ) && ( vn[i].v != tvl->v ) )
! 852: i++;
! 853: if ( i == n )
! 854: break;
! 855: else
! 856: vn[i].n = 1;
! 857: }
! 858: for ( i = 0, tvl = vl2; (i < n) && tvl; i++, tvl = NEXT(tvl) ) {
! 859: while ( ( i < n ) && ( vn[i].v != tvl->v ) )
! 860: i++;
! 861: if ( i == n )
! 862: break;
! 863: else
! 864: vn[i].n = 1;
! 865: }
! 866: vntovl(vn,n,nvlp);
! 867: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to E.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib2 / asir2018 / engine