...in any case, here is a quick and dirty fix for the quartic solver. During the bug hunt I changed almost everything, so it is faster now than the original was, and solved all the examples I fed into it. Among 1 million random equations solved (took less than an hour) there was 1 with larger max error than |maxroot|*1.e-6 (~1.3e-5), and 3 with relative errors around 1.e-7. All others solutions were more accurate.
Code:
cbrt(x) { ; signed real cubic root
Static z = 0.33333333333333333
Return x<0 ? -(-x)**z : x**z
}
Max3(a,b,c) {
Return a < b ? (b > c ? b : c) : (a > c ? a : c)
}
;----------------------------------------------------------------
; C-idea http://van-der-waals.pc.uni-koeln.de/quartic/quintic_C.c
;----------------------------------------------------------------
cubicX(A3,A2,A1,A0, ByRef X) { ; X <- Max Real root of cubic, require X >= 0
Static eps = 1.e-99, o6 = .16666666666666667, p23 = 2.09439510239319549 ; 2pi/3
if (A3 != 0) { ; cubic
W := A2/A3/3
P := A1/A3/3 - W*W, P *= -P*P ; make P >= 0
Q := (A1*W-A0)/A3/2 - W*W*W
DIS := Q*Q - P
if (DIS >= 0) ; DIS > 0: one real solution
D := sqrt(DIS), X := cbrt(Q+D) + cbrt(Q-D) - W
if (DIS < 0 || X < 0) { ; three real solutions, X<0: ROUNDING ERROR!!!
PHI := acos(Q/sqrt(P)) / 3
X := 2*P**o6 * Max3( cos(PHI), cos(PHI+p23), cos(PHI+2*p23)) - W
}
}
else if (A2 != 0) ; QUADRATIC
P := A1/A2/2, X := sqrt(abs(P*P-A0/A2)) - P
else ; LINEAR equation
X := A0/A1
X := X < 0 ? 0 : X ; X < 0: ROUNDING ERROR!!!
y := A1+X*(2*A2+X*3*A3) ; p'(X) = denominator in Newton iteration
z := A0+X*(A1+X*(A2+X*A3)) ; p(X)
U := X - z/y ; Newton iteration, if reduces rounding errors
if (abs(y) > eps && abs(A0+U*(A1+U*(A2+U*A3))) < abs(z) && U >= 0)
X := U
}
quartic(a,b,c,d,f, ByRef x0,ByRef xi0, ByRef x1,ByRef xi1, ByRef x2,ByRef xi2, ByRef x3,ByRef xi3, ByRef Nsol) {
Static eps = 1.e-99
Nsol := 0
b := b/a/4, c := c/a, d := d/a, f := f/a, bb := b*b
p :=-6*bb + c
q := 8*b*bb - 2*b*c + d
r :=-3*bb*bb + bb*c - b*d + f
; biquadratic u^4 + p*u^2 + r = 0: CHEAT (small error: "a" normalized to 1)
q := abs(q)>1.e-13 ? q : q<0 ? -1.e-13 : 1.e-13
cubicX(1,2*p,p*p-4*r,-q*q, z0) ; Max real root z0 >= 0
z0 := z0<eps ? eps : z0 ; NO division by 0 !
xK := sqrt(z0)
xL := q/2/xK
sqm := 4*xL - 2*p - z0
sqp := sqm - 8*xL
xi0 := xi1 := xi2 := xi3 := 0.0
if (sqp >= 0 && sqm >= 0) {
x0 := ( xK + sqrt(sqp))/2
x1 := xK - x0
x2 := (-xK + sqrt(sqm))/2
x3 := -xK - x2
Nsol := 4
}
else if (sqp >= 0 && sqm < 0) {
x0 := (xK + sqrt(sqp))/2
x1 := xK - x0
x2 := x3 := -xK/2
xi3:=-xi2:= sqrt(-sqm/4)
Nsol := 2
}
else if (sqp < 0 && sqm >= 0) {
x0 := (sqrt(sqm)-xK)/2
x1 := -x0-xK
x2 := x3 := xK/2
xi3:=-xi2:= sqrt(-sqp/4)
Nsol := 2
}
else if (sqp < 0 && sqm < 0) {
x2 := x3 := xK/2
x0 := x1 := -x2
xi1:=-xi0:= sqrt(-sqm/4)
xi3:=-xi2:= sqrt(-sqp/4)
Nsol := 0
}
x0 -= b, x1 -= b, x2 -= b, x3 -= b
}
Login
Register
FAQ
Search
