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