Product Of Differences
(library/math/product_of_differences.hpp)

• View this file on GitHub
• Last update: 2023-05-11 13:37:15+09:00
• Include: #include "library/math/product_of_differences.hpp"

Product Of Differences

Depends on

• Multi Point Evaluation (library/polynomial/multi_point_eval.hpp)

Required by

• $\displaystyle \sum _ i i ^ d r ^ i$ (library/math/sum_i^d_r^i.hpp)
• ラグランジュ補間 (library/polynomial/lagrange_interpolation.hpp)

Verified with

• test/src/math/product_of_differences/yuki1938.test.cpp
• test/src/math/sum_i^d_r^i/sum_of_exponential_times_polynomial.test.cpp
• test/src/math/sum_i^d_r^i/sum_of_exponential_times_polynomial_limit.test.cpp
• test/src/polynomial/lagrange_interpolation/cumulative_sum.test.cpp
• test/src/polynomial/lagrange_interpolation/dummy.test.cpp
• test/src/polynomial/lagrange_interpolation/dummy_2.test.cpp

Code

#ifndef SUISEN_PRODUCT_OF_DIFFERNCES
#define SUISEN_PRODUCT_OF_DIFFERNCES

#include <deque>
#include "library/polynomial/multi_point_eval.hpp"

namespace suisen {
    /**
     * O(N(logN)^2)
     * return the vector p of length xs.size() s.t. p[i]=Π[j!=i](x[i]-x[j])
     */
    template <typename FPSType, typename T>
    std::vector<typename FPSType::value_type> product_of_differences(const std::vector<T>& xs) {
        // f(x):=Π_i(x-x[i])
        // => f'(x)=Σ_i Π[j!=i](x-x[j])
        // => f'(x[i])=Π[j!=i](x[i]-x[j])
        const int n = xs.size();
        std::deque<FPSType> dq;
        for (int i = 0; i < n; ++i) dq.push_back(FPSType{ -xs[i], 1 });
        while (dq.size() >= 2) {
            auto f = std::move(dq.front());
            dq.pop_front();
            auto g = std::move(dq.front());
            dq.pop_front();
            dq.push_back(f * g);
        }
        auto f = std::move(dq.front());
        f.diff_inplace();
        return multi_point_eval<FPSType, T>(f, xs);
    }
} // namespace suisen


#endif // SUISEN_PRODUCT_OF_DIFFERNCES

#line 1 "library/math/product_of_differences.hpp"



#include <deque>
#line 1 "library/polynomial/multi_point_eval.hpp"



#include <vector>

namespace suisen {
    template <typename FPSType, typename T>
    std::vector<typename FPSType::value_type> multi_point_eval(const FPSType& f, const std::vector<T>& xs) {
        int n = xs.size();
        if (n == 0) return {};
        std::vector<FPSType> seg(2 * n);
        for (int i = 0; i < n; ++i) seg[n + i] = FPSType{ -xs[i], 1 };
        for (int i = n - 1; i > 0; --i) seg[i] = seg[i * 2] * seg[i * 2 + 1];
        seg[1] = f % seg[1];
        for (int i = 2; i < 2 * n; ++i) seg[i] = seg[i / 2] % seg[i];
        std::vector<typename FPSType::value_type> ys(n);
        for (int i = 0; i < n; ++i) ys[i] = seg[n + i].size() ? seg[n + i][0] : 0;
        return ys;
    }
} // namespace suisen


#line 6 "library/math/product_of_differences.hpp"

namespace suisen {
    /**
     * O(N(logN)^2)
     * return the vector p of length xs.size() s.t. p[i]=Π[j!=i](x[i]-x[j])
     */
    template <typename FPSType, typename T>
    std::vector<typename FPSType::value_type> product_of_differences(const std::vector<T>& xs) {
        // f(x):=Π_i(x-x[i])
        // => f'(x)=Σ_i Π[j!=i](x-x[j])
        // => f'(x[i])=Π[j!=i](x[i]-x[j])
        const int n = xs.size();
        std::deque<FPSType> dq;
        for (int i = 0; i < n; ++i) dq.push_back(FPSType{ -xs[i], 1 });
        while (dq.size() >= 2) {
            auto f = std::move(dq.front());
            dq.pop_front();
            auto g = std::move(dq.front());
            dq.pop_front();
            dq.push_back(f * g);
        }
        auto f = std::move(dq.front());
        f.diff_inplace();
        return multi_point_eval<FPSType, T>(f, xs);
    }
} // namespace suisen

Back to top page