algorithm - a^b^c^... mod m を見つける

Question

私は計算したい：

a ^{b ^{c ^{d ^{. ^{. ^.}}}}}mod m

この数値は大きすぎますが、 a 、 b 、 c 、 ... および m は単純な 32 ビット int に収まるため、効率的な方法をご存知ですか。

何か案は？

警告:この質問は、 ^b mod mを見つけることとは異なります。

^{また、 a b ^c}^{は (a b} ) ^cと同じではないことに注意してください。^{後者は a bc}と同じです。べき乗は右結合です。

score 24 · Accepted Answer

a ^{b ^c} mod m = a ^{b ^c mod n} mod m、ここで n = φ(m)オイラーの全関数.

m が素数の場合、n = m-1 です。

編集: Nabb が指摘したように、これは a が m と互いに素である場合にのみ有効です。したがって、最初にこれを確認する必要があります。

score 19 · Accepted Answer

答えには、正しさの完全な正式な数学的証明は含まれていません。ここでは不要と判断しました。その上、SOでは非常に読みにくいでしょう（たとえば、MathJaxではありません）。
(ほんの少しだけ) 特定の素因数分解アルゴリズムを使用します。それは最良の選択肢ではありませんが、十分です。

tl;dr

計算したいa^x mod m。function を使用しますmodpow(a,x,m)。以下で説明します。

