この記事は、アドベントカレンダー2021 24 日目の記事です。

mod p でそのあと指数に乗るかもしれなくても大丈夫な modint だれか作ってたりしませんか(O(log p)個の整数を持ちます)

— ぽかーん懐古DP@259家(桃音モモ) (@259_Momone) September 30, 2021

ふむふむ、今日はこれをつくりましょう。

要件定義

• 非負整数 P を指定します。
• Modint<P> は非負整数を 1 つ持ちます。
• 持っている非負整数を P で割った余りを取得できます。
• Modint<P> の Modint<P> 乗ができます。
• 0 の 0 乗は 1 です。

🆗

オイラーの定理

$\text{mod }P$ での値を持つだけでは pow は実現できません。例えば、$P=10^9+7$ のとき、$2^0=1$ ですが、$2^P=2$ と、指数の部分は $\text{mod }P$ で扱えないことがわかります。実際の値を持つわけにもいきませんから、何か周期性を見つけて管理したいです。

オイラーの定理

$\phi (P)=$( $1$ 以上 $P$ 以下の整数のうち $P$ と互いに素なものの個数 ) とする。
$P$ と互いに素な整数 $a$ について、

$$a^{\phi (P)}=1(\mathrm{mod}P)$$

が成り立つ。

つまり、指数の部分は $\text{mod }\phi (P)$ で管理すれば良いということですね。(？)

設計

というわけで、以下のように設計すればよさそうです。

• $\text{mod }P$ での値を持ちます。
• $\text{mod }P$ での pow ができるように $\text{mod }\phi (P)$ での値も持ちます。
• $\text{mod }\phi (P)$ での pow ができるように $\text{mod }\phi (\phi (P))$ での値も持ちます。
• $\text{mod }\phi (\phi (P))$ での pow ができるように $\text{mod }\phi (\phi (\phi (P)))$ での値も持ちます。
  $⋮$

$x>1$ において $\phi (x)<x$ なので有限回で終わりますが、実際何個の値を持つことになるでしょうか？

https://www.wolframalpha.com/input/?i=NestList[EulerPhi%2C+10^9+%2B+7%2C+28]

$P=10^9+7$ で $29$ 個、意外と少ないですね。

実は、これが $O(\log P)$ 個になることが証明できます。

定理

$P$ に $O(\log P)$ 回 $\phi$ を適用すると $1$ になる。

証明

$\phi (x)$ がどのような計算であるかに注目しましょう。
$\phi (x)$ は以下のように計算できます。

1. $x$ を素因数分解します。$270\to \{2,3,3,3,5\}$
2. 重複を取り除きます。$\{2,3,3,3,5\}\to \{2,3,5\}$
3. $x$ を $\frac{p-1}{p}$ 倍したやつが答えです。$270\times \frac{1}{2}\times \frac{2}{3}\times \frac{4}{5}=72$

$x$ が偶数であるとき、$\frac{1}{2}$ 倍が掛かります。
$x$ が奇数であるとき、$\frac{p-1}{p}$ 倍が掛かって偶数になります。

2 回に 1 回は半分になるので、$O(\log P)$ 回で $1$ になることが分かります。
(もっと厳しく、$\lceil {\log }_{2}P\rceil$ 回で $1$ になることが証明できます。)

問題点

しかし、これでは問題が残っています。

1 つは、$0$ の $0$ 乗です。

$P=3,\phi (P)=2$ のとき、$0^0=1,0^6=0$ ですが、$\text{mod }3$ と $\text{mod }2$ での値を持つだけでは $0$ と $6$ が区別できません。

もう 1 つは、互いに素でないときオイラーの定理が成り立たない問題です。

$P=10,\phi (P)=4$ のとき、$2^{\phi (P)}=6$ となって、オイラーの定理が成り立ちません。
$(1\to 2\to 4\to 8\to 6\to 2\to 4\to \cdots )$
これは、1 度 $2$ を掛けると $2$ の倍数になってしまい、$1$ に戻って来れないことにより発生します。しかし、1 度 $2$ の倍数になってしまえば、そこからは $\phi (P)$ 周期の周期性が成り立っています。

解決策

$0$ の $0$ 乗の問題は、指数部分が $0$ か $0$ でないかが分かれば解決します。
オイラーの定理が成り立たない問題は、指数部分が周期に入らない小さい数か、周期に入る大きい数かが分かれば解決します。
そこで、単に $\text{mod }P$ での値を管理するのではなく、$0$ から $P-1$ まではそのまま、$P$ 以上のものは $\text{mod }P$ して $P$ 以上 $2P-1$ 以下に収めるとすればいいでしょう。

