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