![[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.5 ! kondoh 2: * $OpenXM: OpenXM_contrib2/asir2000/builtin/isolv.c,v 1.4 2005/02/08 18:06:05 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;
1.5 ! kondoh 69: default:
1.1 saito 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;
1.4 saito 80: pointer root;
81: LIST listp;
1.3 saito 82: V v;
83: Q eps;
84: NODE n, n0, m0, m, ln0;
85: Num r;
86: Itv iv;
87: BF breal;
1.1 saito 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));
1.3 saito 121: MKLIST((LIST)listp, ln0);
1.1 saito 122: BDY(n) = (pointer)listp;
123: }
124: NEXT(n) = 0;
1.3 saito 125: MKLIST((LIST)listp,n0);
1.1 saito 126: *rp = (pointer)listp;
127: break;
128: case O_LIST:
129: fprintf(stderr,"solve,");
130: error(" : Sorry, not yet implement of multivars");
131: break;
1.5 ! kondoh 132: default:
1.1 saito 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);
1.3 saito 164: MKMAT(root, tnumb, 4);
1.1 saito 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);
1.3 saito 177: BDY(root)[pad][3] = p;
1.1 saito 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);
1.3 saito 183: BDY(root)[pad][3] = p;
1.1 saito 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;
1.3 saito 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;
1.1 saito 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;
1.3 saito 255: Q mid, e;
1.1 saito 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;
1.3 saito 264: BDY(root)[*padp][3] = (P *)sseq->body[0];
1.1 saito 265: subq( BDY(root)[*padp][1], BDY(root)[*padp][0], &e );
266: SGN(e) = 1;
267: while (cmpq(e, eps) > 0) {
1.3 saito 268: binary(*padp, root);
1.1 saito 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);
1.3 saito 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;
1.4 saito 289: p = (P)BDY(root)[indx][3];
1.3 saito 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: }
1.1 saito 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