xが十分に小さい (指数形式でない、または存在しない)場合はp^x | m、それを計算して返します。
素数に分割し、関数 p^x mod mを使用して素数ごとに個別に計算しますmodpow
1. 計算c' = gcd(p^x,m)してt' = totient(m/c')
2. 計算するw = modpow(x.base, x.exponent, t') + t'
3. テーブルpow(p, w - log_p c', m) * c'に保存A
A のすべての要素を乗算し、モジュロ m を返します

ここでpowは、python の pow のように見えるはずです。

主な問題：

現在のベストアンサーは特殊なケースのみについてgcd(a,m) = 1であり、OPはこの仮定を問題として考慮していなかったため、この回答を書くことにしました。オイラーのトティエント定理も使います。ウィキペディアを引用：

オイラーのトーティエント定理:とが互いに素な正の整数の場合
、 φ(n) はオイラーのトーティエント関数.na

numbers are co-primeNabbがコメントで示しているように、この仮定は非常に重要です。したがって、まず、数値が互いに素であることを確認する必要があります。（より明確にするために、と仮定しx = b^(c^...)ます。）なぜなら a^x mod m = ((p1^alpha)^x mod m)*(p2... 、 a = p1^アルファ * p2^... を因数分解aし、別々に計算q1 = (p1^alpha)^x mod m,q2 = (p2^beta)^x mod m...してから、簡単な方法で答えを計算できるから(q1 * q2 * q3 * ... mod m)です。数には多くてもo(log a)素因数があるため、ほとんどのo(log a)計算を実行する必要があります。

実際、すべての素因数に分割する必要はaなく (すべてがm他の指数で発生しない場合)、同じ指数と組み合わせることができますが、今のところ注目に値するものではありません。

ここで、素数(p^z)^x mod mがどこにあるかを見てみましょう。pいくつかの重要な観察に注意してください。

a,bがよりも小さい正の整数でありm、cが何らかの正の整数である場合、 true は文です。

上記の観察を使用して、実際の問題の解決策を得ることができます。簡単に計算できますgcd((p^z)^x, m)。mx*z が大きい場合、で割り切れる回数pです。しましょうm' = m /gcd((p^z)^x, m)。(注意(p^z)^x = p^(z*x)) しましょうc = gcd(p^(zx),m)。w = p^(zx - c) mod m'この数値は互いに素であるため、オイラーの定理を使用して簡単に計算できます (以下を参照) 。その後、上記の観察を使用して、を受け取ることができますp^(zx) mod m。上記の仮定からwc mod m'c = p^(zx) mod m、今のところ答えはでp^(zx) mod m = wcあり、w,cは簡単に計算できます。

したがって、簡単に計算できますa^x mod m。

`a^x mod m`オイラーの定理を使って計算する

ここa,mで、互いに素であると仮定します。計算したい場合はa^x mod m、計算t = totient(m)して通知できますa^x mod m = a^(x mod t) mod m。が大きく、例えばのようにxの特定の表現しか知らない場合に役立ちます。xx = 7^200

例を見てくださいx = b^c。アルゴリズムを時間で二乗することにより、べき乗を計算t = totient(m)しx' = b^c mod tて使用できます。そして（同じアルゴリズムを使用して）後、これは解と同じです。Θ(log c)a^x' mod m

x = b^(c^(d^...)再帰的に解決する場合。最初t1 = totient(m)に、その後t2 = totient(t1)などを計算します。たとえば、x=b^(c^d)。もし、、、t1=totient(m)そして、どこで、a^x mod m = a^(b^(c^d) mod t1)と言うことができます。二乗アルゴリズムによる累乗を使用して計算しているすべてのもの。注: 一部の totient が指数と互いに素でない場合は、主な問題と同じトリックを使用する必要があります (実際には、それが指数であることを忘れて、主な問題のように再帰的に問題を解決する必要があります)。上記の例で、が c と互いに素でない場合、このトリックを使用する必要があります。b^(c^d) mod t1 = b^(c^d mod t2) mod t1t2 = totient(t1)t2

計算する`φ(n)`

単純な事実に注目してください:

もしgcd(a,b)=1、その後φ(ab) = φ(a)*φ(b)
p素数の場合φ(p^k)=(p-1)*p^(k-1)

nしたがって、 (ak. )を因数分解し、事実 2 を使用しn = p1^k1 * p2^k2 * ...て個別に計算できますφ(p1^k1),φ(p2^k2),...。次に、事実 1 を使用してこれを結合します。φ(n)=φ(p1^k1)*φ(p2^k2)*...

totient を繰り返し計算する場合は、エラトステネスのふるいを使用して、表に素数を保存することをお勧めします。定数を減らします。

Pythonの例：（この因数分解アルゴリズムと同じ理由で正しい）

def totient(n) :          # n - unsigned int
    result = 1
    p = 2                 #prime numbers - 'iterator'
    while p**2 <= n :
        if(n%p == 0) :    # * (p-1)
            result *= (p-1)
            n /= p
        while(n%p == 0) : # * p^(k-1)
            result *=  p
            n /= p
        p += 1
    if n != 1 :
        result *= (n-1)
    return result         # in O(sqrt(n))

場合：`a`^{b^c}`mod m`

実際には同じことを何度も行っているため、このケースはこれを一般的に解決する方法を示していると思います。
まず、a主要な勢力に分割する必要があります。最良の表現はペアになります<number, exponent>。
C++11 の例:

std::vector<std::tuple<unsigned, unsigned>> split(unsigned n) {
  std::vector<std::tuple<unsigned, unsigned>> result;
  for(unsigned p = 2; p*p <= n; ++p) {
    unsigned current = 0;
    while(n % p == 0) {
      current += 1;
      n /= p;
     }
    if(current != 0)
     result.emplace_back(p, current);
   }
  if(n != 1)
   result.emplace_back(n, 1);
  return result;
 }

分割後、(p^z)^(b^c) mod m=p^(z*(b^c)) mod mペアごとに計算する必要があります。まず、かどうかを確認する必要がありp^(z*(b^c)) | mます。z,b,cはいの場合、答えは (p^z)^(b^c) ですが、非常に小さい場合にのみ可能です。コード例を表示する必要はないと思います。
最後にp^(z*b^c) > m、答えを計算する必要がある場合。まず、計算する必要がありますc' = gcd(m, p^(z*b^c))。計算できるようになったらt = totient(m')。と (z*b^c - c' mod t)。答えを得るための簡単な方法です。

function modpow(p, z, b, c, m : integers) # (p^z)^(b^c) mod m
    c' = 0
    m' = m
    while m' % p == 0 :
        c' += 1
        m' /= p
    # now m' = m / gcd((p^z)^(b^c), m)
    t = totient(m')
    exponent = z*(b^c)-c' mod t
    return p^c' * (p^exponent mod m')

そして以下のPythonの作業例：

def modpow(p, z, b, c, m) : # (p^z)^(b^c) mod m
    cp = 0
    while m % p == 0 :
        cp += 1
        m /= p              # m = m' now
    t = totient(m)
    exponent = ((pow(b,c,t)*z)%t + t - (cp%t))%t
                            # exponent = z*(b^c)-cp mod t
    return pow(p, cp)*pow(p, exponent, m)

この関数を使用すると、を簡単に計算できます。(p^z)^(b^c) mod mすべての結果 ( mod m) を乗算するだけで、継続的にすべてを計算することもできます。以下の例。(書き間違えていないことを祈ります。) 仮定、b、c だけが十分に大きく ( b^c > log(m)ak. それぞれは割り切れp^(z*b^k)ませんm)、単純なチェックであり、それによってごちゃごちゃにするポイントがわかりません。

def solve(a,b,c,m) : # split and solve
    result = 1
    p = 2            # primes
    while p**2 <= a :
        z = 0
        while a % p == 0 :
                     # calculate z
            a /= p
            z += 1
        if z != 0 :
            result *=  modpow(p,z,b,c,m)
            result %= m
        p += 1
    if a != 1 :      # Possible last prime
        result *= modpow(a, 1, b, c, m)
    return result % m

それが動作するように見えます。
デモとそれは正しいです!

score 1 · Accepted Answer

どの関係a=x^yでも、関係は使用している基数 (基数 2、基数 6、基数 16 など) に関して不変です。
mod N 演算は、基数 N の最下位桁 (LSD) を抽出することと同等であるため、
基数 N の結果 A の LSD は、基数 N の X の LSD によってのみ影響を受ける可能性があり、上位の桁の数字には影響されないためです。(例: 34*56 = 30*50+30*6+50*4+4*5 = 10*(3+50+3*6+5*4)+4*6)

したがって、次のことからLSD(A)=LSD(X^Y)推測できます。

LSD(A)=LSD(LSD(X)^Y)

したがって

A mod N = ((X mod N) ^ Y) mod N

と

(X ^ Y) mod N = ((X mod N) ^ Y) mod N)

したがって、各パワーステップの前に mod を実行できます。これにより、結果が整数の範囲内に保たれます。

これは、a が負ではなく、任意の x^y に対して a^y < MAXINT であると仮定します。

この答えは間違った質問に答えます。(アレックス)

score 0 · Accepted Answer

A^X mod M増加するにつれての動作を見てくださいX。最終的にはサイクルに入る必要があります。サイクルに長さがあり、ステップのP後に開始するとしNます。次に、をX >= N意味しA^X = A^(X+P) = A^(X%P + (-N)%P + N) (mod M)ます。したがって、を計算することで計算できA^B^Cますy=B^C, z = y < N ? y : y%P + (-N)%P + N, return A^z (mod m)。

導出された方程式には、指数 <Mまたはより小さい被除数を持つより小さい指数タワーを含む指数があるため、この戦略を再帰的に累乗ツリーに適用できることに注意してください。

N唯一の問題は、P与えられたAとを効率的に計算できるかどうかですM。過大評価しNても問題ないことに注意してください。Nに設定するだけでM、うまくいきます。Pは少し難しいです。AとMが異なる素数の場合、P=M-1. AがすべてMの素因数を持つ場合、0 とで行き詰まりP=1ます。方法がわからないので、それを理解するための演習として残します。

///Returns equivalent to list.reverse().aggregate(1, acc,item => item^acc) % M
func PowerTowerMod(Link<int> list, int M, int upperB = M)
    requires M > 0, upperB >= M
    var X = list.Item
    if list.Next == null: return X
    var P = GetPeriodSomehow(base: X, mod: M)
    var e = PowerTowerMod(list.Next, P, M)
    if e^X < upperB then return e^X //todo: rewrite e^X < upperB so it doesn't blowup for large x
    return ModPow(X, M + (e-M) % P, M)

score 0 · Accepted Answer

剰余累乗は、この問題を解決する正しい方法です。ここにちょっとしたヒントがあります: ^{a b}^{^c}^{^{^d}}

% m を見つけるには、a % m、次に^{a b} % m、次に a ^b^{^c} % mを計算することから始めなければなりません。 a ^b^{^c}^{^{^d}} % m ... (お分かりでしょう)

^{a b} % mを見つけるには、基本的に 2 つのアイデアが必要です: [B=floor(b/2) とする]

a ^b = (a ^B ) ² b が偶数の場合、または a ^b = (a ^B ) ² *b が奇数の場合。

(X*Y)%m = ((X%m) * (Y%m)) %m
(% = mod)

したがって、
b が偶数の場合 a
b ^% m = (a ^B % m) ² % m
または b が奇数の場合 a
b ^% m = (((a ^B % m) ² ) * (a % m)) % m

^{したがって、 B}の値がわかっている場合は、この値を計算できます。^B

を見つけるには、同様のアプローチを適用して、1 になるまで B を分割します。

^{例: 16 13} % 11を計算するには:

16 ¹³ % 11 = (16 % 11) ¹³ % 11 = 5 ¹³ % 11 = (5 ⁶ % 11) * (5 ⁶ % 11) * (5 % 11) <---- (I) 5 ⁶

を求めるには% 11: 5 ⁶ % 11 = ((5 ³ % 11) * (5 ³ % 11)) % 11 <----(II) 5 ³ %11: 5 ³ % 11 = ((5 ¹ % 11) * (5 ¹ % 11) * (5 % 11)) % 11 = (((5 * 5) % 11) * 5) % 11 = ((25 % 11) * 5) % 11 = (3 * 5) % 11 = 15 % 11 = 4 この値を (II) に代入すると、 5 ^{6が得られます。}

% 11 = (((4 * 4) % 11) * 5) % 11 = ((16 % 11) * 5) % 11 = (5 * 5) % 11 = 25 % 11 = 3
この値を (I ) を与える
5 ¹³ % 11 = ((3 % 11) * (3 % 11) * 5) % 11 = ((9 % 11) * 5) % 11 = 45 % 11 = 4

このように 5 ¹³ % 11 = 4これにより、 ⁵^¹³
% 11 などの形式のものを計算できます...

score 0 · Accepted Answer

タセットの答えは良いですが、大幅な単純化が可能です。

x のべき乗、mod m は前周期的です。x が m と互いに素である場合、x の累乗は周期的ですが、その仮定がなくても、周期の前の部分は長くなく、多くても m の素因数分解の指数の最大値であり、多くても log_2 m です。 . 周期の長さは phi(m) を分割し、実際には lambda(m) を分割します。ここで、lambda はCarmichael の関数であり、最大乗法次数 mod m です。これは、phi(m) よりも大幅に小さくすることができます。Lambda(m) は、phi(m) と同様に、m の素因数分解からすばやく計算できます。Lambda(m) は、m の素因数分解におけるすべての素数ベキ p_i^e_i に対する lambda(p_i^e_i) の GCD であり、素数ベキが奇数の場合、lambda(p_i^e_i) = phi(p_i^e_i) です。lambda(2)=1, lambda(4)=2, lambda(2^n)=2^(n-2) より大きな 2 のべき乗の場合。

modPos(a,n) を {0,1,..,n-1} の a の合同クラスの代表となるように定義します。非負の a の場合、これは単なる a%n です。負の場合、何らかの理由で a%n が負として定義されているため、modPos(a,n) は (a%n)+n です。

modMin(a,n,min) を、最小である mod n に一致する最小の正の整数となるように定義します。正の場合、これを min+modPos(a-min,n) として計算できます。

b^c^... が log_2 m より小さい場合 (この不等式が成り立つかどうかは、対数を再帰的に取ることで確認できます)、単純に a^b^c^... を計算できます。それ以外の場合は a^b^c^ ... mod m = a^modMin(b^c^..., lambda(m), [log_2 m])) mod m = a^modMin(b^c^... mod lambda(m), lambda (m)、[log_2 m])。

たとえば、2^3^4^5 mod 100 を計算したいとします。3^4^5 は 489 桁しかないため、他の方法でも実行できますが、計算したくないほど大きいことに注意してください。それを直接。ただし、ここで示した方法を使用すると、2^3^4^5 mod 100 を手動で計算できます。

3^4^5 > log_2 100 なので、

2^3^4^5 mod 100 
= 2^modMin(3^4^5,lambda(100),6) mod 100 
= 2^modMin(3^4^5 mod lambda(100), lambda(100),6) mod 100
= 2^modMin(3^4^5 mod 20, 20,6) mod 100.

3^4^5 mod 20 を計算してみましょう。4^5 > log_2 20 なので、

3^4^5 mod 20 
= 3^modMin(4^5,lambda(20),4) mod 20 
= 3^modMin(4^5 mod lambda(20),lambda(20),4) mod 20
= 3^modMin(4^5 mod 4, 4, 4) mod 20
= 3^modMin(0,4,4) mod 20
= 3^4 mod 20
= 81 mod 20
= 1

これを前の計算に組み込むことができます。

2^3^4^5 mod 100 
= 2^modMin(3^4^5 mod 20, 20,6) mod 100
= 2^modMin(1,20,6) mod 100
= 2^21 mod 100
= 2097152 mod 100
= 52.

2^(3^4^5 mod 20) mod 100 = 2^1 mod 100 = 2 であることに注意してください。これは正しくありません。ベースの力の前周期部分まで減らすことはできません。

algorithm - a^b^c^... mod m を見つける

a b c d . . . mod m

6 に答える 6

tl;dr

主な問題：

a^x mod mオイラーの定理を使って計算する

計算するφ(n)

場合：abc mod m

Related

Reference

a ^{b ^{c ^{d ^{. ^{. ^.}}}}}mod m

`a^x mod m`オイラーの定理を使って計算する

計算する`φ(n)`

場合：`a`^{b^c}`mod m`