$\sum_{k} r^k\cdot poly(k)$
(FPS/sumofpolyexp.hpp)

View this file on GitHub
Last update: 2024-01-14 04:04:38+09:00
Include: #include "FPS/sumofpolyexp.hpp"

Depends on

interpolate (one point) (FPS/interpolate.hpp)

Verified with

Code

#pragma once
#include "FPS/interpolate.hpp"

template <typename T>
T LimitSumOfPolyExp(vector<T> &f, T r) { // sum_{k=0}^inf r^k*f(k)
    assert(r != 1);
    int d = f.size() - 1;
    vector<T> rs(d + 1);
    rs[0] = 1;
    rep(i, 0, d) rs[i + 1] = rs[i] * r;
    T c, add;
    rep(i, 0, d + 1) {
        add += rs[i] * f[i];
        if ((d - i) & 1)
            c -= nCr<T>(d + 1, i + 1) * rs[d - i] * add;
        else
            c += nCr<T>(d + 1, i + 1) * rs[d - i] * add;
    }
    c /= (-r + 1).pow(d + 1);
    return c;
}

template <typename T>
T SumOfPolyExp(vector<T> &f, T r, ll n) { // sum_{k=0}^{n-1} r^k*f(k)
    n--;
    if (n < 0)
        return 0;
    int d = f.size() - 1;
    vector<T> rs(d + 1), rui(d + 1);
    rs[0] = 1;
    rep(i, 0, d) rs[i + 1] = rs[i] * r;
    rep(i, 0, d + 1) rui[i] = rs[i] * f[i];
    rep(i, 0, d) rui[i + 1] += rui[i];
    if (r == 0)
        return f[0];
    else if (r == 1)
        return Interpolate(rui, n);
    else {
        T c;
        rep(i, 0, d + 1) c +=
            nCr<T>(d + 1, i + 1) * rs[d - i] * rui[i] * ((d - i) & 1 ? -1 : 1);
        c /= T(-r + 1).pow(d + 1);
        vector<T> ys(d + 1);
        T pwr = 1, invr = T(r).inv();
        rep(i, 0, d + 1) ys[i] = (rui[i] - c) * pwr, pwr *= invr;
        return T(r).pow(n) * Interpolate(ys, n) + c;
    }
}

/**
 * @brief $\sum_{k} r^k\cdot poly(k)$
 */

#line 2 "FPS/interpolate.hpp"

template<typename T>T Interpolate(vector<T>& ys,ll t){ // f(0),..,f(d) -> f(t)
    int d=ys.size()-1;
    if(t<=d)return ys[t];
    vector<T> L(d+1,1),R(d+1,1);
    rep(i,0,d)L[i+1]=L[i]*(t-i);
    for(int i=d;i;i--)R[i-1]=R[i]*(t-i);
    T ret;
    rep(i,0,d+1){
        T add=ys[i]*L[i]*R[i]*Fact<T>(i,1)*Fact<T>(d-i,1);
        if((d-i)&1)ret-=add;
        else ret+=add;
    }
    return ret;
}

/**
 * @brief interpolate (one point)
*/
#line 3 "FPS/sumofpolyexp.hpp"

template <typename T>
T LimitSumOfPolyExp(vector<T> &f, T r) { // sum_{k=0}^inf r^k*f(k)
    assert(r != 1);
    int d = f.size() - 1;
    vector<T> rs(d + 1);
    rs[0] = 1;
    rep(i, 0, d) rs[i + 1] = rs[i] * r;
    T c, add;
    rep(i, 0, d + 1) {
        add += rs[i] * f[i];
        if ((d - i) & 1)
            c -= nCr<T>(d + 1, i + 1) * rs[d - i] * add;
        else
            c += nCr<T>(d + 1, i + 1) * rs[d - i] * add;
    }
    c /= (-r + 1).pow(d + 1);
    return c;
}

template <typename T>
T SumOfPolyExp(vector<T> &f, T r, ll n) { // sum_{k=0}^{n-1} r^k*f(k)
    n--;
    if (n < 0)
        return 0;
    int d = f.size() - 1;
    vector<T> rs(d + 1), rui(d + 1);
    rs[0] = 1;
    rep(i, 0, d) rs[i + 1] = rs[i] * r;
    rep(i, 0, d + 1) rui[i] = rs[i] * f[i];
    rep(i, 0, d) rui[i + 1] += rui[i];
    if (r == 0)
        return f[0];
    else if (r == 1)
        return Interpolate(rui, n);
    else {
        T c;
        rep(i, 0, d + 1) c +=
            nCr<T>(d + 1, i + 1) * rs[d - i] * rui[i] * ((d - i) & 1 ? -1 : 1);
        c /= T(-r + 1).pow(d + 1);
        vector<T> ys(d + 1);
        T pwr = 1, invr = T(r).inv();
        rep(i, 0, d + 1) ys[i] = (rui[i] - c) * pwr, pwr *= invr;
        return T(r).pow(n) * Interpolate(ys, n) + c;
    }
}

/**
 * @brief $\sum_{k} r^k\cdot poly(k)$
 */

$\sum_{k} r^k\cdot poly(k)$ (FPS/sumofpolyexp.hpp)

Depends on

Verified with

Code

$\sum_{k} r^k\cdot poly(k)$
(FPS/sumofpolyexp.hpp)