feature image

2021年12月6日 | ブログ記事

指数に乗せられる Modint と Japanese Exponentation

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

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

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

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

要件定義

🆗

オイラーの定理

での値を持つだけでは pow は実現できません。例えば、 のとき、 ですが、 と、指数の部分は で扱えないことがわかります。実際の値を持つわけにもいきませんから、何か周期性を見つけて管理したいです。

オイラーの定理

( 以上 以下の整数のうち と互いに素なものの個数 ) とする。
互いに素な整数 について、

が成り立つ。

つまり、指数の部分は で管理すれば良いということですね。(?)

設計

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

において なので有限回で終わりますが、実際何個の値を持つことになるでしょうか?

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

個、意外と少ないですね。

実は、これが 個になることが証明できます。

定理

を適用すると になる。

証明

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

  1. を素因数分解します。
  2. 重複を取り除きます。
  3. 倍したやつが答えです。

が偶数であるとき、 倍が掛かります。
が奇数であるとき、 倍が掛かって偶数になります。

2 回に 1 回は半分になるので、 回で になることが分かります。
(もっと厳しく、 回で になることが証明できます。)

問題点

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

1 つは、 乗です。

のとき、 ですが、 での値を持つだけでは が区別できません。

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

のとき、 となって、オイラーの定理が成り立ちません。

これは、1 度 を掛けると の倍数になってしまい、 に戻って来れないことにより発生します。しかし、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 しかできないの?

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

実装
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

を計算する問題ですね。

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 さんです!楽しみ〜

tatyam icon
この記事を書いた人
tatyam

競技プログラミングと音ゲーと素数大富豪をしています。 AtCoder銅冠 (https://atcoder.jp/users/tatyam) ICPC 2020 Asia Yokohama Regional 3 位 (good_yamikin) ICPC 2021 Asia Yokohama Regional 3 位 (tonosama)

この記事をシェア

このエントリーをはてなブックマークに追加
共有

関連する記事

2021年12月8日
C++ with JUCEでステレオパンを作ってみた【AdC2021 26日目】
liquid1224 icon liquid1224
2021年4月18日
ベズー係数とN項の拡張ユークリッドの互除法
0214sh7 icon 0214sh7
2022年8月30日
【競プロer向け】母関数を習得しよう!
tatyam icon tatyam
2021年11月25日
【一人カラオケ活用術】安く済ませるボーカルレコーディング【AdC2021 13日目】
liquid1224 icon liquid1224
2021年11月24日
クリスマスデートをしよう!!!
mera icon mera
2021年4月26日
CPCTF2021 作問者writeup by hukuda222
hukuda222 icon hukuda222
記事一覧 タグ一覧 Google アナリティクスについて