Math/PolynomialRestore_double.hpp

• View this file on GitHub
• Last update: 2023-04-05 23:10:22+09:00
• Include: #include "Math/PolynomialRestore_double.hpp"

Depends on

• Math/LinearEquation_double.hpp

Code

#include "./LinearEquation_double.hpp"

// f(x) = a_0 + a_1 * x .. .a_d * x^d
// O(d ^ 3)
template<typename T>
vector<T> PolynomialRestore(vector<T> x,vector<T> y,int d,T EPS=1e-6){
    vector<vector<T>> M(x.size(),vector<T>(d+1));
    for(int i=0;i<(int)x.size();i++){
        T p=1;
        for(int j=0;j<=d;j++){
            M[i][j]=p;
            p*=x[i];
        }
    }
    auto ret=LinearEquation(M,y,EPS);
    return ret;
}

#line 1 "Math/LinearEquation_double.hpp"
// solve A * X = Y
// T : floating point
template<typename T>
vector<T> LinearEquation(vector<vector<T>> A,vector<T> Y,T EPS=1e-6){
    int n=(int)A.size(),m=(int)A[0].size();
    for(int i=0;i<n;i++) A[i].push_back(Y[i]);

    int rank=0;
    for(int j=0;j<m;j++){
        int p=-1;
        T mx=EPS;
        for(int i=rank;i<n;i++)if(chmax(mx,abs(A[i][j]))) p=i;
        if(p==-1) continue;

        swap(A[p],A[rank]);
        T t=A[rank][j];
        for(int k=0;k<=m;k++) A[rank][k]/=t;
        for(int i=0;i<n;i++){
            if(i==rank) continue;
            T s=A[i][j];
            for(int k=0;k<=m;k++) A[i][k]-=s*A[rank][k];
        }
        rank++;
    }

    vector<T> ret;
    for(int i=rank;i<n;i++)if(abs(A[i][m])>EPS) return ret;
    ret.assign(m,0);
    for(int i=0;i<rank;i++) ret[i]=A[i][m];
    return ret;
}
#line 2 "Math/PolynomialRestore_double.hpp"

// f(x) = a_0 + a_1 * x .. .a_d * x^d
// O(d ^ 3)
template<typename T>
vector<T> PolynomialRestore(vector<T> x,vector<T> y,int d,T EPS=1e-6){
    vector<vector<T>> M(x.size(),vector<T>(d+1));
    for(int i=0;i<(int)x.size();i++){
        T p=1;
        for(int j=0;j<=d;j++){
            M[i][j]=p;
            p*=x[i];
        }
    }
    auto ret=LinearEquation(M,y,EPS);
    return ret;
}

Back to top page