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