Math Class

$$F(x):=OGF(f):=\sum_{i=0}^\infty f_ix^i, \\x^{\lfloor i\rfloor}:=\prod_{j=0}^{i-1}(x-j),$$ 如果 $$\\OGF(f)=\sum_{i=0}^\infty a_ix^{\lfloor i\rfloor}$$ 那么 $$\exp(-x)\sum_{i=0}^\infty\frac{F(i)}{i!}x^i=\sum_{i=0}^\infty a_ix^i.$$
证明： $$\exp(-x)\sum_{i=0}^\infty\frac{\sum_{j=0}^i a_j\frac{i!}{(i-j)!}}{i!}x^i=\sum_{k=0}^\infty\frac{1}{k!}(-x)^k\sum_{i=0}^\infty\sum_{j=0}^ia_j\frac1{(i-j)!}x^i \\=\sum_{n=0}^\infty\sum_{i=0}^n\sum_{j=0}^i\frac1{(i-j)!(n-i)!}a_j(-1)^{n-i}x^n \\=\sum_{j=0}^\infty a_j\sum_{n=j}^\infty x^n(-1)^n\sum_{i=j}^n\frac{1}{(i-j)!(n-i)!}(-1)^i \\=\sum_{j=0}^\infty a_j\sum_{n=j}^\infty x^n(-1)^n\sum_{k=0}^{n-j}\frac{1}{k!(n-j-k)!}(-1)^{k+j} \\=\sum_{j=0}^\infty a_jx^j$$

实际只需要$a_{0\sim n}$，所以直接$\exp(-x)\sum_{i=0}^n\frac{F(i)}{i!}x^i$就可。
显然$a_{n+1\sim p-1}=0$所以直接$\exp(x)\sum_{i=0}^na_ix^i$就得出了$F(0\sim p-1)$

有了$F(1\sim p-1)$，接着就是bluestein算法，算法核心是

$$F(w^j)=\sum_{i=0}^{p-2}f_iw^{ij}=\sum_{i=0}^{p-2}f_iw^{C_{i+j}^2-C_i^2-C_j^2}=w^{-C_j^2}\sum_{i=0}^{p-2}f_iw^{-C_i^2}\cdot w^{C_{i+j}^2},\\ f_j=\frac{1}{p-1}\sum_{i=0}^{p-2}F(w^i)w^{-ij}=\frac1{p-1}\sum_{i=0}^{p-2}F(w^i)w^{-C_{i+j}^2+C_i^2+C_j^2}=\frac{1}{p-1}w^{C_j^2}\sum_{i=0}^{p-2}F(w^i)w^{C_i^2}\cdot w^{-C_{i+j}^2}$$