実装

constexpr uint32_t totient(uint32_t x){
    uint32_t ans = x;
    for(uint32_t i = 2; i * i <= x; i++) if(x % i == 0){
        ans /= i;
        ans *= i - 1;
        do x /= i; while(x % i == 0);
    }
    if(x != 1){
        ans /= x;
        ans *= x - 1;
    }
    return ans;
}
template<uint32_t P> struct Modint{
    static_assert(P < 0x80000000, "P must be smaller than 2^31");
    uint32_t a;
    Modint<totient(P)> b;
private:
    static uint32_t mod(uint64_t x){
        if(x < P * 2) return uint32_t(x);
        return uint32_t(x % P) + P;
    }
    static uint32_t mul(uint32_t a, uint32_t b){
        return mod(uint64_t(a) * b);
    }
    static uint32_t pow(uint32_t a, uint32_t b){
        uint32_t ans = 1;
        while(b){
            if(b & 1) ans = mul(ans, a);
            a = mul(a, a);
            b >>= 1;
        }
        return ans;
    }
public:
    Modint(uint64_t x): a(mod(x)), b(x){}
    Modint(uint32_t a, Modint<totient(P)> b): a(a), b(b){}
    uint32_t val() const {
        if(a < P) return a;
        return a - P;
    }
    Modint pow(const Modint& other) const {
        return {pow(a, other.b.a), b.pow(other.b)};
    };
};
template<> struct Modint<1>{
    uint32_t a;
    Modint(uint64_t x): a(bool(x)){}
    uint32_t val() const {
        return 0;
    }
    Modint pow(const Modint& other) const {
        return {a || !other.a};
    }
};

pow しかできないの？

低機能すぎるので機能を増やしましょう。

• 掛け算 : できます。
• 足し算 : できます。
• 引き算 : できません。持っている値が 0 かどうか判定するのに実際の値を持つ必要が出てきますし、そもそも (負の数) 乗ができません。

constexpr uint32_t totient(uint32_t x){
    uint32_t ans = x;
    for(uint32_t i = 2; i * i <= x; i++) if(x % i == 0){
        ans /= i;
        ans *= i - 1;
        do x /= i; while(x % i == 0);
    }
    if(x != 1){
        ans /= x;
        ans *= x - 1;
    }
    return ans;
}
template<uint32_t P> struct Modint{
    static_assert(P < 0x80000000, "P must be smaller than 2^31");
    uint32_t a;
    Modint<totient(P)> b;
private:
    static uint32_t mod(uint64_t x){
        if(x < P * 2) return uint32_t(x);
        return uint32_t(x % P) + P;
    }
    static uint32_t mul(uint32_t a, uint32_t b){
        return mod(uint64_t(a) * b);
    }
    static uint32_t pow(uint32_t a, uint32_t b){
        uint32_t ans = 1;
        while(b){
            if(b & 1) ans = mul(ans, a);
            a = mul(a, a);
            b >>= 1;
        }
        return ans;
    }
public:
    Modint(uint64_t x): a(mod(x)), b(x){}
    Modint(uint32_t a, Modint<totient(P)> b): a(a), b(b){}
    uint32_t val() const {
        if(a < P) return a;
        return a - P;
    }
    Modint& operator*=(const Modint& other){
        a = mul(a, other.a);
        b *= other.b;
        return *this;
    }
    Modint operator*(const Modint& other) const {
        return Modint(*this) *= other;
    }
    Modint& operator+=(const Modint& other){
        a += other.a;
        if(a >= P * 2) a -= P;
        if(a >= P * 2) a -= P;
        b += other.b;
        return *this;
    }
    Modint operator+(const Modint& other) const {
        return Modint(*this) += other;
    }
    Modint pow(const Modint& other) const {
        return {pow(a, other.b.a), b.pow(other.b)};
    };
};
template<> struct Modint<1>{
    uint32_t a;
    Modint(uint64_t x): a(bool(x)){}
    uint32_t val() const {
        return 0;
    }
    Modint& operator*=(const Modint& other){
        a &= other.a;
        return *this;
    }
    Modint operator*(const Modint& other) const {
        return Modint(*this) *= other;
    }
    Modint& operator+=(const Modint& other){
        a |= other.a;
        return *this;
    }
    Modint operator+(const Modint& other) const {
        return Modint(*this) += other;
    }
    Modint pow(const Modint& other) const {
        return {a || !other.a};
    }
};

