用到了一元二次方程,一元三次方程的求解.
class QuarticRealPolynomial
{
public:
static Number computeDiscriminant(Number a, Number b, Number c, Number d, Number e);
static std::vector<Number> computeRealRoots(Number a, Number b, Number c, Number d, Number e);
private:
static std::vector<Number> original(Number a3, Number a2, Number a1, Number a0);
static std::vector<Number> neumark(Number a3, Number a2, Number a1, Number a0);
};
---------------------------------------------
Number QuarticRealPolynomial::computeDiscriminant(Number a, Number b, Number c, Number d, Number e)
{
Number a2 = a * a;
Number a3 = a2 * a;
Number b2 = b * b;
Number b3 = b2 * b;
Number c2 = c * c;
Number c3 = c2 * c;
Number d2 = d * d;
Number d3 = d2 * d;
Number e2 = e * e;
Number e3 = e2 * e;
Number discriminant = (b2 * c2 * d2 - 4.0 * b3 * d3 - 4.0 * a * c3 * d2 + 18 * a * b * c * d3 - 27.0 * a2 * d2 * d2 + 256.0 * a3 * e3) +
e * (18.0 * b3 * c * d - 4.0 * b2 * c3 + 16.0 * a * c2 * c2 - 80.0 * a * b * c2 * d - 6.0 * a * b2 * d2 + 144.0 * a2 * c * d2) +
e2 * (144.0 * a * b2 * c - 27.0 * b2 * b2 - 128.0 * a2 * c2 - 192.0 * a2 * b * d);
return discriminant;
}
std::vector<Number> QuarticRealPolynomial::computeRealRoots(Number a, Number b, Number c, Number d, Number e)
{
if (std::abs(a) < CesiumMath::_EPSILON15)
{
return CubicRealPolynomial::computeRealRoots(b, c, d, e);
}
Number a3 = b / a;
Number a2 = c / a;
Number a1 = d / a;
Number a0 = e / a;
int k = (a3 < 0.0) ? 1 : 0;
k += (a2 < 0.0) ? k + 1 : k;
k += (a1 < 0.0) ? k + 1 : k;
k += (a0 < 0.0) ? k + 1 : k;
switch (k)
{
case 0:
return original(a3, a2, a1, a0);
case 1:
return neumark(a3, a2, a1, a0);
case 2:
return neumark(a3, a2, a1, a0);
case 3:
return original(a3, a2, a1, a0);
case 4:
return original(a3, a2, a1, a0);
case 5:
return neumark(a3, a2, a1, a0);
case 6:
return original(a3, a2, a1, a0);
case 7:
return original(a3, a2, a1, a0);
case 8:
return neumark(a3, a2, a1, a0);
case 9:
return original(a3, a2, a1, a0);
case 10:
return original(a3, a2, a1, a0);
case 11:
return neumark(a3, a2, a1, a0);
case 12:
return original(a3, a2, a1, a0);
case 13:
return original(a3, a2, a1, a0);
case 14:
return original(a3, a2, a1, a0);
case 15:
return original(a3, a2, a1, a0);
default:
std::vector<Number> results;
return results;
}
}
std::vector<Number> QuarticRealPolynomial::original(Number a3, Number a2, Number a1, Number a0)
{
std::vector<Number> resultRoots;
Number a3Squared = a3 * a3;
Number p = a2 - 3.0 * a3Squared / 8.0;
Number q = a1 - a2 * a3 / 2.0 + a3Squared * a3 / 8.0;
Number r = a0 - a1 * a3 / 4.0 + a2 * a3Squared / 16.0 - 3.0 * a3Squared * a3Squared / 256.0;
// Find the roots of the cubic equations: h^6 + 2 p h^4 + (p^2 - 4 r) h^2 - q^2 = 0.
std::vector<Number> cubicRoots = CubicRealPolynomial::computeRealRoots(1.0, 2.0 * p, p * p - 4.0 * r, -q * q);
if (cubicRoots.size() > 0)
{
Number temp = -a3 / 4.0;
// Use the largest positive root.
Number hSquared = cubicRoots[cubicRoots.size() - 1];
if (std::abs(hSquared) < CesiumMath::_EPSILON14) {
// y^4 + p y^2 + r = 0.
std::vector<Number> roots;
QuadraticRealPolynomial::computeRealRoots(1.0, p, r, roots);
if (roots.size() == 2) {
Number root0 = roots[0];
Number root1 = roots[1];
Number y;
if (root0 >= 0.0 && root1 >= 0.0) {
Number y0 = std::sqrt(root0);
Number y1 = std::sqrt(root1);
resultRoots.push_back(temp - y1);
resultRoots.push_back(temp - y0);
resultRoots.push_back(temp + y0);
resultRoots.push_back(temp + y1);
return resultRoots;
} else if (root0 >= 0.0 && root1 < 0.0)
{
y = std::sqrt(root0);
resultRoots.push_back(temp - y);
resultRoots.push_back(temp + y);
return resultRoots;
} else if (root0 < 0.0 && root1 >= 0.0)
{
y = std::sqrt(root1);
resultRoots.push_back(temp - y);
resultRoots.push_back(temp + y);
return resultRoots;
}
}
return resultRoots;
} else if (hSquared > 0.0)
{
Number h = std::sqrt(hSquared);
Number m = (p + hSquared - q / h) / 2.0;
Number n = (p + hSquared + q / h) / 2.0;
// Now solve the two quadratic factors: (y^2 + h y + m)(y^2 - h y + n);
std::vector<Number> roots1;
QuadraticRealPolynomial::computeRealRoots(1.0, h, m, roots1);
std::vector<Number> roots2;
QuadraticRealPolynomial::computeRealRoots(1.0, -h, n, roots2);
if (roots1.size() != 0)
{
roots1[0] += temp;
roots1[1] += temp;
if (roots2.size() != 0)
{
roots2[0] += temp;
roots2[1] += temp;
if (roots1[1] <= roots2[0])
{
resultRoots.push_back(roots1[0]);
resultRoots.push_back(roots1[1]);
resultRoots.push_back(roots2[0]);
resultRoots.push_back(roots2[1]);
return resultRoots;
}
else if (roots2[1] <= roots1[0])
{
resultRoots.push_back(roots2[0]);
resultRoots.push_back(roots2[1]);
resultRoots.push_back(roots1[0]);
resultRoots.push_back(roots1[1]);
return resultRoots;
}
else if (roots1[0] >= roots2[0] && roots1[1] <= roots2[1])
{
resultRoots.push_back(roots2[0]);
resultRoots.push_back(roots1[0]);
resultRoots.push_back(roots1[1]);
resultRoots.push_back(roots2[1]);
return resultRoots;
}
else if (roots2[0] >= roots1[0] && roots2[1] <= roots1[1])
{
resultRoots.push_back(roots1[0]);
resultRoots.push_back(roots2[0]);
resultRoots.push_back(roots2[1]);
resultRoots.push_back(roots1[1]);
return resultRoots;
}
else if (roots1[0] > roots2[0] && roots1[0] < roots2[1])
{
resultRoots.push_back(roots2[0]);
resultRoots.push_back(roots1[0]);
resultRoots.push_back(roots2[1]);
resultRoots.push_back(roots1[1]);
return resultRoots;
}
resultRoots.push_back(roots1[0]);
resultRoots.push_back(roots2[0]);
resultRoots.push_back(roots1[1]);
resultRoots.push_back(roots2[1]);
return resultRoots;
}
return roots1;
}
if (roots2.size() != 0) {
roots2[0] += temp;
roots2[1] += temp;
return roots2;
}
resultRoots.clear();
return resultRoots;
}
}
resultRoots.clear();
return resultRoots;
}
std::vector<Number> QuarticRealPolynomial::neumark(Number a3, Number a2, Number a1, Number a0)
{
std::vector<Number> resultRoots;
Number a1Squared = a1 * a1;
Number a2Squared = a2 * a2;
Number a3Squared = a3 * a3;
Number p = -2.0 * a2;
Number q = a1 * a3 + a2Squared - 4.0 * a0;
Number r = a3Squared * a0 - a1 * a2 * a3 + a1Squared;
std::vector<Number> cubicRoots = CubicRealPolynomial::computeRealRoots(1.0, p, q, r);
if (cubicRoots.size() > 0) {
// Use the most positive root
Number y = cubicRoots[0];
Number temp = (a2 - y);
Number tempSquared = temp * temp;
Number g1 = a3 / 2.0;
Number h1 = temp / 2.0;
Number m = tempSquared - 4.0 * a0;
Number mError = tempSquared + 4.0 * std::abs(a0);
Number n = a3Squared - 4.0 * y;
Number nError = a3Squared + 4.0 * std::abs(y);
Number g2;
Number h2;
if (y < 0.0 || (m * nError < n * mError)) {
Number squareRootOfN = std::sqrt(n);
g2 = squareRootOfN / 2.0;
h2 = squareRootOfN == 0.0 ? 0.0 : (a3 * h1 - a1) / squareRootOfN;
} else {
Number squareRootOfM = std::sqrt(m);
g2 = squareRootOfM == 0.0 ? 0.0 : (a3 * h1 - a1) / squareRootOfM;
h2 = squareRootOfM / 2.0;
}
Number G;
Number g;
if (g1 == 0.0 && g2 == 0.0) {
G = 0.0;
g = 0.0;
} else if (CesiumMath::sign(g1) == CesiumMath::sign(g2)) {
G = g1 + g2;
g = y / G;
} else {
g = g1 - g2;
G = y / g;
}
Number H;
Number h;
if (h1 == 0.0 && h2 == 0.0) {
H = 0.0;
h = 0.0;
} else if (CesiumMath::sign(h1) == CesiumMath::sign(h2))
{
H = h1 + h2;
h = a0 / H;
} else {
h = h1 - h2;
H = a0 / h;
}
// Now solve the two quadratic factors: (y^2 + G y + H)(y^2 + g y + h);
std::vector<Number> roots1;
QuadraticRealPolynomial::computeRealRoots(1.0, G, H, roots1);
std::vector<Number> roots2;
QuadraticRealPolynomial::computeRealRoots(1.0, g, h, roots2);
if (roots1.size() != 0)
{
if (roots2.size() != 0) {
if (roots1[1] <= roots2[0])
{
resultRoots.push_back(roots1[0]);
resultRoots.push_back(roots1[1]);
resultRoots.push_back(roots2[0]);
resultRoots.push_back(roots2[1]);
return resultRoots;
} else if (roots2[1] <= roots1[0])
{
resultRoots.push_back(roots2[0]);
resultRoots.push_back(roots2[1]);
resultRoots.push_back(roots1[0]);
resultRoots.push_back(roots1[1]);
return resultRoots;
} else if (roots1[0] >= roots2[0] && roots1[1] <= roots2[1])
{
resultRoots.push_back(roots2[0]);
resultRoots.push_back(roots1[0]);
resultRoots.push_back(roots1[1]);
resultRoots.push_back(roots2[1]);
return resultRoots;
} else if (roots2[0] >= roots1[0] && roots2[1] <= roots1[1])
{
resultRoots.push_back(roots1[0]);
resultRoots.push_back(roots2[0]);
resultRoots.push_back(roots2[1]);
resultRoots.push_back(roots1[1]);
return resultRoots;
} else if (roots1[0] > roots2[0] && roots1[0] < roots2[1])
{
resultRoots.push_back(roots2[0]);
resultRoots.push_back(roots1[0]);
resultRoots.push_back(roots2[1]);
resultRoots.push_back(roots1[1]);
return resultRoots;
} else
{
resultRoots.push_back(roots1[0]);
resultRoots.push_back(roots2[0]);
resultRoots.push_back(roots1[1]);
resultRoots.push_back(roots2[1]);
return resultRoots;
}
}
return roots1;
}
if (roots2.size() != 0)
{
return roots2;
}
}
resultRoots.clear();
return resultRoots;
}