Annotation of OpenXM_contrib2/asir2018/builtin/isolv.c, Revision 1.1
1.1 ! noro 1: /*
! 2: * $OpenXM$
! 3: */
! 4:
! 5: #include "ca.h"
! 6: #include "parse.h"
! 7: #include "version.h"
! 8:
! 9: #if defined(INTERVAL)
! 10:
! 11: static void Solve(NODE, Obj *);
! 12: static void NSolve(NODE, Obj *);
! 13:
! 14: void Solve1(P, Q, pointer *);
! 15: void Sturm(P, VECT *);
! 16: void boundbody(P, Q *);
! 17: void binary(int , MAT);
! 18: void separate(Q, Q, VECT, Q, Q, int, int, MAT, int *);
! 19: void ueval(P, Q, Q *);
! 20:
! 21: struct ftab isolv_tab[] = {
! 22: {"solve", Solve, 2},
! 23: {"nsolve", NSolve, 2},
! 24: {0,0,0},
! 25: };
! 26:
! 27: static void
! 28: Solve(arg, rp)
! 29: NODE arg;
! 30: Obj *rp;
! 31: {
! 32: pointer p, Eps;
! 33: pointer root;
! 34: V v;
! 35: Q eps;
! 36:
! 37: p = (pointer)ARG0(arg);
! 38: if ( !p ) {
! 39: *rp = 0;
! 40: return;
! 41: }
! 42: Eps = (pointer)ARG1(arg);
! 43: asir_assert(Eps, O_N, "solve");
! 44: if ( NID(Eps) != N_Q ) {
! 45: fprintf(stderr,"solve arg 2 is required a rational number");
! 46: error(" : invalid argument");
! 47: return;
! 48: }
! 49: DUPQ((Q)Eps, eps);
! 50: SGN(eps) = 1;
! 51: switch (OID(p)) {
! 52: case O_N:
! 53: *rp = 0;
! 54: break;
! 55: case O_P:
! 56: Pvars(arg, &root);
! 57: if (NEXT(BDY((LIST)root)) != 0) {
! 58: fprintf(stderr,"solve arg 1 is univariate polynormial");
! 59: error(" : invalid argument");
! 60: break;
! 61: }
! 62: Solve1((P)p, eps, &root);
! 63: *rp = (Obj)root;
! 64: break;
! 65: case O_LIST:
! 66: fprintf(stderr,"solve,");
! 67: error(" : Sorry, not yet implement of multivars");
! 68: break;
! 69: default:
! 70: *rp = 0;
! 71: }
! 72: }
! 73:
! 74: static void
! 75: NSolve(arg, rp)
! 76: NODE arg;
! 77: Obj *rp;
! 78: {
! 79: pointer p, Eps;
! 80: pointer root;
! 81: LIST listp;
! 82: V v;
! 83: Q eps;
! 84: NODE n, n0, m0, m, ln0;
! 85: Num r;
! 86: Itv iv;
! 87: BF breal;
! 88:
! 89: p = (pointer)ARG0(arg);
! 90: if ( !p ) {
! 91: *rp = 0;
! 92: return;
! 93: }
! 94: Eps = (pointer)ARG1(arg);
! 95: asir_assert(Eps, O_N, "solve");
! 96: if ( NID(Eps) != N_Q ) {
! 97: fprintf(stderr,"solve arg 2 is required a rational number");
! 98: error(" : invalid argument");
! 99: return;
! 100: }
! 101: DUPQ((Q)Eps, eps);
! 102: SGN(eps) = 1;
! 103: switch (OID(p)) {
! 104: case O_N:
! 105: *rp = 0;
! 106: break;
! 107: case O_P:
! 108: Pvars(arg, &root);
! 109: if (NEXT(BDY((LIST)root)) != 0) {
! 110: fprintf(stderr,"solve arg 1 is univariate polynormial");
! 111: error(" : invalid argument");
! 112: break;
! 113: }
! 114: Solve1((P)p, eps, &root);
! 115: for (m0 = BDY((LIST)root), n0 = 0; m0; m0 = NEXT(m0)) {
! 116: m = BDY((LIST)BDY(m0));
! 117: miditvp(BDY(m), &r);
! 118: ToBf(r, &breal);
! 119: NEXTNODE( n0, n );
! 120: MKNODE(ln0, breal, NEXT(m));
! 121: MKLIST(listp, ln0);
! 122: BDY(n) = (pointer)listp;
! 123: }
! 124: NEXT(n) = 0;
! 125: MKLIST(listp,n0);
! 126: *rp = (pointer)listp;
! 127: break;
! 128: case O_LIST:
! 129: fprintf(stderr,"solve,");
! 130: error(" : Sorry, not yet implement of multivars");
! 131: break;
! 132: default:
! 133: *rp = 0;
! 134: }
! 135: }
! 136:
! 137: void
! 138: Solve1(inp, eps, rt)
! 139: P inp;
! 140: Q eps;
! 141: pointer *rt;
! 142: {
! 143: P p;
! 144: Q up, low, a;
! 145: DCP fctp, onedeg, zerodeg;
! 146: LIST listp;
! 147: VECT sseq;
! 148: MAT root;
! 149: int chvu, chvl, pad, tnumb, numb, i, j;
! 150: Itv iv;
! 151: NODE n0, n, ln0, ln;
! 152:
! 153: boundbody(inp, &up);
! 154: if (!up) {
! 155: *rt = 0;
! 156: return;
! 157: }
! 158: Sturm(inp, &sseq);
! 159: DUPQ(up,low);
! 160: SGN(low) = -1;
! 161: chvu = stumq(sseq, up);
! 162: chvl = stumq(sseq, low);
! 163: tnumb = abs(chvu - chvl);
! 164: MKMAT(root, tnumb, 4);
! 165: pad = -1;
! 166:
! 167: fctrp(CO,inp,&fctp);
! 168: for (fctp = NEXT(fctp), i = 0; fctp; fctp = NEXT(fctp)) {
! 169: p = COEF(fctp);
! 170: onedeg = DC(p);
! 171: if ( !cmpq(DEG(onedeg), ONE) ) {
! 172: pad++;
! 173: if ( !NEXT(onedeg) ) {
! 174: BDY(root)[pad][0] = 0;
! 175: BDY(root)[pad][1] = 0;
! 176: BDY(root)[pad][2] = DEG(fctp);
! 177: BDY(root)[pad][3] = p;
! 178: } else {
! 179: divq((Q)COEF(NEXT(onedeg)),(Q)COEF(onedeg),&a);
! 180: BDY(root)[pad][0] = a;
! 181: BDY(root)[pad][1] = BDY(root)[pad][0];
! 182: BDY(root)[pad][2] = DEG(fctp);
! 183: BDY(root)[pad][3] = p;
! 184: }
! 185: continue;
! 186: }
! 187: boundbody(p, &up);
! 188: Sturm(p, &sseq);
! 189: DUPQ(up,low);
! 190: SGN(low) = -1;
! 191: chvu = stumq(sseq, up);
! 192: chvl = stumq(sseq, low);
! 193: numb = abs(chvu - chvl);
! 194: separate(DEG(fctp), eps, sseq, up, low, chvu, chvl, root, &pad);
! 195: }
! 196: for (i = 0; i < pad; i++) {
! 197: for (j = i; j <= pad; j++) {
! 198: if (cmpq(BDY(root)[i][0], BDY(root)[j][0]) > 0) {
! 199: a = BDY(root)[i][0];
! 200: BDY(root)[i][0] = BDY(root)[j][0];
! 201: BDY(root)[j][0] = a;
! 202: a = BDY(root)[i][1];
! 203: BDY(root)[i][1] = BDY(root)[j][1];
! 204: BDY(root)[j][1] = a;
! 205: a = BDY(root)[i][2];
! 206: BDY(root)[i][2] = BDY(root)[j][2];
! 207: BDY(root)[j][2] = a;
! 208: a = BDY(root)[i][3];
! 209: BDY(root)[i][3] = BDY(root)[j][3];
! 210: BDY(root)[j][3] = a;
! 211: }
! 212: }
! 213: }
! 214: for (i = 0; i < pad; i++) {
! 215: while(cmpq(BDY(root)[i][1], BDY(root)[i+1][0]) > 0 ) {
! 216: binary(i, root);
! 217: binary(i+1, root);
! 218: if ( cmpq(BDY(root)[i][0], BDY(root)[i+1][1]) > 0 ) {
! 219: a = BDY(root)[i][0];
! 220: BDY(root)[i][0] = BDY(root)[i+1][0];
! 221: BDY(root)[i+1][0] = a;
! 222: a = BDY(root)[i][1];
! 223: BDY(root)[i][1] = BDY(root)[i+1][1];
! 224: BDY(root)[i+1][1] = a;
! 225: a = BDY(root)[i][2];
! 226: BDY(root)[i][2] = BDY(root)[i+1][2];
! 227: BDY(root)[i+1][2] = a;
! 228: a = BDY(root)[i][3];
! 229: BDY(root)[i][3] = BDY(root)[i+1][3];
! 230: BDY(root)[i+1][3] = a;
! 231: break;
! 232: }
! 233: }
! 234: }
! 235: for (i = 0, n0 = 0; i <= pad; i++) {
! 236: istoitv(BDY(root)[i][0], BDY(root)[i][1], &iv);
! 237: NEXTNODE(n0,n);
! 238: MKNODE(ln, BDY(root)[i][2], 0); MKNODE(ln0, iv, ln);
! 239: MKLIST(listp, ln0);BDY(n) = (pointer)listp;
! 240: }
! 241: NEXT(n) = 0;
! 242: MKLIST(listp,n0);
! 243: *rt = (pointer)listp;
! 244: }
! 245:
! 246: void
! 247: separate(mult, eps, sseq, up, low, upn, lown, root, padp)
! 248: VECT sseq;
! 249: Q mult, eps, up, low;
! 250: int upn, lown;
! 251: MAT root;
! 252: int *padp;
! 253: {
! 254: int de, midn;
! 255: Q mid, e;
! 256: P p;
! 257:
! 258: de = abs(lown - upn);
! 259: if (de == 0) return;
! 260: if (de == 1) {
! 261: (*padp)++;
! 262: BDY(root)[*padp][0] = up;
! 263: BDY(root)[*padp][1] = low;
! 264: BDY(root)[*padp][3] = (P *)sseq->body[0];
! 265: subq( BDY(root)[*padp][1], BDY(root)[*padp][0], &e );
! 266: SGN(e) = 1;
! 267: while (cmpq(e, eps) > 0) {
! 268: binary(*padp, root);
! 269: subq( BDY(root)[*padp][1], BDY(root)[*padp][0], &e);
! 270: SGN(e) = 1;
! 271: }
! 272: BDY(root)[*padp][2] = mult;
! 273: return;
! 274: }
! 275: addq(up, low, &mid);
! 276: divq(mid, TWO, &mid);
! 277: midn = stumq(sseq, mid);
! 278: separate(mult, eps, sseq, low, mid, lown, midn, root, padp);
! 279: separate(mult, eps, sseq, mid, up, midn, upn, root, padp);
! 280: }
! 281:
! 282: void
! 283: binary(indx, root)
! 284: int indx;
! 285: MAT root;
! 286: {
! 287: Q a, b, c, d, e;
! 288: P p;
! 289: p = (P)BDY(root)[indx][3];
! 290: addq(BDY(root)[indx][0], BDY(root)[indx][1], &c);
! 291: divq(c, TWO, &d);
! 292: ueval(p, BDY(root)[indx][1], &a);
! 293: ueval(p, d, &b);
! 294: if (SGN(a) == SGN(b)){
! 295: BDY(root)[indx][1] = d;
! 296: } else {
! 297: BDY(root)[indx][0] = d;
! 298: }
! 299: }
! 300:
! 301: void
! 302: Sturm(p, ret)
! 303: P p;
! 304: VECT *ret;
! 305: {
! 306: P g1,g2,q,r,s, *t;
! 307: Q a,b,c,d,h,l,m,x;
! 308: V v;
! 309: VECT seq;
! 310: int i,j;
! 311:
! 312: v = VR(p);
! 313: t = (P *)ALLOCA((deg(v,p)+1)*sizeof(P));
! 314: g1 = t[0] = p; diffp(CO,p,v,(P *)&a); ptozp((P)a,1,&c,&g2); t[1] = g2;
! 315: for ( i = 1, h = ONE, x = ONE; ; ) {
! 316: if ( NUM(g2) ) break;
! 317: subq(DEG(DC(g1)),DEG(DC(g2)),&d);
! 318: l = (Q)LC(g2);
! 319: if ( SGN(l) < 0 ) {
! 320: chsgnq(l,&a); l = a;
! 321: }
! 322: addq(d,ONE,&a); pwrq(l,a,&b); mulp(CO,(P)b,g1,(P *)&a);
! 323: divsrp(CO,(P)a,g2,&q,&r);
! 324: if ( !r ) break;
! 325: chsgnp(r,&s); r = s; i++;
! 326: if ( NUM(r) ) {
! 327: t[i] = r; break;
! 328: }
! 329: pwrq(h,d,&m); g1 = g2;
! 330: mulq(m,x,&a); divsp(CO,r,(P)a,&g2); t[i] = g2;
! 331: x = (Q)LC(g1);
! 332: if ( SGN(x) < 0 ) {
! 333: chsgnq(x,&a); x = a;
! 334: }
! 335: pwrq(x,d,&a); mulq(a,h,&b); divq(b,m,&h);
! 336: }
! 337: MKVECT(seq,i+1);
! 338: for ( j = 0; j <= i; j++ ) seq->body[j] = (pointer)t[j];
! 339: *ret = seq;
! 340: }
! 341:
! 342: int
! 343: stumq(s, val)
! 344: VECT s;
! 345: Q val;
! 346: {
! 347: int len, i, j, c;
! 348: P *ss;
! 349: Q a, a0;
! 350:
! 351: len = s->len;
! 352: ss = (P *)s->body;
! 353: for ( j = 0; j < len; j++ ){
! 354: ueval(ss[j],val,&a0);
! 355: if (a0) break;
! 356: }
! 357: for ( i = j++, c =0; i < len; i++) {
! 358: ueval( ss[i], val, &a);
! 359: if ( a ) {
! 360: if( (SGN(a) > 0 && SGN(a0) < 0) || (SGN(a) < 0 && SGN(a0) > 0) ){
! 361: c++;
! 362: a0 = a;
! 363: }
! 364: }
! 365: }
! 366: return c;
! 367: }
! 368:
! 369: void
! 370: boundbody(p, q)
! 371: P p;
! 372: Q *q;
! 373: {
! 374: Q t, max, tmp;
! 375: DCP dc;
! 376:
! 377: if ( !p )
! 378: *q = 0;
! 379: else if ( p->id == O_N )
! 380: *q = 0;
! 381: else {
! 382: NEWQ(tmp);
! 383: SGN(tmp)=1;
! 384: for ( dc = DC(p), max=0; dc; dc = NEXT(dc) ) {
! 385: t = (Q)COEF(dc);
! 386: NM(tmp)=NM(t);
! 387: DN(tmp)=DN(t);
! 388: if ( cmpq(tmp, max) > 0 ) DUPQ(tmp, max);
! 389: }
! 390: addq(ONE, max, q);
! 391: }
! 392: }
! 393:
! 394: void
! 395: ueval(p, q, rp)
! 396: P p;
! 397: Q q;
! 398: Q *rp;
! 399: {
! 400: Q d, d1, a, b, t;
! 401: Q deg, da;
! 402: Q nm, dn;
! 403: DCP dc;
! 404:
! 405: if ( !p ) *rp = 0;
! 406: else if ( NUM(p) ) *rp = (Q)p;
! 407: else {
! 408: if ( q ) {
! 409: NTOQ( DN(q), 1, dn );
! 410: NTOQ( NM(q), SGN(q), nm );
! 411: } else {
! 412: dn = 0;
! 413: nm = 0;
! 414: }
! 415: if ( !dn ) {
! 416: dc = DC(p); t = (Q)COEF(dc);
! 417: for ( d = DEG(dc), dc = NEXT(dc); dc; d = DEG(dc), dc= NEXT(dc) ) {
! 418: subq(d, DEG(dc), &d1); pwrq(nm, d1, &a);
! 419: mulq(t,a,&b); addq(b,(Q)COEF(dc),&t);
! 420: }
! 421: if ( d ) {
! 422: pwrq(nm,d,&a); mulq(t,a,&b); t = b;
! 423: }
! 424: *rp = t;
! 425: } else {
! 426: dc = DC(p); t = (Q)COEF(dc);
! 427: for ( d=deg= DEG(dc), dc = NEXT(dc); dc; d = DEG(dc), dc= NEXT(dc) ) {
! 428: subq(d, DEG(dc), &d1); pwrq(nm, d1, &a);
! 429: mulq(t,a,&b);
! 430: subq(deg, DEG(dc), &d1); pwrq(dn, d1, &a);
! 431: mulq(a, (Q)COEF(dc), &da);
! 432: addq(b,da,&t);
! 433: }
! 434: if ( d ) {
! 435: pwrq(nm,d,&a); mulq(t,a,&b); t = b;
! 436: }
! 437: *rp = t;
! 438: }
! 439: }
! 440: }
! 441: #endif
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to isolv.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib2 / asir2018 / builtin