問題を解こう

ABC228 E - Integer Sequence Fair

$M^{K^N}\mathrm{mod}998244353$ を計算する問題ですね。

int main(){
    uint64_t n, k, m;
    cin >> n >> k >> m;
    Modint<998244353> N = n, K = k, M = m;
    cout << M.pow(K.pow(N)).val() << endl;
}

提出結果

Xmas Contest 2018 J - Japanese Exponentation

まさにこのためにあるような構造体ですね。
なかなか見た目のすごい問題ですが、この構造体があればあとは構文解析をするだけです。
日本語表記を BNF で形式的に表してみましょう。

<expression> ::= <億> | <expression> の <expression> 乗
<億> ::= <万> | <千> 億 <万>
<万> ::= <千> | <千> 万 <千>
<千> ::= <百> | <二> 千 <百>
<百> ::= <十> | <二> 百 <十>
<十> ::= <一> | <二> 十 <一>
<一> ::= [空文字列] | 〇 | … | 九
<二> ::= [空文字列] | 二 | … | 九

<expression> が空文字列を受理してしまいますが、この方が簡単なので気にしないことにしましょう。
<expression> を読むには、一番後ろの文字が「乗」であるかを見て分岐すれば良さそうです。
<億> を読むには、後ろから <万> を読み、その後、一番後ろの文字が「億」であるかを見て分岐すれば良さそうです。
他も同じようにできますね。
あとはこれに従って関数を書くだけです。
1 文字 3 byte なのがポイントですね。

using M = Modint<1000000009>;
string get_char(string_view s){
    if(s.empty()) return "";
    return {s.end() - 3, s.end()};
}
void next(string_view& s){
    s = s.substr(0, s.size() - 3);
}
const vector<string> digit = {"〇", "一", "二", "三", "四", "五", "六", "七", "八", "九"};
M 一(string_view& s){
    const auto p = find(digit.begin(), digit.end(), get_char(s));
    if(p == digit.end()) return 0;
    next(s);
    return p - digit.begin();
}
M 二(string_view& s){
    const auto p = find(digit.begin(), digit.end(), get_char(s));
    if(p == digit.end()) return 1;
    next(s);
    return p - digit.begin();
}
M 十(string_view& s){
    M a = 一(s);
    if(get_char(s) != "十") return a;
    next(s);
    return 二(s) * 10 + a;
}
M 百(string_view& s){
    M a = 十(s);
    if(get_char(s) != "百") return a;
    next(s);
    return 二(s) * 100 + a;
}
M 千(string_view& s){
    M a = 百(s);
    if(get_char(s) != "千") return a;
    next(s);
    return 二(s) * 1000 + a;
}
M 万(string_view& s){
    M a = 千(s);
    if(get_char(s) != "万") return a;
    next(s);
    return 千(s) * 10000 + a;
}
M 億(string_view& s){
    M a = 万(s);
    if(get_char(s) != "億") return a;
    next(s);
    return 千(s) * 1'0000'0000 + a;
}
M expression(string_view& s){
    if(get_char(s) != "乗") return 億(s);
    next(s);
    M b = expression(s);
    assert(get_char(s) == "の");
    next(s);
    M a = expression(s);
    return a.pow(b);
}
string encode(uint32_t x){
    function<string(uint32_t)> f = [&](uint32_t x){
        if(x >= 100000000) return f(x / 100000000) + "億" + f(x % 100000000);
        if(x >= 10000) return f(x / 10000) + "万" + f(x % 10000);
        if(x >= 2000) return f(x / 1000) + "千" + f(x % 1000);
        if(x >= 1000) return "千" + f(x % 1000);
        if(x >= 200) return f(x / 100) + "百" + f(x % 100);
        if(x >= 100) return "百" + f(x % 100);
        if(x >= 20) return f(x / 10) + "十" + f(x % 10);
        if(x >= 10) return "十" + f(x % 10);
        return vector<string>{"", "一", "二", "三", "四", "五", "六", "七", "八", "九"}[x];
    };
    string ans = f(x);
    if(ans.empty()) ans = "〇";
    return ans;
}
int main(){
    string s;
    cin >> s;
    string_view p = s;
    cout << encode(expression(p).val()) << endl;
}

提出結果

おわり

結構大変だったけど構文解析はおもしろいな〜

明日の担当は @Kashiwade さんです！楽しみ〜