一元四次方程求解C++实现
用到了一元二次方程,一元三次方程的求解.
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;
}
一元四次方程求解C++实现相关推荐
- Fast Planner——KinodynamicAstar::estimateHeuristic()中一元三次方程和一元四次方程求解
Fast Planner的代码中,前端路径搜索时需要评估路径的启发成本(函数KinodynamicAstar::estimateHeuristic),涉及到一元四次方程和一元三次方程的求解计算.Kin ...
- python求一元三次方程的根_关于二次、三次、四次方程求解方法讨论
高次方程求解的一般方法是将高次方程通过配方求解,然后进行次数降解,高次方程转化为容易求解的低次方程. 一元二次方程 求解高次方程,一元二次方程是最为简单的方程.关于一元二次方程 ,通过配方法可以求解: ...
- logit方程怎么写_一元四次方程的常规解法
第一次写知乎文章还是有点小激动的,我不熟悉公式编辑,就用我的卑微MathType好了 这篇文章初中生也可以听懂的! 正文开始 我们考虑标准一元四次方程 这里a≠0,我们第一个想到的应该是配方法,我们令 ...
- 一元三次方程重根判别式_一元四次方程的常规解法
第一次写知乎文章还是有点小激动的,我不熟悉公式编辑,就用我的卑微MathType好了 这篇文章初中生也可以听懂的! 正文开始 我们考虑标准一元四次方程 这里a≠0,我们第一个想到的应该是配方法,我们令 ...
- 问题三十七：C++怎么解一元四次方程？（3）——怎么解一元四次方程
37.3 怎么解一元四次方程? 用"费拉里"方法求解:将四次方程化为两个二次方程,然后求解二次方程. --------------------------------------- ...
- 一元四次方程求根公式
求解一元四次方程 VC++ 代码链接如下: https://github.com/hanford77/SolveEquation
- zcmu-2116一元三次方程求解
2116: 一元三次方程求解 Time Limit: 1 Sec Memory Limit: 128 MB Submit: 65 Solved: 23 [Submit][Status][Web B ...
- 1814: 一元三次方程求解
//很久之前写的,记录一下~ 1814: 一元三次方程求解 Time Limit: 1 Sec Memory Limit: 128 MB Submit: 45 Solved: 28 [Submit][ ...
- 【luogu 1024 一元三次方程求解】二分思想
题目出自luogu 1024 一元三次方程求解 描述: 有形如:ax3+bx2+cx+d=0 这样的一个一元三次方程.给出该方程中各项的系数(a,b,c,d 均为实数),并约定该方程存在三个不同实根( ...
最新文章
- 异常org.xml.sax.SAXParseException; lineNumber: 5; columnNumber: 11; 注释中不允许出现字符串 --。的原因...
- C# Winform程序中使用TeeChart实现简单的图表展示
- Laravel 实现定时任务
- 工作147:外部that
- python字符串垂直输出加循环_将漂亮的soup嵌套循环垂直输出到datafram中
- paper reading：[第一代GCN] Spectral Networks and Deep Locally Connected Networks on Graphs
- Spring JDK动态代理详解
- Java实现微信公众号授权登录
- 普通学历，大一大二要不要打ACM？
- 一起认识国产又好用的uni-app
- 阿里和CVTE秋招面试题
- 【GNSS】gfzrnx-用法
- dB、dbm、dbw、W 相互关系
- 在vue中实现picker样式_vue中van-picker的多列联动数据格式如何设置以及调用
- 能挽救这条船的,唯有你图片_这是科技如何帮助普通民众挽救生命
- context capture如进行空三迁移
- 最短增广路Isap算法 网络流
- 疫情期北京融资信息分析---疫情对北京社会经济影响分析---科技战疫·大数据公益挑战赛---2020北京数据开放创新应用大赛
- 统计推断——假设检验——两变量关联性分析
- Nginx - 记一次Nginx端口转发失败案例