Lagrange interpolation
(number/lagrange_interpolation.hpp)

• View this file on GitHub
• Last update: 2021-08-30 04:35:37+09:00
• Include: #include "number/lagrange_interpolation.hpp"

Depends on

• utils/macros.hpp

Code

#include <cassert>
#include <utility>
#include <vector>
#include "../utils/macros.hpp"

/**
 * @brief Lagrange interpolation
 * @note O(n^2)
 * @note works on any field
 * @note the error is very big
 */
std::vector<long double> lagrange_interpolate(const std::vector<std::pair<long double, long double> > & points) {
    // given ((x_0, y_1), \dots, (x_{n - 1}, y_{n - 1}))
    int n = points.size();

    // functions
    auto mul_monomial = [&](std::vector<long double> & f, long double a) {  // f \gets f \cdot (x - a)
        assert (f[n] == 0);
        REP_R (i, n + 1) {
            f[i] = (i >= 1 ? f[i - 1] : 0) - a * f[i];
        }
    };
    auto div_monomial = [&](std::vector<long double> & f, long double a) {  // f \gets f / (x - a)
        long double last = 0;
        REP_R (i, n + 1) {
            std::swap(f[i], last);
            last += a * f[i];
        }
    };
    auto apply = [&](std::vector<long double> & f, long double x) {  // f(x)
        long double y = 0;
        REP_R (i, n + 1) {
            y = y * x + f[i];
        }
        return y;
    };

    // let g = \prod_i (x - x_i)
    // let g_i = g / (x - x_i)
    // assert for all j \ne i. g_i(x_j) = 0
    std::vector<long double> g(n + 1);
    g[0] = 1;
    for (auto point : points) {
        mul_monomial(g, point.first);
    }

    // let c_i = y_i / g_i(x_i)
    // let f = \sum_i c_i g_i / (x - x_i)
    // assert for all i. f(x_i) = y_i
    std::vector<long double> f(n + 1);
    for (auto point : points) {
        div_monomial(g, point.first);
        long double c_i = point.second / apply(g, point.first);
        REP (j, n + 1) {
            f[j] += c_i * g[j];
        }
        mul_monomial(g, point.first);
    }
    assert (f.back() == 0);
    f.pop_back();
    return f;
}

#line 1 "number/lagrange_interpolation.hpp"
#include <cassert>
#include <utility>
#include <vector>
#line 2 "utils/macros.hpp"
#define REP(i, n) for (int i = 0; (i) < (int)(n); ++ (i))
#define REP3(i, m, n) for (int i = (m); (i) < (int)(n); ++ (i))
#define REP_R(i, n) for (int i = (int)(n) - 1; (i) >= 0; -- (i))
#define REP3R(i, m, n) for (int i = (int)(n) - 1; (i) >= (int)(m); -- (i))
#define ALL(x) std::begin(x), std::end(x)
#line 5 "number/lagrange_interpolation.hpp"

/**
 * @brief Lagrange interpolation
 * @note O(n^2)
 * @note works on any field
 * @note the error is very big
 */
std::vector<long double> lagrange_interpolate(const std::vector<std::pair<long double, long double> > & points) {
    // given ((x_0, y_1), \dots, (x_{n - 1}, y_{n - 1}))
    int n = points.size();

    // functions
    auto mul_monomial = [&](std::vector<long double> & f, long double a) {  // f \gets f \cdot (x - a)
        assert (f[n] == 0);
        REP_R (i, n + 1) {
            f[i] = (i >= 1 ? f[i - 1] : 0) - a * f[i];
        }
    };
    auto div_monomial = [&](std::vector<long double> & f, long double a) {  // f \gets f / (x - a)
        long double last = 0;
        REP_R (i, n + 1) {
            std::swap(f[i], last);
            last += a * f[i];
        }
    };
    auto apply = [&](std::vector<long double> & f, long double x) {  // f(x)
        long double y = 0;
        REP_R (i, n + 1) {
            y = y * x + f[i];
        }
        return y;
    };

    // let g = \prod_i (x - x_i)
    // let g_i = g / (x - x_i)
    // assert for all j \ne i. g_i(x_j) = 0
    std::vector<long double> g(n + 1);
    g[0] = 1;
    for (auto point : points) {
        mul_monomial(g, point.first);
    }

    // let c_i = y_i / g_i(x_i)
    // let f = \sum_i c_i g_i / (x - x_i)
    // assert for all i. f(x_i) = y_i
    std::vector<long double> f(n + 1);
    for (auto point : points) {
        div_monomial(g, point.first);
        long double c_i = point.second / apply(g, point.first);
        REP (j, n + 1) {
            f[j] += c_i * g[j];
        }
        mul_monomial(g, point.first);
    }
    assert (f.back() == 0);
    f.pop_back();
    return f;
}

Back to top page