拉格朗日插值法（c++实现，生成的多项式的系数由分数表示）

最新推荐文章于 2024-02-11 23:21:46 发布

无畏的气球兵

最新推荐文章于 2024-02-11 23:21:46 发布

阅读量876

点赞数 23

文章标签： c++ 数学建模学习

本文链接：https://blog.csdn.net/qq_62604677/article/details/135434707

版权

前言

今天清理电脑空间，无意间发现很久之前写的一个拉格朗日插值法的代码。现在回过头看，当时写的纯纯就是s**t，不过还是有一定纪念价值，所以写个文章纪念一下。

我依稀记得，事情的起因是知乎上刷到了这样一个问题：

有争议的数字找规律题，数字4、7、11、18、29、47、（？），请问括号内应填几？ - 知乎

下面的回答清一色是这样的：

这个眼花缭乱的多项式在当初那个只学了工科阉割版线代并且没学数学建模的我看来，非常恐怖，而且我想了很久也没想明白这是怎么得出来的，好在有些比较“良心”的回答提到了这个多项式是用拉格朗日插值法得出的，于是我学习了一下，并写出了这个代码。

本来是不想自己写的，结果我在网上搜到的代码都是直接用浮点数计算出结果，没有把系数用分数表示这么“炫酷”，索性自己写了一个。

由于当时刚学完c++，不太熟悉，所以写的很烂，而且核心代码（就是实现拉格朗日插值的那部分），当时也没写注释，我看了半天实在看不懂，所以就稍微加了一点注释直接发出来了。

拉格朗日插值法介绍

我也懒得自己写了，反正也不难，只贴两个链接。

如何直观地理解拉格朗日插值法？ - 知乎

实现代码

#include <bits/stdc++.h>
typedef unsigned long long ull;
using namespace std;

vector<int> expansion(vector<int> a){
    //将形如(x-a0)(x-a1)...(x-an)的式子展开为x^n+b1x^(n-1)+b2x^(n-2)...+bn的形式，算法没仔细想，可能时间复杂度高
    ull len=a.size();
    vector<int> ret(len+1,0);

    for(int i=0;i<len-1;++i){
        if(i==0){
            ret[len]=a[0]*a[1];
            ret[len-1]=a[0]+a[1];
            ret[len-2]=1;
        }else{
            for(ull j=len-2-i;j<=len;++j){
                if(j==len){
                    ret[j]*=a[i+1];
                }else{
                    ret[j]*=a[i+1];
                    ret[j]+=ret[j+1];
                }
            }
        }
    }
    return ret;
}
int gcd(int a,int b){
    //求a，b最大公约数
    if(a>=b){
        if(a%b==0){
            return b;
        }
        int c=1;
        while(c){
            c=a%b;
            a=b;
            b=c;
        }
    }else{
        return gcd(b,a);
    }
    return a;
}

class fraction{
    //用来实现分数的类 
private:
    int numerator;           //分子
    int denominator;         //分母
public:
    fraction(int n=0,int d=1){
    	//构造函数，根据输入自动约分 
        //默认为0
        numerator=n;
        denominator=d;
        if(numerator==0) return;
        //用最大公约数实现约分
        int g=gcd(numerator,denominator);
        numerator/=g;
        denominator/=g;
    }
    fraction(const fraction &b){
        //拷贝构造函数
        //其实这个函数根本没必要，可能当时不太熟悉c++吧 
        this->numerator=b.numerator;
        this->denominator=b.denominator;
    }
    void set_value(int n,int d){
		// 更改这个分数的值，并根据输入自动约分
        numerator=n;
        denominator=d;
        if(numerator==0) return;
        //用最大公约数实现约分
        int g=gcd(numerator,denominator);
        numerator/=g;
        denominator/=g;
    }
    fraction operator+(const fraction &b){
    	//重载了加号运算符  
        int lcm=this->denominator*b.denominator/gcd(this->denominator,b.denominator);
        int num=this->numerator*lcm/ this->denominator+b.numerator*lcm/b.denominator;
        int den=lcm;
        return fraction(num,den);
    }
    fraction multiply(const int &b){
    	
    	//相乘并自动约分 
        if(b==0){
            this->numerator=0;
            return *this;
        }
        this->numerator= this->numerator*b;
        int g=gcd(this->numerator,this->denominator);
        this->numerator/=g;
        this->denominator/=g;
    }
    bool judge() {
    	//判断该分数正负 
        if((numerator>=0 && denominator>=0) || (numerator<0 && denominator<0)) return true;
        else return false;
    }
    void display(){
    	// 打印这个分数 
        if(numerator==0) cout<<0;
        else if(denominator==0) cout<<numerator<<'/'<<denominator;
        else if(denominator==1) cout<<numerator/denominator;
        else if(denominator<0) cout<<'-'<<numerator<<'/'<<-denominator;
        else cout<<numerator<<'/'<<denominator;
    }

};

int main() {
    int n;
    cout<<"请输入n的大小："<<endl;
    cin>>n;
    vector<int> a(n);
    vector<int> b(n);
// 输入的这两个向量a，b分别代表：f(a[i]) == b[i] (0 <= i <= n)
    cout<<"请输入向量a的n个元素：(直接输入-1使这向量n为{1,2,3。。。n}）"<<endl;
    for(int i=0;i<n;++i){
        cin>>a[i];
        if(i==0 && a[0]==-1){
            for(int j=0;j<n;++j){
                a[j]=j+1;
            }
            break;
        }
    }
    cout<<"请输入向量b的n个元素："<<endl;
    for(int i=0;i<n;++i){
        cin>>b[i];
    }
    vector<fraction> l(n);
    fraction temp(0,1);
    vector<fraction> f(n,temp);
    for(int i=0;i<n;++i){
        vector<int> temp1;
        int den=1;
        for(int j=0;j<n;++j){
            if(j!=i){
                temp1.push_back(-a[j]);
                den*=(a[i]-a[j]);
            }
        }
        //最后temp1有n-1个元素
        vector<int> temp2(expansion(temp1));
        //temp2有n个元素
        for(int j=0;j<n;++j){

            l[j].set_value(temp2[j],den);
            l[j].multiply(b[i]);
            f[j]=f[j]+l[j];
        }
    }
    cout<<"函数为：f(x)=";
    for(int i=0;i<n;++i){
        f[i].display();
        if(i==n-1){
            cout<<endl;
            break;
        }
        if(f[i+1].judge()){
            cout<<"x^"<<n-i-1<<'+';
        }else{
            cout<<"x^"<<n-i-1;
        }

    }

    return 0;
}