ぽか

これはなに

ARC182-C の tester 解にある変形をわたしがどうやって思いつきましたか記事です（？）

259-momone.hatenablog.com

考えた道筋記事でもあるかもしれません（？？？）

6 \( (=\pi(M))\) 変数の fps（解説とだいたい同じです）

\(f=1+x _ 2+x _ 3+{x _ 2} ^ 2+x _ 5+x _ 2x _ 3+\dotsb{}\) みたいなのを考えると、\(\displaystyle\dfrac{f-f ^ {N+1}}{1-f}=\sum _ {x=1} ^ \infty\Biggl\lbrace\)長さが \(N\) 以下で積が \(x\) である列の個数 \(\displaystyle\times\prod _ {p\in\mathbb P} {x _ p} ^ {v _ p(x)}\Biggr\rbrace\) とできます。

ここで、\(\displaystyle\prod _ p {x _ p} ^ {v _ p(x)}\) を \(\displaystyle\prod _ p(v _ p(x)+1)\) に変換できれば問題を解くことができます。

この変換は

\(\displaystyle\prod _ p x _ p\) をかけて
すべての \(x _ p\) で \(1\) 回ずつ微分して
すべての \(x _ p\) に \(1\) を代入する

とできます。

微分するパートが難しそうです。ここで二重数（双対数）を用いた自動微分について想いを馳せました。

二重数

ざっくり言うと、\(\mathbb R\) に対する \(\mathbb R[\varepsilon]/(\varepsilon ^ 2)\) です（でどころが実数でない場合についても考えることができます）。

なにがうれしいかというと、\(\mathbb R\) 係数多項式 \(P\) があったとき、これに \(\mathbb R[\varepsilon]/(\varepsilon ^ 2)\) の元 \(a+\varepsilon\) を放り込むと \(P(a)+P ^ {\prime}(a)\varepsilon\) が得られることです（なんでそうなるかは各自で確かめてください）。

\(\sin,\cos,\log\) などの初等関数に対しても二重数に関するオーバーロードを定義しておくことで、これらの組み合わせで作られた関数の微分係数を自動的に（数式処理をすることなく）得ることができて、便利です。

今回の問題では、\(\varepsilon _ 2,\varepsilon _ 3,\varepsilon _ 5,\varepsilon _ 7,\varepsilon _ {11},\varepsilon _ {13}\) を \(\mathbb F _ {998244353}\) に入れることで微分するパートを処理することができます。

微分して \(1\) を代入した結果の値が得たかったものなので、\(x _ p\) に \(1+\varepsilon _ p\) を代入したときの \(\displaystyle\prod _ p\varepsilon _ p\) の係数が求める値です。これでもとの解説に戻ってきました。

のこり

解説中では \(\displaystyle\sum f ^ n\) の計算でわちゃわちゃ式変形がされていますが、単にダブリングしてしまったほうが式変形は少ない気がします（人それぞれ）。わたしは \(\dfrac{1-f ^ {N+1}}{1-f}-1\) をダブリングしました（初項がきれいになるので（？））。

また、\(\pi(M)\) 種の二重数の積が \(O(4 ^ {\pi(M)})\) で行われていますが、\(\displaystyle\prod _ {p\in S}\varepsilon _ p\) と \(\displaystyle\prod _ {p\in T}\varepsilon _ p\) をかけるのは \(S\cap T=\emptyset\) のときに限っていいので \(O(3 ^ {\pi(M)})\) とできます。集合べき級数を使うことでもうすこし速くできるらしいです？

コンテスト中の実装はこんな感じでした

#include <bits/extc++.h>
#include <atcoder/modint>

namespace atcoder {
template <int m>
std::ostream &operator<<(std::ostream &os, const atcoder::static_modint<m> &mint) {
    os << mint.val();
    return os;
}
} // namespace atcoder

int main() {
    using namespace std;
    using modint = atcoder::static_modint<998244353>;
    using diff_type = array<modint, 64>;
    unsigned long N;
    cin >> N;
    unsigned M;
    cin >> M;
    diff_type base{};
    for(unsigned i{1}; i <= M; ++i){
        auto x{i};
        array<unsigned, 6> factor_count{};
        for(auto& [j, p] : vector<pair<unsigned, unsigned>>{{0, 2}, {1, 3}, {2, 5}, {3, 7}, {4, 11}, {5, 13}}){
            while(x % p == 0){
                ++factor_count[j];
                x /= p;
            }
        }
        for(unsigned j{}; j < 64; ++j){
            modint tmp{1};
            for(unsigned b{}; b < 6; ++b)
                if(1 & (j >> b))
                    tmp *= factor_count[b];
            base[j] += tmp;
        }
    }
    
    const auto mul{[](const diff_type& lhs, const diff_type& rhs){
        diff_type result{};
        for(unsigned i{}; i < 64; ++i)
            for(unsigned j{i}; j < 64; ++j |= i)
                result[j] += lhs[i] * rhs[i ^ j];
        return result;
    }};

    ++N;
    // prefix + ∑ coef × base^i
    diff_type prefix{}, coef{};
    ranges::fill(coef, modint::raw(1));
    while(N){
        if(N & 1){
            for(unsigned i{}; i < 64; ++i)
                prefix[i] += coef[i];
            coef = mul(coef, base);
        }
        const auto tmp{mul(coef, base)};
        for(unsigned i{}; i < 64; ++i)
            coef[i] += tmp[i];
        base = mul(base, base);
        N /= 2;
    }
    cout << prefix.back() - 1 << endl;
    return 0;
}