#include <bits/stdc++.h>
using ul = std::uint32_t;
using li = std::int32_t;
using ll = std::int64_t;
using ull = std::uint64_t;
using ulll = __uint128_t;
using llf = long double;
using lf = double;
using vul = std::vector<ul>;
using vvul = std::vector<vul>;
using pulb = std::pair<ul, bool>;
using vpulb = std::vector<pulb>;
using vvpulb = std::vector<vpulb>;
using vb = std::vector<bool>;
using vull = std::vector<ull>;
const ull base = 4179340454199820289ull;
ull plus(ull a, ull b)
{
    return a + b < base ? a + b : a + b - base;
}
ull minus(ull a, ull b)
{
    return a < b ? a + base - b : a - b;
}
ull mul(ull a, ull b)
{
    ll tmp = ll(a) * ll(b) - ll(llf(a) * llf(b) / llf(base)) * ll(base);
    return tmp < 0 ? tmp + ll(base) : (tmp < base ? tmp : tmp - ll(base));
}
void exgcd(ll a, ll b, ll& x, ll& y)
{
    if (b) {
        exgcd(b, a % b, y, x);
        y -= x * (a / b);
    } else {
        x = 1;
        y = 0;
    }
}
ull inverse(ull a)
{
    ll x, y;
    exgcd(a, base, x, y);
    return x < 0 ? x + ll(base) : x;
}
ull pow(ull a, ull b)
{
    ull ret = 1;
    ull temp = a;
    while (b) {
        if (b & 1) {
            ret = mul(ret, temp);
        }
        temp = mul(temp, temp);
        b >>= 1;
    }
    return ret;
}
class polynomial : public vull {
public:
    void clearzero() {
        while (size() && !back()) {
            pop_back();
        }
    }
    polynomial()=default;
    polynomial(const vull& a): vull(a) { }
    polynomial(vull&& a): vull(std::move(a)) { }
    ul degree() const {
        return size() - 1;
    }
};
class ntt_t {
public:
    static const ul lgsz = 21;
    static const ul sz = 1 << lgsz;
    static const ull g = 3;
    ull w[sz + 1];
    ull leninv[lgsz + 1];
    ntt_t() {
        ull wn = pow(g, (base - 1) >> lgsz);
        w[0] = 1;
        for (ul i = 1; i <= sz; ++i) {
            w[i] = mul(w[i - 1], wn);
        }
        leninv[lgsz] = inverse(sz);
        for (ul i = lgsz - 1; ~i; --i) {
            leninv[i] = plus(leninv[i + 1], leninv[i + 1]);
        }
    }
    void operator()(vull& v, ul& n, bool inv) {
        ul lgn = 0;
        while ((1 << lgn) < n) {
            ++lgn;
        }
        n = 1 << lgn;
        v.resize(n, 0);
        for (ul i = 0, j = 0; i != n; ++i) {
            if (i < j) {
                std::swap(v[i], v[j]);
            }
            ul k = n >> 1;
            while (k & j) {
                j &= ~k;
                k >>= 1;
            }
            j |= k;
        }
        for (ul lgmid = 0; (1 << lgmid) != n; ++lgmid) {
            ul mid = 1 << lgmid;
            ul len = mid << 1;
            for (ul i = 0; i != n; i += len) {
                for (ul j = 0; j != mid; ++j) {
                    ull t0 = v[i + j];
                    ull t1 = mul(w[inv ? (len - j << lgsz - lgmid - 1) : (j << lgsz - lgmid - 1)], v[i + j + mid]);
                    v[i + j] = plus(t0, t1);
                    v[i + j + mid] = minus(t0, t1);
                }
            }
        }
        if (inv) {
            for (ul i = 0; i != n; ++i) {
                v[i] = mul(v[i], leninv[lgn]);
            }
        }
    }
} ntt;
polynomial& operator*=(polynomial& a, const polynomial& b)
{
    if (!b.size() || !a.size()) {
        a.resize(0);
        return a;
    }
    polynomial temp = b;
    ul npmp1 = a.size() + b.size() - 1;
    if (ull(a.size()) * ull(b.size()) <= ull(npmp1) * ull(50)) {
        temp.resize(0);
        temp.resize(npmp1, 0);
        for (ul i = 0; i != a.size(); ++i) {
            for (ul j = 0; j != b.size(); ++j) {
                temp[i + j] = plus(temp[i + j], mul(a[i], b[j]));
            }
        }
        a = temp;
        a.clearzero();
        return a;
    }
    ntt(a, npmp1, false);
    ntt(temp, npmp1, false);
    for (ul i = 0; i != npmp1; ++i) {
        a[i] = mul(a[i], temp[i]);
    }
    ntt(a, npmp1, true);
    a.clearzero();
    return a;
}
polynomial operator*(const polynomial& a, const polynomial& b)
{
    polynomial ret = a;
    return ret *= b;
}
const ul maxp = 5e5;
ul T;
ul n, p;
ul F[maxp];
ul f[maxp];
ul fac[maxp];
ul fiv[maxp];
polynomial tmp1;
polynomial tmp2;
std::vector<ul> facpm1;
ul wpow[maxp - 1];
void exgcd(li a, li b, li& x, li& y)
{
    if (b) {
        exgcd(b, a % b, y, x);
        y -= a / b * x;
    } else {
        x = 1;
        y = 0;
    }
}
ul pow(ul a, ul b, ul mod)
{
    ul ret = 1;
    while (b) {
        if (b & 1) {
            ret = ull(ret) * ull(a) % mod;
        }
        a = ull(a) * ull(a) % mod;
        b >>= 1;
    }
    return ret;
}
ull c2(ul x)
{
    return ull(x) * ull(x - 1) / 2;
}
int main()
{
    std::scanf("%u", &T);
    for (ul Case = 1; Case <= T; ++Case) {
        std::scanf("%u%u", &n, &p);
        fac[0] = 1;
        for (ul i = 1; i <= p - 1; ++i) {
            fac[i] = ull(fac[i - 1]) * ull(i) % p;
        }
        if (true) {
            li x, y;
            exgcd(fac[p - 1], p, x, y);
            if (x < 0) {
                x += p;
            }
            fiv[p - 1] = x;
        }
        for (ul i = p - 1; i >= 1; --i) {
            fiv[i - 1] = ull(fiv[i]) * ull(i) % p;
        }
        for (ul i = 0; i <= n; ++i) {
            std::scanf("%u", &F[i]);
        }
        tmp1.resize(n + 1);
        tmp2.resize(n + 1);
        for (ul i = 0; i <= n; ++i) {
            tmp1[i] = ull(F[i]) * ull(fiv[i]) % p;
            tmp2[i] = (i & 1) ? p - fiv[i] : fiv[i];
        }
        tmp1 *= tmp2;
        tmp1.resize(n + 1);
        for (ul i = 0; i <= n; ++i) {
            tmp1[i] %= p;
        }
        tmp1.resize(p);
        tmp2.resize(p);
        for (ul i = 0; i <= p - 1; ++i) {
            tmp2[i] = fiv[i];
        }
        tmp1 *= tmp2;
        tmp1.resize(p);
        for (ul i = 0; i <= p - 1; ++i) {
            F[i] = ull(tmp1[i] % p) * ull(fac[i]) % p;
        }
        ul pm1 = p - 1;
        facpm1.resize(0);
        for (ul i = 2; pm1 != 1; ++i) {
            while (pm1 % i == 0) {
                pm1 /= i;
                facpm1.push_back(i);
            }
        }
        ul w;
        pm1 = p - 1;
        for (w = 1; ; ++w) {
            bool flag = true;
            for (ul i : facpm1) {
                if (pow(w, pm1 / i, p) == 1) {
                    flag = false;
                    break;
                }
            }
            if (flag) {
                break;
            }
        }
        wpow[0] = 1;
        for (ul i = 1; i <= p - 2; ++i) {
            wpow[i] = ull(wpow[i - 1]) * ull(w) % p;
        }
        tmp1.resize(p - 1);
        for (ul i = 0; i <= p - 2; ++i) {
            tmp1[i] = ull(F[wpow[i]]) * ull(wpow[c2(i) % pm1]) % p;
        }
        tmp2.resize(p - 1 + p - 1);
        for (ul i = 0; i <= p + p - 3; ++i) {
            tmp2[p + p - 3 - i] = wpow[(pm1 - c2(i) % pm1) % pm1];
        }
        tmp1 *= tmp2;
        tmp1.resize(p - 1 + p - 1);
        for (ul i = 0; i <= n; ++i) {
            ul pos = p + p - 3 - i;
            f[i] = (p - ull(tmp1[pos] % p) * ull(wpow[c2(i) % pm1]) % p) % p;
        }
        std::printf("Case #%u:\n", Case);
        for (ul i = 0; i <= n; ++i) {
            if (i) {
                std::putchar(' ');
            }
            std::printf("%u", f[i]);
        }
        std::printf("\n");
    }
    return 0;
}