newton interpolation

#include <vector>
#include <iostream>

using namespace std;

void print(vector<double> s)
{
    
    for(double x:s)
    {
    
        cout<<x<<" ";
    }
    cout<<endl;
}

vector<double> mul(vector<double> a,vector<double> b)    // Polynomial multiplication 
{
    
    int n=a.size(), m=b.size();
    vector<double> s(n+m-1,0);
    for(int i=0;i<n;++i)
    {
    
        for(int j=0;j<m;++j)
        {
    
            s[i+j]+=a[i]*b[j];
        }
    }
    return s;
}

vector<double> add(vector<double> a,vector<double> b)  
{
    
    int n=a.size(), m=b.size();
    vector<double> s(max(n,m),0);
    for(int i=0;i<n;++i){
    
        s[i]+=a[i];
    }
    for(int i=0;i<m;++i){
    
        s[i]+=b[i];
    }
    return s;
}

void print_chashang_table(vector<vector<double>> A)
{
    
    int n=A.size(), m=A[0].size();
    for(int i=0;i<n;++i)
    {
    
        for(int j=0;j<=i;++j)
        {
    
            printf("% 5.2f ",A[i][j]);
        }
        printf("\n");
    }
}

vector<double> Newton(vector<vector<double>> Point){
        
    int n=Point.size()-1;
    if(n==-1){
    
        cout<<" error "<<endl;
    }
    vector<vector<double>> diffCoef_Table;
    for(int i=0;i<n+1;++i){
               //n+1 Item difference quotient , Initialize the difference quotient table 
        vector<double> V(i+1,0);
        diffCoef_Table.emplace_back(V);
        diffCoef_Table[i][0]=Point[i][1];    // The first column of the difference quotient table is the function value 
    }
    for(int k=1;k<=n;++k){
            // Recursively calculate the difference quotient table 
        for(int i=k;i<=n;++i){
    
            diffCoef_Table[i][k]=(diffCoef_Table[i][k-1]-diffCoef_Table[i-1][k-1])/(Point[i][0]-Point[i-k][0]);
        }
    }
    print_chashang_table(diffCoef_Table);
    vector<double> ans(n+1,0);
    ans[0]=Point[0][1];
    for(int k=1;k<=n;++k){
    
        vector<double>s={
    diffCoef_Table[k][k]};
        for(int j=0;j<=k-1;++j){
    
            s=mul(s,{
    -Point[j][0],1});
        }
        ans=add(ans,s);
    }
    return ans;
}
int main(){
    
    vector<vector<double>> P={
    {
    0,0},{
    2,0},{
    1,1},{
    5,3}};
    vector<double> s=Newton(P);
    cout<<endl;
    print(s);
    system("pause");
    return 0;
}