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
|
#include "CubicUtilities.h"
#include "DataTypes.h"
#include "QuadraticUtilities.h"
void coefficients(const double* cubic, double& A, double& B, double& C, double& D) {
A = cubic[6]; // d
B = cubic[4] * 3; // 3*c
C = cubic[2] * 3; // 3*b
D = cubic[0]; // a
A -= D - C + B; // A = -a + 3*b - 3*c + d
B += 3 * D - 2 * C; // B = 3*a - 6*b + 3*c
C -= 3 * D; // C = -3*a + 3*b
}
// cubic roots
const double PI = 4 * atan(1);
static bool is_unit_interval(double x) {
return x > 0 && x < 1;
}
// from SkGeometry.cpp (and Numeric Solutions, 5.6)
int cubicRoots(double A, double B, double C, double D, double t[3]) {
if (approximately_zero(A)) { // we're just a quadratic
return quadraticRoots(B, C, D, t);
}
double a, b, c;
{
double invA = 1 / A;
a = B * invA;
b = C * invA;
c = D * invA;
}
double a2 = a * a;
double Q = (a2 - b * 3) / 9;
double R = (2 * a2 * a - 9 * a * b + 27 * c) / 54;
double Q3 = Q * Q * Q;
double R2MinusQ3 = R * R - Q3;
double adiv3 = a / 3;
double* roots = t;
double r;
if (R2MinusQ3 < 0) // we have 3 real roots
{
double theta = acos(R / sqrt(Q3));
double neg2RootQ = -2 * sqrt(Q);
r = neg2RootQ * cos(theta / 3) - adiv3;
if (is_unit_interval(r))
*roots++ = r;
r = neg2RootQ * cos((theta + 2 * PI) / 3) - adiv3;
if (is_unit_interval(r))
*roots++ = r;
r = neg2RootQ * cos((theta - 2 * PI) / 3) - adiv3;
if (is_unit_interval(r))
*roots++ = r;
}
else // we have 1 real root
{
double A = fabs(R) + sqrt(R2MinusQ3);
A = cube_root(A);
if (R > 0) {
A = -A;
}
if (A != 0) {
A += Q / A;
}
r = A - adiv3;
if (is_unit_interval(r))
*roots++ = r;
}
return (int)(roots - t);
}
|
