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