Repository navigation
Expand file tree
/
Copy pathcubic.go
More file actions
81 lines (75 loc) · 2.56 KB
/
Copy pathcubic.go
File metadata and controls
81 lines (75 loc) · 2.56 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
// SPDX-FileCopyrightText: 2025 Dominik Honnef and contributors
//
// SPDX-License-Identifier: MIT
package polyroot
import (
"math"
)
// cubicFirstRootBetween is a special case of [rootsBetweenWithBuffer] that only
// returns the first root it finds. Additionally, for unbounded ranges, it
// applies logic specific to cubics to find a bound when finding the roots of
// the derivative fails.
func (poly *Polynomial) cubicFirstRootBetween(deriv *Polynomial, x0, x1, xError float64) (float64, bool) {
if r1, r2, n := deriv.quadraticRoots(); n > 0 {
xa := r1
xb := xa
if n > 1 {
xb = r2
}
if math.IsInf(x0, 0) || math.IsInf(x1, 0) {
// FIXME(dh): this discards one of x0 or x1 if only one of them is inf
ya := poly.Evaluate(xa)
if math.Signbit(poly.coeffs[3]) != math.Signbit(ya) {
yb := poly.Evaluate(xb)
return findRoot(poly, deriv, xb, math.Inf(1), yb, xError), true
} else {
return findRoot(poly, deriv, math.Inf(-1), xa, ya, xError), true
}
} else {
// This loop expands into the same kind of control flow as the highly
// nested implementation of CubicRoots in Yuksel's code, except that it
// executes a handful of superfluous comparisons. It is, however, much
// more readable.
possibleEndpoints := [3]float64{xa, xb, x1}
xl := x0
yl := poly.Evaluate(x0)
for _, xu := range possibleEndpoints {
if xu > xl && xu <= x1 {
yu := poly.Evaluate(xu)
if math.Signbit(yl) != math.Signbit(yu) {
return findRoot(poly, deriv, xl, xu, yl, xError), true
}
xl = xu
yl = yu
}
}
}
} else {
// If the derivative doesn't have any roots, then the cubic has an
// inflection point that isn't a stationary point.
if math.IsInf(x0, 0) || math.IsInf(x1, 0) {
// The roots of a cubic end with a -c/3d term, so is the idea here
// that the roots cannot be smaller or larger than that term,
// depending on the signs involved, and that's how we find one of
// the missing bounds(?)
xInf := -poly.coeffs[2] / (poly.coeffs[3] * 3)
yInf := poly.Evaluate(xInf)
if !isFin(xInf) {
return 0, false
}
// FIXME(dh): handle cases where only one of x0 or x1 is an infinity
if math.Signbit(poly.coeffs[3]) != math.Signbit(yInf) {
return findRoot(poly, deriv, xInf, math.Inf(1), yInf, xError), true
} else {
return findRoot(poly, deriv, math.Inf(-1), xInf, yInf, xError), true
}
} else {
y0 := poly.Evaluate(x0)
y1 := poly.Evaluate(x1)
if math.Signbit(y0) != math.Signbit(y1) {
return findRoot(poly, deriv, x0, x1, y0, xError), true
}
}
}
return 0, false
}