Interpolation using Neville's Algorithm in C++

gistfile1.txt

class poly
{
    const std::vector<double> xs;
    const std::vector<double> ys;
    std::vector<double> data;

    auto& qs(int i, int j){return data[i * xs.size() + j];}
public:

    poly(const std::vector<double> xs, const std::vector<double> ys): xs(xs), ys(ys), data(xs.size() * xs.size())
    {
        assert(xs.size() == ys.size());
        //std::copy(ys.begin(), ys.end(), data.begin());
        for(int i = 0; i < xs.size(); i++)
        {
            qs(i, 0) = ys[i];
        }
    }

    double neville(size_t i, size_t j, double x);
    double operator()(double x);
};

double poly::neville(size_t i, size_t j, double x)
{
    for(int i = 1; i < std::ssize(xs); i++)
    {
        for(int j = i; j < std::ssize(xs); j++)
        {
            qs(j, i) = ((x - xs[j-i]) * qs(j,i-1) - (x - xs[j]) * qs(j-1,i-1)) / (xs[j] - xs[j-i]);
        }
    }
    return qs(xs.size() - 1, xs.size() - 1);
}

double poly::operator()(double x)
{
    return neville(0, xs.size() - 1, x);
}