Annotation of OpenXM_contrib2/asir2000/engine/E.c, Revision 1.9
1.2 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
1.3 noro 26: * e-mail at risa-admin@sec.flab.fujitsu.co.jp of the detailed specification
1.2 noro 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: *
1.9 ! noro 48: * $OpenXM: OpenXM_contrib2/asir2000/engine/E.c,v 1.8 2001/10/09 01:36:09 noro Exp $
1.2 noro 49: */
1.1 noro 50: #include "ca.h"
51:
1.8 noro 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,Q q,int k,P *gp,P *hp)
1.1 noro 53: {
1.9 ! noro 54: P g1,h1,a1,b1;
! 55: N qn;
! 56: Q q2;
! 57:
! 58: divin((NM(q)),2,&qn); NTOQ(qn,1,q2);
! 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);
1.1 noro 61: }
62:
1.8 noro 63: void henmvmain(VL vl,VN vn,P f,P fi0,P fi1,P gi0,P gi1,P l0,P l1,Q mod,Q mod2,int k,P *fr0,P *fr1)
1.1 noro 64: {
1.9 ! noro 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: invl((Q)LC(fi1),mod,(Q *)&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: }
1.1 noro 145: END:
1.9 ! noro 146: restore(vl,f0,vn,&t); cm2p(mod,mod2,t,fr0);
! 147: restore(vl,f1,vn,&t); cm2p(mod,mod2,t,fr1);
1.1 noro 148: }
149:
150: /*
1.9 ! noro 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)) )
1.1 noro 153: */
154:
1.8 noro 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,Q *qrp)
1.1 noro 156: {
1.9 ! noro 157: N qn;
! 158: Q q,qq,q2;
! 159: int n,i;
! 160: UM wg0,wg1,wf0,wf1;
! 161: P f0,f1,g0,g1,m,m1,s,c,d,q1,r1,rm,rm1,a,a1,a2;
! 162:
! 163: n = UDEG(f);
! 164: wf0 = W_UMALLOC(n); wf1 = W_UMALLOC(n);
! 165: wg0 = W_UMALLOC(n); wg1 = W_UMALLOC(n);
! 166: cpyum(fi0,wf0); cpyum(fi1,wf1);
! 167: eucum(p,wf0,wf1,wg1,wg0);
! 168: umtop(VR(f),fi0,&f0); umtop(VR(f),fi1,&f1);
! 169: umtop(VR(f),wg0,&g0); umtop(VR(f),wg1,&g1);
1.1 noro 170:
1.9 ! noro 171: STOQ(p,q); divin(NM(q),2,&qn); NTOQ(qn,1,q2);
! 172: for ( i = 1; i < k; i++ ) {
1.1 noro 173: #if 0
1.9 ! noro 174: mulp(CO,i1,g0,&m); mulp(CO,i0,g1,&m1); addp(CO,m,m1,&a);
! 175: if ( NUM(a) ) {
! 176: for ( STOQ(p,q), j = 1; j < k; j++ ) {
! 177: mulq(q,q,&qq); q = qq;
! 178: }
! 179: f0 = i0; f1 = i1;
! 180: invl(a,q,&qq);
! 181: mulp(CO,g0,qq,&s); g0 = s; mulp(CO,g1,qq,&s); g1 = s;
! 182: break;
! 183: }
1.1 noro 184: #endif
1.9 ! noro 185: /* c = ((f - f0*f1)/q) mod q;
! 186: q1 = (c*g1) / f1;
! 187: r1 = (c*g1) mod f1;
! 188: f1 += (r1 mod q)*q;
! 189: f0 += ((c*g0 + q1*f0) mod q)*q;
! 190:
! 191: d = ((1 - (f1*g0 + f0*g1))/q) mod q;
! 192: q1 = (d*g0) / f1;
! 193: r1 = (d*g0) mod f1;
! 194: g1 += (r1 mod q)*q;
! 195: g0 += ((c*g0 + q1*f0) mod q)*q;
! 196: q = q^2;
! 197: */
! 198:
! 199: /* c = ((f - f0*f1)/q) mod q */
! 200: mulp(CO,f0,f1,&m); subp(CO,f,m,&s);
! 201: divcp(s,q,&m); cm2p(q,q2,m,&c);
! 202:
! 203: /* q1 = (c*g1) / f1; r1 = (c*g1) mod f1; */
! 204: mulp(CO,c,g1,&m); cm2p(q,q2,m,&s);
! 205: udivpwm(q,s,f1,&q1,&r1);
! 206:
! 207: /* f1 = f1 + (r1 mod q)*q; */
! 208: cm2p(q,q2,r1,&rm);
! 209: mulpq(rm,(P)q,&rm1); addp(CO,f1,rm1,&a);
! 210: f1 = a;
! 211:
! 212: /* a1 = (c*g0 + q1*f0) mod q; */
! 213: mulp(CO,c,g0,&m); mulp(CO,q1,f0,&m1); addp(CO,m,m1,&a);
! 214: cm2p(q,q2,a,&a1);
! 215:
! 216: /* f0 = f0 + a1*q; */
! 217: mulpq(a1,(P)q,&a2);
! 218: addp(CO,f0,a2,&a);
! 219: f0 = a;
! 220:
! 221: /* d = ((1 - (f1*g0 + f0*g1))/q) mod q; */
! 222: mulp(CO,f1,g0,&m); mulp(CO,f0,g1,&m1); addp(CO,m,m1,&a);
! 223: subp(CO,(P)ONE,a,&s);
! 224: divcp(s,q,&m); cm2p(q,q2,m,&d);
! 225:
! 226: /* q1 = (d*g0) / f1; r1 = (d*g0) mod f1; */
! 227: mulp(CO,d,g1,&m); udivpwm(q,m,f1,&q1,&r1);
! 228:
! 229: /* g1 = g1 + (r1 mod q )*q; */
! 230: cm2p(q,q2,r1,&rm);
! 231: mulpq(rm,(P)q,&rm1); addp(CO,g1,rm1,&a);
! 232: g1 = a;
! 233:
! 234: /* g0 = g0 + (c*g0 + q1*f0 mod q)*q; */
! 235: mulp(CO,d,g0,&m); mulp(CO,q1,f0,&m1); addp(CO,m,m1,&a);
! 236: cm2p(q,q2,a,&a1); mulpq(a1,(P)q,&a2);
! 237: addp(CO,g0,a2,&a);
! 238: g0 = a;
! 239:
! 240: /* q = q^2; */
! 241: mulq(q,q,&qq);
! 242: q = qq;
! 243: divin(NM(q),2,&qn); NTOQ(qn,1,q2);
! 244: }
! 245: *fr0p = f0; *fr1p = f1; *gr0p = g0; *gr1p = g1; *qrp = q;
1.1 noro 246: }
247:
1.8 noro 248: void henzq1(P g,P h,Q bound,P *gcp,P *hcp,Q *qp)
1.1 noro 249: {
1.9 ! noro 250: V v;
! 251: Q f,q,q1;
! 252: Q u,t,a,b,s;
! 253: P c,c1;
! 254: P tg,th,tr;
! 255: UM wg,wh,ws,wt,wm;
! 256: int n,m,modres,mod,index,i;
! 257: P gc0,hc0;
! 258: P z,zz,zzz;
! 259:
! 260:
! 261: v = VR(g); n=deg(v,g); m=deg(v,h);
! 262: norm(g,&a); norm(h,&b);
! 263: STOQ(m,u); pwrq(a,u,&t);
! 264: STOQ(n,u); pwrq(b,u,&s);
! 265: mulq(t,s,&u);
! 266:
! 267: factorial(n+m,&t); mulq(u,t,&s);
! 268: addq(s,s,&f);
! 269:
! 270: wg = W_UMALLOC(m+n); wh = W_UMALLOC(m+n);
! 271: wt = W_UMALLOC(m+n); ws = W_UMALLOC(m+n);
! 272: wm = W_UMALLOC(m+n);
! 273:
! 274: for ( q = ONE, t = 0, c = 0, index = 0; ; ) {
! 275: mod = get_lprime(index++);
! 276: if ( !rem(NM((Q)LC(g)),mod) || !rem(NM((Q)LC(h)),mod) )
! 277: continue;
! 278: ptomp(mod,g,&tg); ptomp(mod,h,&th);
! 279: srchump(mod,tg,th,&tr);
! 280: if ( !tr )
! 281: continue;
! 282: else
! 283: modres = CONT((MQ)tr);
! 284:
! 285: mptoum(tg,wg); mptoum(th,wh);
! 286: eucum(mod,wg,wh,ws,wt); /* wg * ws + wh * wt = 1 */
! 287: DEG(wm) = 0; COEF(wm)[0] = modres;
! 288: mulum(mod,ws,wm,wt);
! 289: for ( i = DEG(wt); i >= 0; i-- )
! 290: if ( ( COEF(wt)[i] * 2 ) > mod )
! 291: COEF(wt)[i] -= mod;
! 292: chnrem(mod,v,c,q,wt,&c1,&q1);
! 293: if ( !ucmpp(c,c1) ) {
! 294: mulp(CO,c,g,&z); divsrp(CO,z,h,&zz,&zzz);
! 295: if ( NUM(zzz) ) {
! 296: q = q1; c = c1;
! 297: break;
! 298: }
! 299: }
! 300: q = q1; c = c1;
! 301:
! 302: if ( cmpq(f,q) < 0 )
! 303: break;
! 304: }
! 305: ptozp(c,1,&s,&gc0);
! 306: /* gc0*g = zz*h + zzz -> c*g + (-zz)*h = zzz */
! 307: mulp(CO,gc0,g,&z);
! 308: divsrp(CO,z,h,&zz,&zzz);
! 309: ptozp(zz,1,&s,(P *)&t);
! 310: if ( INT((Q)s) )
! 311: chsgnp(zz,&hc0);
! 312: else {
! 313: NTOQ(DN((Q)s),1,q); mulq(q,(Q)zzz,&q1); zzz = (P)q1;
! 314: mulp(CO,gc0,(P)q,&c); gc0 = c; mulp(CO,zz,(P)q,&c); chsgnp(c,&hc0);
! 315: }
! 316: if ( !INT((Q)zzz) ) {
! 317: NTOQ(DN((Q)zzz),1,q); NTOQ(NM((Q)zzz),1,q1); zzz = (P)q1;
! 318: mulp(CO,gc0,(P)q,&c); gc0 = c; mulp(CO,hc0,(P)q,&c); hc0 = c;
! 319: }
! 320: for ( index = 0; ; ) {
! 321: mod = get_lprime(index++);
! 322: if ( !rem(NM((Q)zzz),mod) ||
! 323: !rem(NM((Q)LC(g)),mod) ||
! 324: !rem(NM((Q)LC(h)),mod) )
! 325: continue;
! 326: for ( STOQ(mod,q); cmpq(q,bound) < 0; ) {
! 327: mulq(q,q,&q1); q = q1;
! 328: }
! 329: *qp = q;
! 330: invl((Q)zzz,q,&q1);
! 331: mulp(CO,gc0,(P)q1,hcp); mulp(CO,hc0,(P)q1,gcp);
! 332: return;
! 333: }
1.1 noro 334: }
335:
1.8 noro 336: void addm2p(VL vl,Q mod,Q mod2,P n1,P n2,P *nr)
1.1 noro 337: {
1.9 ! noro 338: P t;
1.1 noro 339:
1.9 ! noro 340: addp(vl,n1,n2,&t);
! 341: if ( !t )
! 342: *nr = 0;
! 343: else
! 344: cm2p(mod,mod2,t,nr);
1.1 noro 345: }
346:
1.8 noro 347: void subm2p(VL vl,Q mod,Q mod2,P n1,P n2,P *nr)
1.1 noro 348: {
1.9 ! noro 349: P t;
1.1 noro 350:
1.9 ! noro 351: subp(vl,n1,n2,&t);
! 352: if ( !t )
! 353: *nr = 0;
! 354: else
! 355: cm2p(mod,mod2,t,nr);
1.1 noro 356: }
357:
1.8 noro 358: void mulm2p(VL vl,Q mod,Q mod2,P n1,P n2,P *nr)
1.1 noro 359: {
1.9 ! noro 360: P t;
1.1 noro 361:
1.9 ! noro 362: mulp(vl,n1,n2,&t);
! 363: if ( !t )
! 364: *nr = 0;
! 365: else
! 366: cm2p(mod,mod2,t,nr);
1.1 noro 367: }
368:
1.8 noro 369: void cmp(Q mod,P p,P *pr)
1.1 noro 370: {
1.9 ! noro 371: P t;
! 372: DCP dc,dcr,dcr0;
1.1 noro 373:
1.9 ! noro 374: if ( !p )
! 375: *pr = 0;
! 376: else if ( NUM(p) )
! 377: remq((Q)p,mod,(Q *)pr);
! 378: else {
! 379: for ( dc = DC(p), dcr0 = 0; dc; dc = NEXT(dc) ) {
! 380: cmp(mod,COEF(dc),&t);
! 381: if ( t ) {
! 382: NEXTDC(dcr0,dcr);
! 383: DEG(dcr) = DEG(dc);
! 384: COEF(dcr) = t;
! 385: }
! 386: }
! 387: if ( !dcr0 )
! 388: *pr = 0;
! 389: else {
! 390: NEXT(dcr) = 0;
! 391: MKP(VR(p),dcr0,*pr);
! 392: }
! 393: }
1.1 noro 394: }
395:
1.8 noro 396: void cm2p(Q mod,Q m,P p,P *pr)
1.1 noro 397: {
1.9 ! noro 398: P t;
! 399: DCP dc,dcr,dcr0;
1.1 noro 400:
1.9 ! noro 401: if ( !p )
! 402: *pr = 0;
! 403: else if ( NUM(p) )
! 404: rem2q((Q)p,mod,m,(Q *)pr);
! 405: else {
! 406: for ( dc = DC(p), dcr0 = 0; dc; dc = NEXT(dc) ) {
! 407: cm2p(mod,m,COEF(dc),&t);
! 408: if ( t ) {
! 409: NEXTDC(dcr0,dcr);
! 410: DEG(dcr) = DEG(dc);
! 411: COEF(dcr) = t;
! 412: }
! 413: }
! 414: if ( !dcr0 )
! 415: *pr = 0;
! 416: else {
! 417: NEXT(dcr) = 0;
! 418: MKP(VR(p),dcr0,*pr);
! 419: }
! 420: }
1.1 noro 421: }
422:
1.8 noro 423: void addm2q(Q mod,Q mod2,Q n1,Q n2,Q *nr)
1.1 noro 424: {
1.9 ! noro 425: Q t;
1.1 noro 426:
1.9 ! noro 427: addq(n1,n2,&t);
! 428: if ( !t )
! 429: *nr = 0;
! 430: else
! 431: rem2q(t,mod,mod2,nr);
1.1 noro 432: }
433:
1.8 noro 434: void subm2q(Q mod,Q mod2,Q n1,Q n2,Q *nr)
1.1 noro 435: {
1.9 ! noro 436: Q t;
1.1 noro 437:
1.9 ! noro 438: subq(n1,n2,&t);
! 439: if ( !t )
! 440: *nr = 0;
! 441: else
! 442: rem2q(t,mod,mod2,nr);
1.1 noro 443: }
444:
1.8 noro 445: void mulm2q(Q mod,Q mod2,Q n1,Q n2,Q *nr)
1.1 noro 446: {
1.9 ! noro 447: Q t;
1.1 noro 448:
1.9 ! noro 449: mulq(n1,n2,&t);
! 450: if ( !t )
! 451: *nr = 0;
! 452: else
! 453: rem2q(t,mod,mod2,nr);
1.1 noro 454: }
455:
1.8 noro 456: void rem2q(Q n,Q m,Q m2,Q *nr)
1.1 noro 457: {
1.9 ! noro 458: N q,r,s;
! 459: int sgn;
1.1 noro 460:
1.9 ! noro 461: divn(NM(n),NM(m),&q,&r);
! 462: if ( !r )
! 463: *nr = 0;
! 464: else {
! 465: sgn = cmpn(r,NM(m2));
! 466: if ( sgn > 0 ) {
! 467: subn(NM(m),r,&s); NTOQ(s,-1*SGN(n),*nr);
! 468: } else
! 469: NTOQ(r,SGN(n),*nr);
! 470: }
1.6 noro 471: }
472:
473: /*
1.9 ! noro 474: extract d-homogeneous part with respect to vl - {v}
1.6 noro 475: */
476:
1.8 noro 477: void exthpc_generic(VL vl,P p,int d,V v,P *pr)
1.6 noro 478: {
1.9 ! noro 479: P w,x,t,t1,a,xd;
! 480: V v0;
! 481: DCP dc;
! 482:
! 483: if ( d < 0 || !p )
! 484: *pr = 0;
! 485: else if ( NUM(p) )
! 486: if ( d == 0 )
! 487: *pr = p;
! 488: else
! 489: *pr = 0;
! 490: else if ( v == VR(p) )
! 491: exthpc(vl,v,p,d,pr);
! 492: else {
! 493: v0 = VR(p);
! 494: for ( MKV(v0,x), dc = DC(p), w = 0; dc; dc = NEXT(dc) ) {
! 495: exthpc_generic(vl,COEF(dc),d-QTOS(DEG(dc)),v,&t);
! 496: pwrp(vl,x,DEG(dc),&xd);
! 497: mulp(vl,xd,t,&t1); addp(vl,w,t1,&a); w = a;
! 498: }
! 499: *pr = w;
! 500: }
1.1 noro 501: }
502:
1.8 noro 503: void exthp(VL vl,P p,int d,P *pr)
1.1 noro 504: {
1.9 ! noro 505: P t,t1,a,w,x,xd;
! 506: DCP dc;
1.1 noro 507:
1.9 ! noro 508: if ( d < 0 )
! 509: *pr = 0;
! 510: else if ( NUM(p) )
! 511: if ( d == 0 )
! 512: *pr = p;
! 513: else
! 514: *pr = 0;
! 515: else {
! 516: for ( MKV(VR(p),x), dc = DC(p), w = 0; dc; dc = NEXT(dc) ) {
! 517: exthp(vl,COEF(dc),d - QTOS(DEG(dc)),&t);
! 518: pwrp(vl,x,DEG(dc),&xd);
! 519: mulp(vl,xd,t,&t1); addp(vl,w,t1,&a); w = a;
! 520: }
! 521: *pr = w;
! 522: }
1.1 noro 523: }
524:
1.8 noro 525: void exthpc(VL vl,V v,P p,int d,P *pr)
1.1 noro 526: {
1.9 ! noro 527: P t,t1,a,w,x,xd;
! 528: DCP dc;
1.1 noro 529:
1.9 ! noro 530: if ( v != VR(p) )
! 531: exthp(vl,p,d,pr);
! 532: else if ( d < 0 )
! 533: *pr = 0;
! 534: else {
! 535: for ( MKV(v,x), dc = DC(p), w = 0; dc; dc = NEXT(dc) ) {
! 536: exthp(vl,COEF(dc),d,&t);
! 537: pwrp(vl,x,DEG(dc),&xd);
! 538: mulp(vl,xd,t,&t1); addp(vl,w,t1,&a); w = a;
! 539: }
! 540: *pr = w;
! 541: }
1.1 noro 542: }
543:
1.8 noro 544: void cbound(VL vl,P p,Q *b)
1.1 noro 545: {
1.9 ! noro 546: Q n,e,t,m;
! 547: int k;
1.1 noro 548:
1.9 ! noro 549: cmax(p,&n);
! 550: addq(n,n,&m);
1.1 noro 551:
1.9 ! noro 552: k = geldb(vl,p);
! 553: STOQ(3,t); STOQ(k,e);
1.1 noro 554:
1.9 ! noro 555: pwrq(t,e,&n);
! 556: mulq(m,n,b);
1.1 noro 557: }
558:
1.8 noro 559: int geldb(VL vl,P p)
1.1 noro 560: {
1.9 ! noro 561: int m;
1.1 noro 562:
1.9 ! noro 563: for ( m = 0; vl; vl = NEXT(vl) )
! 564: m += getdeg(vl->v,p);
! 565: return ( m );
1.1 noro 566: }
567:
1.8 noro 568: int getdeg(V v,P p)
1.1 noro 569: {
1.9 ! noro 570: int m,t;
! 571: DCP dc;
! 572:
! 573: if ( !p || NUM(p) )
! 574: return ( 0 );
! 575: else if ( v == VR(p) )
! 576: return ( deg(v,p) );
! 577: else {
! 578: for ( dc = DC(p), m = 0; dc; dc = NEXT(dc) ) {
! 579: t = getdeg(v,COEF(dc));
! 580: m = MAX(m,t);
! 581: }
! 582: return ( m );
! 583: }
1.1 noro 584: }
585:
1.8 noro 586: void cmax(P p,Q *b)
1.1 noro 587: {
1.9 ! noro 588: DCP dc;
! 589: Q m,m1;
! 590: N tn;
! 591:
! 592: if ( NUM(p) ) {
! 593: tn = NM((Q)p);
! 594: NTOQ(tn,1,*b);
! 595: } else {
! 596: for ( dc = DC(p), m = 0; dc; dc = NEXT(dc) ) {
! 597: cmax(COEF(dc),&m1);
! 598: if ( cmpq(m1,m) > 0 )
! 599: m = m1;
! 600: }
! 601: *b = m;
! 602: }
1.1 noro 603: }
604:
1.8 noro 605: int nextbin(VN vn,int n)
1.1 noro 606: {
1.9 ! noro 607: int tmp,i,carry;
1.1 noro 608:
1.9 ! noro 609: if ( n == 0 )
! 610: return ( 1 );
1.1 noro 611:
1.9 ! noro 612: for ( i = n - 1, carry = 1; i >= 0; i-- ) {
! 613: tmp = vn[i].n + carry;
! 614: vn[i].n = tmp % 2;
! 615: carry = tmp / 2;
! 616: }
! 617: return ( carry );
1.1 noro 618: }
619:
1.8 noro 620: void mulsgn(VN vn,VN vnt,int n,VN vn1)
1.1 noro 621: {
1.9 ! noro 622: int i;
1.1 noro 623:
1.9 ! noro 624: for ( i = 0; vn[i].v; i++ )
! 625: vn1[i].n = vn[i].n;
! 626: for ( i = 0; i < n; i++ )
! 627: if ( vnt[i].n )
! 628: vn1[(int)vnt[i].v].n *= -1;
1.1 noro 629: }
630:
1.8 noro 631: void next(VN vn)
1.1 noro 632: {
1.9 ! noro 633: int i,m,n,tmp,carry;
1.1 noro 634:
1.9 ! noro 635: for ( m = 0, i = 0; vn[i].v; i++ )
! 636: m = MAX(m,ABS(vn[i].n));
! 637: if ( m == 0 ) {
! 638: vn[--i].n = 1;
! 639: return;
! 640: }
! 641: for ( carry = 1, m++, i--; ( i >= 0 ) && carry; i-- ) {
! 642: tmp = vn[i].n + carry;
! 643: vn[i].n = tmp % m;
! 644: carry = tmp / m;
! 645: }
! 646: if ( ( i == -1 ) && carry ) {
! 647: for ( i = 0; vn[i].v; i++ )
! 648: vn[i].n = 0;
! 649: vn[--i].n = m;
! 650: } else {
! 651: for ( n = 0, i = 0; vn[i].v; i++ )
! 652: n = MAX(n,ABS(vn[i].n));
! 653: if ( n < m - 1 )
! 654: vn[--i].n = m - 1;
! 655: }
1.1 noro 656: }
1.9 ! noro 657:
1.8 noro 658: void clctv(VL vl,P p,VL *nvlp)
1.1 noro 659: {
1.9 ! noro 660: int i,n;
! 661: VL tvl;
! 662: VN tvn;
! 663:
! 664: if ( !p || NUM(p) ) {
! 665: *nvlp = 0;
! 666: return;
! 667: }
! 668:
! 669: for ( n = 0, tvl = vl; tvl; tvl = NEXT(tvl), n++ );
! 670: tvn = (VN) ALLOCA((n+1)*sizeof(struct oVN));
! 671: for ( i = 0, tvl = vl; tvl; tvl = NEXT(tvl), i++ ) {
! 672: tvn[i].v = tvl->v;
! 673: tvn[i].n = 0;
! 674: }
1.1 noro 675:
1.9 ! noro 676: markv(tvn,n,p);
! 677: vntovl(tvn,n,nvlp);
1.1 noro 678: }
679:
1.8 noro 680: void markv(VN vn,int n,P p)
1.1 noro 681: {
1.9 ! noro 682: V v;
! 683: DCP dc;
! 684: int i;
! 685:
! 686: if ( NUM(p) )
! 687: return;
! 688: v = VR(p);
! 689: for ( i = 0, v = VR(p); i < n; i++ )
! 690: if ( v == vn[i].v ) {
! 691: vn[i].n = 1;
! 692: break;
! 693: }
! 694: for ( dc = DC(p); dc; dc = NEXT(dc) )
! 695: markv(vn,n,COEF(dc));
1.1 noro 696: }
1.9 ! noro 697:
1.8 noro 698: void vntovl(VN vn,int n,VL *vlp)
1.1 noro 699: {
1.9 ! noro 700: int i;
! 701: VL tvl,tvl0;
1.1 noro 702:
1.9 ! noro 703: for ( i = 0, tvl0 = 0; ; ) {
! 704: while ( ( i < n ) && ( vn[i].n == 0 ) ) i++;
! 705: if ( i == n )
! 706: break;
! 707: else {
! 708: if ( !tvl0 ) {
! 709: NEWVL(tvl0);
! 710: tvl = tvl0;
! 711: } else {
! 712: NEWVL(NEXT(tvl));
! 713: tvl = NEXT(tvl);
! 714: }
! 715: tvl->v = vn[i++].v;
! 716: }
! 717: }
! 718: if ( tvl0 )
! 719: NEXT(tvl) = 0;
! 720: *vlp = tvl0;
1.1 noro 721: }
722:
1.8 noro 723: int dbound(V v,P f)
1.1 noro 724: {
1.9 ! noro 725: int m,t;
! 726: DCP dc;
1.1 noro 727:
1.9 ! noro 728: if ( !f )
! 729: return ( -1 );
! 730: else if ( v != VR(f) )
! 731: return homdeg(f);
! 732: else {
! 733: for ( dc = DC(f), m = 0; dc; dc = NEXT(dc) ) {
! 734: t = homdeg(COEF(dc));
! 735: m = MAX(m,t);
! 736: }
! 737: return ( m );
! 738: }
1.1 noro 739: }
740:
1.8 noro 741: int homdeg(P f)
1.1 noro 742: {
1.9 ! noro 743: int m,t;
! 744: DCP dc;
1.1 noro 745:
1.9 ! noro 746: if ( !f )
! 747: return ( -1 );
! 748: else if ( NUM(f) )
! 749: return( 0 );
! 750: else {
! 751: for ( dc = DC(f), m = 0; dc; dc = NEXT(dc) ) {
! 752: t = QTOS(DEG(dc))+homdeg(COEF(dc));
! 753: m = MAX(m,t);
! 754: }
! 755: return ( m );
! 756: }
1.1 noro 757: }
758:
1.8 noro 759: int minhomdeg(P f)
1.1 noro 760: {
1.9 ! noro 761: int m,t;
! 762: DCP dc;
1.1 noro 763:
1.9 ! noro 764: if ( !f )
! 765: return ( -1 );
! 766: else if ( NUM(f) )
! 767: return( 0 );
! 768: else {
! 769: for ( dc = DC(f), m = homdeg(f); dc; dc = NEXT(dc) ) {
! 770: t = QTOS(DEG(dc))+minhomdeg(COEF(dc));
! 771: m = MIN(m,t);
! 772: }
! 773: return ( m );
! 774: }
1.1 noro 775: }
776:
1.8 noro 777: void adjc(VL vl,P f,P a,P lc0,Q q,P *fr,P *ar)
1.1 noro 778: {
1.9 ! noro 779: P m,m1;
! 780: Q t;
1.1 noro 781:
1.9 ! noro 782: invl((Q)LC(f),q,&t);
! 783: mulq((Q)lc0,t,(Q *)&m);
! 784: mulpq(f,m,&m1); cmp(q,m1,fr);
! 785: invl((Q)m,q,&t);
! 786: mulpq(a,(P)t,&m1);
! 787: cmp(q,m1,ar);
1.1 noro 788: }
789:
790: #if 1
1.8 noro 791: void affinemain(VL vl,P p,V v0,int n,P *pl,P *pr)
1.1 noro 792: {
1.9 ! noro 793: P x,t,m,c,s,a;
! 794: DCP dc;
! 795: Q d;
! 796:
! 797: if ( !p )
! 798: *pr = 0;
! 799: else if ( NUM(p) )
! 800: *pr = p;
! 801: else if ( VR(p) != v0 ) {
! 802: MKV(VR(p),x);
! 803: for ( c = 0, dc = DC(p); dc; dc = NEXT(dc) ) {
! 804: affinemain(vl,COEF(dc),v0,n,pl,&t);
! 805: if ( DEG(dc) ) {
! 806: pwrp(vl,x,DEG(dc),&s); mulp(vl,s,t,&m);
! 807: addp(vl,m,c,&a); c = a;
! 808: } else {
! 809: addp(vl,t,c,&a); c = a;
! 810: }
! 811: }
! 812: *pr = c;
! 813: } else {
! 814: dc = DC(p);
! 815: c = COEF(dc);
! 816: for ( d = DEG(dc), dc = NEXT(dc);
! 817: dc; d = DEG(dc), dc = NEXT(dc) ) {
! 818: mulp(vl,pl[QTOS(d)-QTOS(DEG(dc))],c,&m);
! 819: addp(vl,m,COEF(dc),&c);
! 820: }
! 821: if ( d ) {
! 822: mulp(vl,pl[QTOS(d)],c,&m); c = m;
! 823: }
! 824: *pr = c;
! 825: }
1.1 noro 826: }
827: #endif
828:
829: #if 0
1.8 noro 830: affine(VL vl,P p,VN vn,P *r)
1.1 noro 831: {
1.9 ! noro 832: int n,d,d1,i;
! 833: Q *q;
! 834: Q **bc;
! 835:
! 836: if ( !p || NUM(p) )
! 837: *r = p;
! 838: else {
! 839: for ( i = 0, d = 0; vn[i].v; i++ )
! 840: d1 = getdeg(vn[i].v,p), d = MAX(d,d1);
! 841: W_CALLOC(d+1,Q *,bc);
! 842: for ( i = 0; i <= d; i++ )
! 843: W_CALLOC(i+1,Q,q), mkbc(i,q), bc[i] = q;
! 844: afmain(vl,bc,p,vn,r);
! 845: }
1.1 noro 846: }
847:
1.8 noro 848: afmain(VL vl,Q **bc,P p,VN vn,P *r)
1.1 noro 849: {
1.9 ! noro 850: P t,s,u;
! 851: P *c,*rc;
! 852: Q *q;
! 853: DCP dc;
! 854: int n,i,j;
! 855:
! 856: if ( !p || NUM(p) || !vn[0].v )
! 857: *r = p;
! 858: else if ( vn[0].v != VR(p) ) {
! 859: for ( i = 0; vn[i].v && (vn[i].v != VR(p)); i++ );
! 860: if ( vn[i].v )
! 861: afmain(vl,bc,p,vn+i,r);
! 862: else {
! 863: n = QTOS(DEG(DC(p))); W_CALLOC(n+1,P,c);
! 864: for ( dc = DC(p); dc; dc = NEXT(dc) )
! 865: afmain(vl,bc,COEF(dc),vn,&c[QTOS(DEG(dc))]);
! 866: plisttop(c,VR(p),n,r);
! 867: }
! 868: } else {
! 869: n = QTOS(DEG(DC(p))); W_CALLOC(n+1,P,c);
! 870: W_CALLOC(n+1,P,rc); W_CALLOC(n+1,Q,q);
! 871: for ( dc = DC(p); dc; dc = NEXT(dc) )
! 872: afmain(vl,bc,COEF(dc),vn+1,&c[QTOS(DEG(dc))]);
! 873: if ( !vn[0].n )
! 874: bcopy(c,rc,sizeof(P)*(n+1));
! 875: else {
! 876: for ( q[0] = ONE, STOQ(vn[0].n,q[1]), i = 2; i <= n; i++ )
! 877: mulq(q[i-1],q[1],&q[i]);
! 878: for ( j = 0; j <= n; rc[j] = t, j++ )
! 879: for ( i = j, t = 0; i <= n; i++ )
! 880: if ( c[i] )
! 881: mulq(bc[i][j],q[i-j],&s), mulp(vl,s,c[i],&u),
! 882: addp(CO,u,t,&s), t = s;
! 883: }
! 884: plisttop(rc,VR(p),n,r);
! 885: }
1.1 noro 886: }
887: #endif
888:
1.8 noro 889: void restore(VL vl,P f,VN vn,P *fr)
1.1 noro 890: {
1.9 ! noro 891: int i;
! 892: P vv,g,g1,t;
! 893: Q s;
! 894:
! 895: g = f;
! 896: for ( i = 0; vn[i].v; i++ ) {
! 897: MKV(vn[i].v,t);
! 898: if ( vn[i].n ) {
! 899: STOQ(-vn[i].n,s);
! 900: addp(vl,t,(P)s,&vv);
! 901: } else
! 902: vv = t;
! 903:
! 904: substp(vl,g,vn[i].v,vv,&g1); g = g1;
! 905: }
! 906: *fr = g;
1.1 noro 907: }
908:
1.8 noro 909: void mergev(VL vl,VL vl1,VL vl2,VL *nvlp)
1.1 noro 910: {
1.9 ! noro 911: int i,n;
! 912: VL tvl;
! 913: VN vn;
! 914:
! 915: if ( !vl1 ) {
! 916: *nvlp = vl2; return;
! 917: } else if ( !vl2 ) {
! 918: *nvlp = vl1; return;
! 919: }
! 920: for ( i = 0, tvl = vl; tvl; i++, tvl = NEXT(tvl) );
! 921: n = i;
! 922: W_CALLOC(n,struct oVN,vn);
! 923: for ( i = 0, tvl = vl; tvl; tvl = NEXT(tvl), i++ )
! 924: vn[i].v = tvl->v;
! 925: for ( i = 0, tvl = vl1; (i < n) && tvl; i++, tvl = NEXT(tvl)) {
! 926: while ( ( i < n ) && ( vn[i].v != tvl->v ) )
! 927: i++;
! 928: if ( i == n )
! 929: break;
! 930: else
! 931: vn[i].n = 1;
! 932: }
! 933: for ( i = 0, tvl = vl2; (i < n) && tvl; i++, tvl = NEXT(tvl) ) {
! 934: while ( ( i < n ) && ( vn[i].v != tvl->v ) )
! 935: i++;
! 936: if ( i == n )
! 937: break;
! 938: else
! 939: vn[i].n = 1;
! 940: }
! 941: vntovl(vn,n,nvlp);
1.1 noro 942: }
943:
944: #if 0
1.8 noro 945: void substvp(VL vl,P f,VN vn,P *g)
1.1 noro 946: {
1.9 ! noro 947: V v;
! 948: int i;
! 949: P h,h1;
! 950: Q t;
! 951:
! 952: h = f;
! 953: for ( i = 0; v = vn[i].v; i++ ) {
! 954: STOQ(vn[i].n,t);
! 955: substp(vl,h,v,(P)t,&h1); h = h1;
! 956: }
! 957: *g = h;
1.1 noro 958: }
959:
1.8 noro 960: void affine(VL vl,P f,VN vn,P *fr)
1.1 noro 961: {
1.9 ! noro 962: int i,j,n;
! 963: P vv,g,g1,t,u;
! 964: Q s;
! 965: int *dlist;
! 966: P **plist;
! 967:
! 968: for ( n = 0; vn[n].v; n++);
! 969: dlist = (int *)ALLOCA((n+1)*sizeof(int));
! 970: plist = (P **)ALLOCA((n+1)*sizeof(P *));
! 971: for ( i = 0; vn[i].v; i++ ) {
! 972: if ( !vn[i].n )
! 973: continue;
! 974: dlist[i] = getdeg(vn[i].v,f);
! 975: plist[i] = (P *)ALLOCA((dlist[i]+1)*sizeof(P));
! 976:
! 977: MKV(vn[i].v,t);
! 978: if ( vn[i].n ) {
! 979: STOQ(vn[i].n,s); addp(vl,t,(P)s,&vv);
! 980: } else
! 981: vv = t;
! 982:
! 983: for ( j = 0, t = (P)ONE; j < dlist[i]; j++ ) {
! 984: plist[i][j] = t;
! 985: mulp(vl,t,vv,&u);
! 986: t = u;
! 987: }
! 988: plist[i][j] = t;
! 989: }
! 990:
! 991: g = f;
! 992: for ( i = 0; vn[i].v; i++ ) {
! 993: if ( !vn[i].n )
! 994: continue;
! 995: affinemain(vl,g,vn[i].v,dlist[i],plist[i],&g1); g = g1;
! 996: }
! 997: *fr = g;
1.1 noro 998: }
999: #endif
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 / asir2000 / engine