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