version 1.3, 2019/10/17 03:03:12 |
version 1.4, 2019/12/24 10:26:38 |
|
|
/* |
/* |
* $OpenXM: OpenXM_contrib2/asir2018/builtin/isolv.c,v 1.2 2019/06/04 07:11:23 kondoh Exp $ |
* $OpenXM: OpenXM_contrib2/asir2018/builtin/isolv.c,v 1.3 2019/10/17 03:03:12 kondoh Exp $ |
*/ |
*/ |
|
|
#include "ca.h" |
#include "ca.h" |
|
|
|
|
#if defined(INTERVAL) |
#if defined(INTERVAL) |
|
|
|
/* in Q.c */ |
|
void dupz(Z, Z*); |
|
|
static void Solve(NODE, Obj *); |
static void Solve(NODE, Obj *); |
static void NSolve(NODE, Obj *); |
static void NSolve(NODE, Obj *); |
|
|
|
|
NODE arg; |
NODE arg; |
Obj *rp; |
Obj *rp; |
{ |
{ |
pointer p, Eps; |
//pointer p, Eps; |
pointer root; |
pointer root; |
V v; |
V v; |
Q eps; |
Q eps; |
|
Q Eps; |
|
P p; |
|
|
|
|
*rp = 0; |
p = (P)ARG0(arg); |
#if 0 |
|
p = (pointer)ARG0(arg); |
|
if ( !p ) { |
if ( !p ) { |
*rp = 0; |
*rp = 0; |
return; |
return; |
} |
} |
Eps = (pointer)ARG1(arg); |
Eps = (Q)ARG1(arg); |
asir_assert(Eps, O_N, "solve"); |
asir_assert(Eps, O_N, "solve"); |
if ( NID(Eps) != N_Q ) { |
if ( NID(Eps) != N_Q ) { |
fprintf(stderr,"solve arg 2 is required a rational number"); |
fprintf(stderr,"solve arg 2 is required a rational number"); |
error(" : invalid argument"); |
error(" : invalid argument"); |
return; |
return; |
} |
} |
DUPQ((Q)Eps, eps); |
//DUPQ((Q)Eps, eps); |
SGN(eps) = 1; |
//SGN(eps) = 1; |
|
NEWQ(eps); mpq_abs(BDY(eps),BDY(Eps)); |
|
|
switch (OID(p)) { |
switch (OID(p)) { |
case O_N: |
case O_N: |
*rp = 0; |
*rp = 0; |
|
|
error(" : invalid argument"); |
error(" : invalid argument"); |
break; |
break; |
} |
} |
Solve1((P)p, eps, &root); |
Solve1(p, eps, &root); |
*rp = (Obj)root; |
*rp = (Obj)root; |
break; |
break; |
case O_LIST: |
case O_LIST: |
|
|
default: |
default: |
*rp = 0; |
*rp = 0; |
} |
} |
#endif |
|
} |
} |
|
|
static void |
static void |
NSolve(NODE arg, Obj *rp) |
NSolve(NODE arg, Obj *rp) |
{ |
{ |
pointer p, Eps; |
P p; |
|
Q Eps; |
|
Q eps; |
pointer root; |
pointer root; |
LIST listp; |
LIST listp; |
V v; |
V v; |
Q eps; |
|
NODE n, n0, m0, m, ln0; |
NODE n, n0, m0, m, ln0; |
Num r; |
Num r; |
Itv iv; |
Itv iv; |
BF breal; |
BF breal; |
|
|
*rp = 0; |
p = (P)ARG0(arg); |
#if 0 |
|
|
|
p = (pointer)ARG0(arg); |
|
if ( !p ) { |
if ( !p ) { |
*rp = 0; |
*rp = 0; |
return; |
return; |
} |
} |
Eps = (pointer)ARG1(arg); |
Eps = (Q)ARG1(arg); |
asir_assert(Eps, O_N, "solve"); |
asir_assert(Eps, O_N, "solve"); |
if ( NID(Eps) != N_Q ) { |
if ( NID(Eps) != N_Q ) { |
fprintf(stderr,"solve arg 2 is required a rational number"); |
fprintf(stderr,"solve arg 2 is required a rational number"); |
error(" : invalid argument"); |
error(" : invalid argument"); |
return; |
return; |
} |
} |
DUPQ((Q)Eps, eps); |
//DUPQ((Q)Eps, eps); |
SGN(eps) = 1; |
//SGN(eps) = 1; |
|
NEWQ(eps); mpq_abs(BDY(eps),BDY(Eps)); |
|
|
switch (OID(p)) { |
switch (OID(p)) { |
case O_N: |
case O_N: |
*rp = 0; |
*rp = 0; |
Line 125 NSolve(NODE arg, Obj *rp) |
|
Line 129 NSolve(NODE arg, Obj *rp) |
|
error(" : invalid argument"); |
error(" : invalid argument"); |
break; |
break; |
} |
} |
Solve1((P)p, eps, &root); |
Solve1(p, eps, &root); |
for (m0 = BDY((LIST)root), n0 = 0; m0; m0 = NEXT(m0)) { |
for (m0 = BDY((LIST)root), n0 = 0; m0; m0 = NEXT(m0)) { |
m = BDY((LIST)BDY(m0)); |
m = BDY((LIST)BDY(m0)); |
miditvp(BDY(m), &r); |
miditvp(BDY(m), &r); |
|
|
return; |
return; |
} |
} |
Sturm(inp, &sseq); |
Sturm(inp, &sseq); |
DUPQ(up,low); |
//DUPQ(up,low); |
SGN(low) = -1; |
//SGN(low) = -1; |
|
NEWQ(low); mpq_neg(BDY(low),BDY(up)); |
chvu = stumq(sseq, up); |
chvu = stumq(sseq, up); |
chvl = stumq(sseq, low); |
chvl = stumq(sseq, low); |
tnumb = abs(chvu - chvl); |
tnumb = abs(chvu - chvl); |
|
|
for (fctp = NEXT(fctp), i = 0; fctp; fctp = NEXT(fctp)) { |
for (fctp = NEXT(fctp), i = 0; fctp; fctp = NEXT(fctp)) { |
p = COEF(fctp); |
p = COEF(fctp); |
onedeg = DC(p); |
onedeg = DC(p); |
if ( !cmpq(DEG(onedeg), ONE) ) { |
if ( !cmpz(DEG(onedeg), ONE) ) { |
pad++; |
pad++; |
if ( !NEXT(onedeg) ) { |
if ( !NEXT(onedeg) ) { |
BDY(root)[pad][0] = 0; |
BDY(root)[pad][0] = 0; |
|
|
} |
} |
boundbody(p, &up); |
boundbody(p, &up); |
Sturm(p, &sseq); |
Sturm(p, &sseq); |
DUPQ(up,low); |
//DUPQ(up,low); |
SGN(low) = -1; |
//SGN(low) = -1; |
|
NEWQ(low); mpq_neg(BDY(low),BDY(up)); |
chvu = stumq(sseq, up); |
chvu = stumq(sseq, up); |
chvl = stumq(sseq, low); |
chvl = stumq(sseq, low); |
numb = abs(chvu - chvl); |
numb = abs(chvu - chvl); |
separate(DEG(fctp), eps, sseq, up, low, chvu, chvl, root, &pad); |
separate((Q)DEG(fctp), eps, sseq, up, low, chvu, chvl, root, &pad); |
} |
} |
for (i = 0; i < pad; i++) { |
for (i = 0; i < pad; i++) { |
for (j = i; j <= pad; j++) { |
for (j = i; j <= pad; j++) { |
|
|
int de, midn; |
int de, midn; |
Q mid, e; |
Q mid, e; |
P p; |
P p; |
|
Q tmp_e; |
|
|
de = abs(lown - upn); |
de = abs(lown - upn); |
if (de == 0) return; |
if (de == 0) return; |
|
|
BDY(root)[*padp][0] = up; |
BDY(root)[*padp][0] = up; |
BDY(root)[*padp][1] = low; |
BDY(root)[*padp][1] = low; |
BDY(root)[*padp][3] = (P *)sseq->body[0]; |
BDY(root)[*padp][3] = (P *)sseq->body[0]; |
subq( BDY(root)[*padp][1], BDY(root)[*padp][0], &e ); |
subq( BDY(root)[*padp][1], BDY(root)[*padp][0], &tmp_e ); |
SGN(e) = 1; |
//SGN(e) = 1; |
|
absq(tmp_e, &e); |
while (cmpq(e, eps) > 0) { |
while (cmpq(e, eps) > 0) { |
binary(*padp, root); |
binary(*padp, root); |
subq( BDY(root)[*padp][1], BDY(root)[*padp][0], &e); |
subq( BDY(root)[*padp][1], BDY(root)[*padp][0], &tmp_e); |
SGN(e) = 1; |
//SGN(e) = 1; |
|
absq(tmp_e, &e); |
} |
} |
BDY(root)[*padp][2] = mult; |
BDY(root)[*padp][2] = mult; |
return; |
return; |
} |
} |
addq(up, low, &mid); |
addq(up, low, &mid); |
divq(mid, TWO, &mid); |
divq(mid, (Q)TWO, &mid); |
midn = stumq(sseq, mid); |
midn = stumq(sseq, mid); |
separate(mult, eps, sseq, low, mid, lown, midn, root, padp); |
separate(mult, eps, sseq, low, mid, lown, midn, root, padp); |
separate(mult, eps, sseq, mid, up, midn, upn, root, padp); |
separate(mult, eps, sseq, mid, up, midn, upn, root, padp); |
|
|
P p; |
P p; |
p = (P)BDY(root)[indx][3]; |
p = (P)BDY(root)[indx][3]; |
addq(BDY(root)[indx][0], BDY(root)[indx][1], &c); |
addq(BDY(root)[indx][0], BDY(root)[indx][1], &c); |
divq(c, TWO, &d); |
divq(c, (Q)TWO, &d); |
ueval(p, BDY(root)[indx][1], &a); |
ueval(p, BDY(root)[indx][1], &a); |
ueval(p, d, &b); |
ueval(p, d, &b); |
if (SGN(a) == SGN(b)){ |
if (sgnq(a) == sgnq(b)){ |
BDY(root)[indx][1] = d; |
BDY(root)[indx][1] = d; |
} else { |
} else { |
BDY(root)[indx][0] = d; |
BDY(root)[indx][0] = d; |
|
|
v = VR(p); |
v = VR(p); |
t = (P *)ALLOCA((deg(v,p)+1)*sizeof(P)); |
t = (P *)ALLOCA((deg(v,p)+1)*sizeof(P)); |
g1 = t[0] = p; diffp(CO,p,v,(P *)&a); ptozp((P)a,1,&c,&g2); t[1] = g2; |
g1 = t[0] = p; diffp(CO,p,v,(P *)&a); ptozp((P)a,1,&c,&g2); t[1] = g2; |
for ( i = 1, h = ONE, x = ONE; ; ) { |
NEWQ(h); |
|
NEWQ(x); |
|
dupz(ONE, (Z *)&h); |
|
dupz(ONE, (Z *)&x); |
|
for ( i = 1; ; ) { |
if ( NUM(g2) ) break; |
if ( NUM(g2) ) break; |
subq(DEG(DC(g1)),DEG(DC(g2)),&d); |
subz(DEG(DC(g1)),DEG(DC(g2)),(Z*)&d); |
l = (Q)LC(g2); |
l = (Q)LC(g2); |
if ( SGN(l) < 0 ) { |
if ( sgnq(l) < 0 ) { |
chsgnq(l,&a); l = a; |
chsgnq(l,&a); l = a; |
} |
} |
addq(d,ONE,&a); pwrq(l,a,&b); mulp(CO,(P)b,g1,(P *)&a); |
addq(d,(Q)ONE,&a); pwrq(l,a,&b); mulp(CO,(P)b,g1,(P *)&a); |
divsrp(CO,(P)a,g2,&q,&r); |
divsrp(CO,(P)a,g2,&q,&r); |
if ( !r ) break; |
if ( !r ) break; |
chsgnp(r,&s); r = s; i++; |
chsgnp(r,&s); r = s; i++; |
|
|
pwrq(h,d,&m); g1 = g2; |
pwrq(h,d,&m); g1 = g2; |
mulq(m,x,&a); divsp(CO,r,(P)a,&g2); t[i] = g2; |
mulq(m,x,&a); divsp(CO,r,(P)a,&g2); t[i] = g2; |
x = (Q)LC(g1); |
x = (Q)LC(g1); |
if ( SGN(x) < 0 ) { |
if ( sgnq(x) < 0 ) { |
chsgnq(x,&a); x = a; |
chsgnq(x,&a); x = a; |
} |
} |
pwrq(x,d,&a); mulq(a,h,&b); divq(b,m,&h); |
pwrq(x,d,&a); mulq(a,h,&b); divq(b,m,&h); |
Line 370 stumq(VECT s, Q val) |
|
Line 383 stumq(VECT s, Q val) |
|
for ( i = j++, c =0; i < len; i++) { |
for ( i = j++, c =0; i < len; i++) { |
ueval( ss[i], val, &a); |
ueval( ss[i], val, &a); |
if ( a ) { |
if ( a ) { |
if( (SGN(a) > 0 && SGN(a0) < 0) || (SGN(a) < 0 && SGN(a0) > 0) ){ |
if( (sgnq(a) > 0 && sgnq(a0) < 0) || (sgnq(a) < 0 && sgnq(a0) > 0) ){ |
c++; |
c++; |
a0 = a; |
a0 = a; |
} |
} |
|
|
else if ( p->id == O_N ) |
else if ( p->id == O_N ) |
*q = 0; |
*q = 0; |
else { |
else { |
NEWQ(tmp); |
//NEWQ(tmp); |
SGN(tmp)=1; |
//SGN(tmp)=1; |
for ( dc = DC(p), max=0; dc; dc = NEXT(dc) ) { |
NEWQ(max); |
|
for ( dc = DC(p); dc; dc = NEXT(dc) ) { |
t = (Q)COEF(dc); |
t = (Q)COEF(dc); |
NM(tmp)=NM(t); |
if ( cmpq(t, max) > 0 ) { //DUPQ(tmp, max); |
DN(tmp)=DN(t); |
mpq_abs(BDY(max),BDY(t)); |
if ( cmpq(tmp, max) > 0 ) DUPQ(tmp, max); |
} |
} |
} |
addq(ONE, max, q); |
addq((Q)ONE, max, q); |
|
//mpq_clear(max); |
} |
} |
} |
} |
|
|
|
|
Q q; |
Q q; |
Q *rp; |
Q *rp; |
{ |
{ |
Q d, d1, a, b, t; |
Z d, d1; |
Q deg, da; |
Z deg; |
Q nm, dn; |
Q da; |
|
Q a, b, t; |
|
Z nm, dn; |
DCP dc; |
DCP dc; |
|
|
if ( !p ) *rp = 0; |
if ( !p ) *rp = 0; |
else if ( NUM(p) ) *rp = (Q)p; |
else if ( NUM(p) ) *rp = (Q)p; |
else { |
else { |
if ( q ) { |
if ( q ) { |
NTOQ( DN(q), 1, dn ); |
//NTOQ( DN(q), 1, dn ); |
NTOQ( NM(q), SGN(q), nm ); |
dnq(q, &dn); |
|
//NTOQ( NM(q), sgnq(q), nm ); |
|
nmq(q, &nm); |
} else { |
} else { |
dn = 0; |
dn = 0; |
nm = 0; |
nm = 0; |
|
|
if ( !dn ) { |
if ( !dn ) { |
dc = DC(p); t = (Q)COEF(dc); |
dc = DC(p); t = (Q)COEF(dc); |
for ( d = DEG(dc), dc = NEXT(dc); dc; d = DEG(dc), dc= NEXT(dc) ) { |
for ( d = DEG(dc), dc = NEXT(dc); dc; d = DEG(dc), dc= NEXT(dc) ) { |
subq(d, DEG(dc), &d1); pwrq(nm, d1, &a); |
subz(d, DEG(dc), &d1); pwrq((Q)nm, (Q)d1, &a); |
mulq(t,a,&b); addq(b,(Q)COEF(dc),&t); |
mulq(t,a,&b); addq(b,(Q)COEF(dc),&t); |
} |
} |
if ( d ) { |
if ( d ) { |
pwrq(nm,d,&a); mulq(t,a,&b); t = b; |
pwrq((Q)nm,(Q)d,&a); mulq(t,a,&b); t = b; |
} |
} |
*rp = t; |
*rp = t; |
} else { |
} else { |
dc = DC(p); t = (Q)COEF(dc); |
dc = DC(p); t = (Q)COEF(dc); |
for ( d=deg= DEG(dc), dc = NEXT(dc); dc; d = DEG(dc), dc= NEXT(dc) ) { |
for ( d=deg= DEG(dc), dc = NEXT(dc); dc; d = DEG(dc), dc= NEXT(dc) ) { |
subq(d, DEG(dc), &d1); pwrq(nm, d1, &a); |
subz(d, DEG(dc), &d1); pwrq((Q)nm, (Q)d1, &a); |
mulq(t,a,&b); |
mulq(t,a,&b); |
subq(deg, DEG(dc), &d1); pwrq(dn, d1, &a); |
subz(deg, DEG(dc), &d1); pwrq((Q)dn, (Q)d1, &a); |
mulq(a, (Q)COEF(dc), &da); |
mulq(a, (Q)COEF(dc), &da); |
addq(b,da,&t); |
addq(b,da,&t); |
} |
} |
if ( d ) { |
if ( d ) { |
pwrq(nm,d,&a); mulq(t,a,&b); t = b; |
pwrz(nm,d,(Z*)&a); mulq(t,a,&b); t = b; |
} |
} |
*rp = t; |
*rp = t; |
} |
} |
} |
} |
#endif |
|
} |
} |
#endif |
#endif |