P1998 题解 / DGF 广义求逆

飞雨烟雁 · 2025-03-27 19:34:40 · 题解

前置知识是 DGF ln / exp，推荐阅读 DGF 的计算理论：牛顿迭代与特殊求逆运算。

下设数论函数 f 对应的 DGF 为 F（g,h 同理），并记 F 的第 n 项为 F_n。

给定 DGF F 和正整数 n,m，我们应该如何快速计算 F+F^2+\cdots+F^m 的前 n 项呢？

熟悉等比数列求和的读者会说，F+\cdots+F^m=(F^{m+1}-1)/(F-1)-1，然后采用 LOJ #6713 的方法就行了。不过我们会发现，如果 F_1=1 的话，F-1 就不可逆，最后会卡在求逆上。而且，因为 DGF 本质是多元多项式，我们不能像一元多项式一样分子分母除个 x 再求逆。那么，我们有没有什么不用倍增的快速做法呢？

有的兄弟，有的。这就是广义求逆。

适当泛化问题，给定 DGF F,G，若已知存在 H 使得 F=GH，且 F_1=G_1=0,G_2\neq 0，那要如何求出 H？

根据 Dirichlet 卷积的式子：

稍微整理一下就是：

依此式递推即可，时间复杂度为 \Theta(n\log n)。

这个做法需要得知 F,G 的前 2n 项，才能算出 H 的前 n 项。看似浪费了很多已知信息，实则我们已经利用了所有信息。因为想要唯一确定 h(n)，就必须知道 f(2n) 是多少。即便是递推式中未出现的 f 的奇数值也很重要，它可以告诉我们这样的 H 是否存在。不过，如果我们已知 H 是存在的，那奇数值确实没啥用（

代码如下：

#include <ctime>
#include <cstdio>
#include <iostream>
#define ll long long
using namespace std;

int Read(){
    int res = 0; char c = getchar();
    while(c < '0' || c > '9') c = getchar();
    while(c >= '0' && c <= '9') res = res * 10 + (c ^ 48), c = getchar();
    return res;
}

const int Mx = 2e6 + 10, Mod = 1e9 + 7;

inline int FastPow(ll a, int b){
    int res = 1;
    while(b){
        if(b & 1) res = res * a % Mod;
        b >>= 1, a = a * a % Mod;
    }
    return res;
}
ll Inverse(ll x){ return FastPow(x % Mod, Mod - 2);}

void Print(const int *F, int n){
    int ans = 0;
    for(int i = 1; i <= n; ++i) ans ^= ((F[i] < 0) ? (F[i] + Mod) : F[i]);
    cout << ans;
}

bool Vis[Mx];
int Prime[Mx], tot;
int Lns[Mx], Inv[30];
void Sieve(){
    for(int i = 2; i < Mx; ++i){
        if(!Vis[i]) Prime[++tot] = i, Lns[i] = 1;
        for(int j = 1; j <= tot && Prime[j] * i < Mx; ++j){
            Vis[i * Prime[j]] = true;
            Lns[i * Prime[j]] = 1 + Lns[i];
            if(i % Prime[j] == 0) break;
        }
    }
    for(int i = 1; i < 30; ++i) Inv[i] = Inverse(i);
}

int Temp[Mx], Temp2[Mx], Temp3[Mx], Temp4[Mx];
void Mul(const int *F, const int *G, int *H, int n){
    for(int i = 1; i <= n; ++i) Temp[i] = 1ll * F[i] * G[1] % Mod;
    for(int i = 1; i <= n; ++i) for(int j = (i << 1); j <= n; j += i)
        Temp[j] = (Temp[j] + 1ll * F[i] * G[j / i]) % Mod;
    for(int i = 1; i <= n; ++i) H[i] = Temp[i];
}

void Div(const int *F, const int *G, int *H, int n){
    for(int i = 1; i <= n; ++i) Temp[i] = F[i];
    ll Invg = Inverse(G[1]);
    for(int i = 1; i <= n; ++i){
        Temp[i] = Temp[i] * Invg % Mod;
        for(int j = (i << 1); j <= n; j += i) Temp[j] = (Temp[j] - 1ll * Temp[i] * G[j / i]) % Mod;
    }
    for(int i = 1; i <= n; ++i) H[i] = Temp[i];
}

void Derivate(const int *F, int *G, int n){ for(int i = 1; i <= n; ++i) G[i] = 1ll * F[i] * Lns[i] % Mod;}
void Integrate(const int *F, int *G, int n){ for(int i = 1; i <= n; ++i) G[i] = 1ll * F[i] * Inv[Lns[i]] % Mod;}

void Logarithm(const int *F, int *G, int n){ Derivate(F, Temp2, n), Div(Temp2, F, Temp3, n), Integrate(Temp3, G, n);}

void Exponent(const int *F, int *G, int n){
    Derivate(F, Temp2, n);
    for(int i = 2; i <= n; i++) Temp3[i] = 0;
    Temp3[1] = 1;
    for(int i = 1; i <= n; i++){
        if(i != 1) Temp3[i] = 1ll * Temp3[i] * Inv[Lns[i]] % Mod;
        for(int j = (i << 1); j <= n; j += i) Temp3[j] = (Temp3[j] + 1ll * Temp3[i] * Temp2[j / i]) % Mod;
    }
    for(int i = 1; i <= n; ++i) G[i] = Temp3[i];
}

void Pow(const int *F, int *G, int n, int k){
    Logarithm(F, Temp4, n);
    for(int i = 1; i <= n; ++i) Temp4[i] = 1ll * Temp4[i] * k % Mod;
    Exponent(Temp4, G, n);
}

int F[Mx], G[Mx], H[Mx];
int main(){

    Sieve();

    int n = Read(), m = Read();

    for(int i = 1; i <= n; ++i) F[i] = Read();

    Pow(F, G, 2 * n, m + 1);
    G[1] = F[1] = 0;

    for(int i = 1; i <= n; ++i) H[i] = G[2 * i];
    ll ivs = Inverse(F[2]);
    for(int t = 1; t <= n; ++t){
        H[t] = H[t] * ivs % Mod;
        for(int d = 3; d * t <= 2 * n; ++d) if(d * t % 2 == 0){
            H[d * t / 2] = (H[d * t / 2] - 1ll * F[d] * H[t]) % Mod;
        }
    }

    --H[1];
    Print(H, n);

    return 0;
}