基于二分和倍增的一元n次方程求近似实根算法
一元n次方程求根算法
update
2025.3.15 贴上了时间复杂度估算、代码和数据
2025.3.16 发现算法不能解决单调奇次幂函数，已修复
2025.5.19 学神经网络梯度下降的时候发现一阶导数为0的点不一定是函数驻点（没正经学数学导致的），如\(f(x)=x^3\)在\(x=0\)时一阶导数为0，但那是拐点不是驻点，增加单调性判断以保证为驻点
2025.5.20 优化驻点判断方案，判断算法时间复杂度降为\(O(1)\)
前言
这是我于2025.3.11在椅子上躺着玩手机突然想到的，由于思路过于简单感觉几百年前就有人想到了，当deepseek找类似文章没找到，自己又不会找，并且想记录一下，于是有了这篇博客
该算法有很大的精度缺陷，且没有时间复杂度的严格推导（因为我不会），要是能提出改进建议或者彻底否定我将感激不尽
算法内容
根据阿贝尔-鲁菲尼定理，五次及以上的方程没有求根公式，换句话说不存在\(O(1)\)时间复杂度的通解算法，因此考虑多项式复杂度
对于n元一次方程\(\sum_{i=0}^{i≤n}k_ix^i=0\)（其中\(k_i\)为方程\(i\)次项系数），它的根（最多有\(n\)个）为函数\(f(x)=\sum_{i=0}^{i≤n}k_ix^i\)与\(x\)轴的交点横坐标
函数分为两种考虑：
1.无驻点
2.有驻点
对于无驻点的函数，显然函数整体满足单调性，直接从0点开始暴力倍增二分即可，只有一个根
对于有驻点的函数，易发现函数与\(x\)轴交点分三种情况：
1.在函数的驻点上
2.在函数两个相邻的驻点之间
3.在函数最左（右）的驻点的左（右）边
假设我们现在知道函数的每一个驻点，三种情况下根就很容易求出来了。
对于第一种情况，设该驻点坐标为\((x,f(x))\)，显然当\(|f(x)|<eps\)（这里\(eps\)为最大能接受的精度误差）时\(x\)即为其中一个根
对于第二种情况，根据驻点的定义，两相邻的驻点之间的函数一定满足单调性，因此容易想到二分
对于第三种情况，以最左边的驻点为例，由于左边不存在驻点，因此一定满足单调性。因此我们可以先通过从小到大倍增确定二分的左端点，再进行二分
那么，怎么求出每一个驻点呢？
众所周知，函数的驻点处导数为0，而\(f'(x)=0\)，不就是对导函数对应的方程进行求根吗？不过学过微积分的都知道一阶导数为0的点不一定是驻点，不过我们可以发现原函数的驻点两边单调性发生变化，对应到一阶导函数，该点左右正负性一定不同。也就是说当且仅当一阶导数为0且该点不为一阶导函数驻点时该点为函数的驻点，因此可以在求导函数根的时候舍去与驻点重合的根（注意原函数求根时不要舍去 ）。导函数可以在\(O(n)\)时间复杂度求出，每一次求导，函数都会降一次幂，直到最终为一次函数，可以通过一次函数求根公式在\(O(1)\)时间复杂度求出，将其返回即是上一层递归计算所需的函数顶点。
算法时间复杂度估算
假设要解的是\(n\)次方程，那么要向下递归\(n\)次，每次都要循环\(n\)次，在循环中又套了个二分，然后二分的过程中要实时求函数值，设\(X\)为最大的相邻驻点的\(y\)坐标差，那么时间复杂度差不多就是\(O(n^3logX)\)的规模吧
听起来挺吓人但在极大数据下TLE之前已经炸double了（难道真有算一千次方程的吗）
code(史)
细节很多，代码比较史
列一下要注意的地方：
1.必须按照从左到右的顺序求根，要不然求完还得排序
2.变量lb1,lb2,rb1,rb2必须先赋值false，具体原因注释写了
3.注意去掉重根
4.导函数无实数根不一定原函数无实数根，判断无实数根时加个判断是原函数的条件，如果导函数无实数根说明原函数是单调的
#include <bits/stdc++.h>
//#define int int64_t
//#define int __int128
//#define MOD (1000000007)
#define endl '\n'
#define debug_endl cout<<endl;
#define debug cout<<"debug"<<endl;
#define deltaX2 (0)//我也不知道这个东西有没有必要总之设了个这么个参数防止二分死循环，但由于eps的存在似乎不必要所以设置为0了
using namespace std;
typedef long double LD;
vector<LD> k;
int n;
LD deltaX;
LD eps;
LD f(LD x,vector<LD> tk){//求函数值
LD res=0;
for(int i=0;i<(int)tk.size();++i){
res+=tk[i]*pow(x,i);
}
return res;
}
vector<LD> solve(vector<LD> tk){//系数为当前要求解的方程系数
vector<LD> res;//方程的根
if(tk.size()==2){//x=-b/k，一次方程
res.emplace_back(-tk[0]/tk[1]);
return res;
}
vector<LD> nexk;//导函数系数
for(int i=1;i<(int)tk.size();++i){//根据幂法则求导，常数项舍去从1开始
nexk.emplace_back(tk[i]*i);
}
vector<LD> x=solve(nexk);//向下递归求导函数方程根，获得函数驻点
if(x.empty()){//该函数无驻点，单调，那肯定有根，直接倍增
LD fz=f(0,tk),df=f(deltaX,tk);
if(fabsl(fz)<eps){
res.emplace_back(0);
}
else if(df>fz){//单调递增
LD l,r;
if(fz>0){//往左找
r=0;
for(int i=1;;i<<=1){//从小到大倍增左端点
if(f(-i,tk)<0){//直到左端点与驻点在x轴异侧
l=-LD(i);
break;
}
}
}
else{
l=0;
for(int i=1;;i<<=1){//从小到大倍增右端点
if(f(i,tk)>0){//直到右端点与驻点在x轴异侧
r=LD(i);
break;
}
}
}
LD mid;
while(fabsl(r-l)>=eps){//二分求根
mid=(l+r)/LD(2.0);
if(fabsl(f(mid,tk))<eps){
break;
}
if(f(mid,tk)>0){
r=mid-deltaX2;
}
else{
l=mid+deltaX2;
}
}
res.emplace_back(mid);
}
else{//单调递减
LD l,r;
if(fz<0){//往左找
r=0;
for(int i=1;;i<<=1){//从小到大倍增左端点
if(f(-i,tk)>0){//直到左端点与驻点在x轴异侧
l=-LD(i);
break;
}
}
}
else{
l=0;
for(int i=1;;i<<=1){//从小到大倍增右端点
if(f(i,tk)<0){//直到右端点与驻点在x轴异侧
r=LD(i);
break;
}
}
}
LD mid;
while(fabsl(r-l)>=eps){//二分求根
mid=(l+r)/LD(2.0);
if(fabsl(f(mid,tk))<eps){
break;
}
if(f(mid,tk)>0){
l=mid+deltaX2;
}
else{
r=mid-deltaX2;
}
}
res.emplace_back(mid);
}
return res;
}
//先处理左边界驻点
LD lft=x.front();
LD flft=f(lft,tk);
bool lb1=false,lb2=false;//lb1:当前是单调递增 lb2:当前是单调递减
//赋值是必要的，由于c++的优化机制，或运算如果前面的=1直接跳过后半部分，也就是说lb2可能没有赋值，如果直接调用返回的是个整数，直接完蛋
if(fabsl(flft)<eps&&int(tk.size())==n+1){//驻点就是根且当前为原方程求解
if(res.empty()||lft!=res.back())//防止重根
res.emplace_back(lft);
}
else if((lb1=(flft>0&&f(lft-deltaX,tk)<flft))||(lb2=(flft<0&&f(lft-deltaX,tk)>flft))){//如果驻点在x轴上方，开口必须向下才有根，反之亦然
LD l,r=lft;
for(int i=1;;i<<=1){//从小到大倍增左端点
if((lb1&&f(lft-i,tk)<0)||(lb2&&f(lft-i,tk)>0)){//直到左端点与驻点在x轴异侧
l=lft-LD(i);
break;
}
}
LD mid;
while(fabsl(r-l)>=eps){//二分求根
mid=(l+r)/LD(2.0);
if(fabsl(f(mid,tk))<eps){
break;
}
if(f(mid,tk)>0){
(lb1?r=mid-deltaX2:l=mid+deltaX2);//画图很好理解
}
else{
(lb1?l=mid+deltaX2:r=mid-deltaX2);
}
}
if(res.empty()||mid!=res.back())//防止重根
res.emplace_back(mid);
}
//处理中间的
for(int i=0;i<(int)x.size()-1;++i){
LD fl=f(x[i],tk),fr=f(x[i+1],tk);
if(fabsl(fl)<eps&&int(tk.size())==n+1){//驻点就是根
if(res.empty()||fl!=res.back())//防止重根
res.emplace_back(fl);
continue;
}
if((fl>0&&fr>0)||(fl<0&&fr<0)){//两个驻点都在x轴同侧，不存在根
continue;
}
bool b=(fl<fr);//单调递增还是单调递减
LD l=x[i],r=x[i+1];
LD mid;
while(fabsl(r-l)>=eps){//二分求根
mid=(l+r)/LD(2.0);
if(fabsl(f(mid,tk))<eps){
break;
}
if(f(mid,tk)>0){
(b?r=mid-deltaX2:l=mid+deltaX2);
}
else{
(b?l=mid+deltaX2:r=mid-deltaX2);
}
}
if(res.empty()||mid!=res.back())//防止重根
res.emplace_back(mid);
}
//最后处理右边界驻点
LD rit=x.back();
LD frit=f(rit,tk);
bool rb1=false,rb2=false;
if(fabsl(frit)<eps&&int(tk.size())==n+1){//驻点就是根
if(res.empty()||rit!=res.back())//防止重根
res.emplace_back(rit);
}
else if((rb2=(frit>0&&f(rit+deltaX,tk)<frit))||(rb1=(frit<0&&f(rit+deltaX,tk)>frit))){
LD l=rit,r;
for(int i=1;;i<<=1){//从小到大倍增右端点
if((rb1&&f(rit+i,tk)>0)||(rb2&&f(rit+i,tk)<0)){//直到右端点与驻点在x轴异侧
r=rit+LD(i);
break;
}
}
LD mid;
while(fabsl(r-l)>=eps){//二分求根
mid=(l+r)/LD(2.0);
if(fabsl(f(mid,tk))<eps){
break;
}
if(f(mid,tk)>0){
(rb1?r=mid-deltaX2:l=mid+deltaX2);
}
else{
(rb1?l=mid+deltaX2:r=mid-deltaX2);
}
}
if(res.empty()||mid!=res.back())//防止重根
res.emplace_back(mid);
}
if(res.empty()&&int(tk.size())==n+1){//无实数根
cout<<"无实数根";
exit(0);
}
return res;
}
signed main(){
cout<<"输入精度（eps）：";
cin>>eps;
cout<<"输入用于判断单调性的微小增量△x：";
cin>>deltaX;
cin>>n;
for(int i=0;i<=n;++i){//从n次项到0次项依次输入系数
LD x;
cin>>x;
k.emplace_back(x);
}
reverse(k.begin(),k.end());//方便计算次幂和下标对齐
vector<LD> ans=solve(k);
for(LD answer:ans){
cout<<answer<<' ';
}
return 0;
}
验证数据
对于无驻点的单调函数，AT_past202203_g的样例是个不错的hack数据（事实上我就是从那发现算法的不足的）
期望输出:\(1.5\)
对于有驻点的函数，来一组四次无理方程求解，这是我自己因式分解乘出来构造的，\(eps\)精度差不多\(1^{-10}\)就行,\(\Delta x\)\(0.001\)即可
期望输出::
运行结果是正确的，图我就不贴了懒，可以自己运行试一下
浙公网安备 33010602011771号