C言語でべき乗と累乗を高速に計算する方法を教えて！
こういった悩みにお答えします．
本記事の信頼性
- リアルタイムシステムの研究歴12年．
- 東大教員の時に，英語でOS（Linuxカーネル）の授業．
- 2012年9月～2013年8月にアメリカのノースカロライナ大学チャペルヒル校（UNC）コンピュータサイエンス学部で客員研究員として勤務．C言語でリアルタイムLinuxの研究開発．
- プログラミング歴15年以上，習得している言語: C/C++，Python，Solidity/Vyper，Java，Ruby，Go，Rust，D，HTML/CSS/JS/PHP，MATLAB，Assembler (x64，ARM)．
- 東大教員の時に，C++言語で開発した「LLVMコンパイラの拡張」，C言語で開発した独自のリアルタイムOS「Mcube Kernel」をGitHubにオープンソースとして公開．
- 2020年1月～現在はアメリカのノースカロライナ州チャペルヒルにあるGuarantee Happiness LLCのCTOとしてECサイト開発やWeb/SNSマーケティングの業務．2022年6月～現在はアメリカのノースカロライナ州チャペルヒルにあるJapanese Tar Heel, Inc.のCEO兼CTO．
- 最近は自然言語処理AIとイーサリアムに関する有益な情報発信に従事．
- （AI全般を含む）自然言語処理AIの論文の日本語訳や，AIチャットボット（ChatGPT，Auto-GPT，Gemini（旧Bard）など）の記事を50本以上執筆．アメリカのサンフランシスコ（広義のシリコンバレー）の会社でプロンプトエンジニア・マネージャー・Quality Assurance（QA）の業務委託の経験あり．
- （スマートコントラクトのプログラミングを含む）イーサリアムや仮想通貨全般の記事を200本以上執筆．イギリスのロンドンの会社で仮想通貨の英語の記事を日本語に翻訳する業務委託の経験あり．
こういった私から学べます．
C言語を独学で習得することは難しいです．
私にC言語の無料相談をしたいあなたは，公式LINE「ChishiroのC言語」の友だち追加をお願い致します．
私のキャパシティもあり，一定数に達したら終了しますので，今すぐ追加しましょう！
独学が難しいあなたは，元東大教員がおすすめするC言語を学べるオンラインプログラミングスクール5社で自分に合うスクールを見つけましょう．後悔はさせません！
べき乗と累乗を高速に計算
べき乗と累乗を高速に計算する方法を紹介します．
ここで，\(a^n\)におけるべき乗と累乗の違いは以下になります．（本記事ではaとnを自然数と仮定します．つまり，aが自然数の場合の累乗です．）
- べき乗：\(a^n\)（aとnは任意の数）
- 累乗：\(a^n\)（aは任意の数，nは自然数）
べき乗と累乗の計算でaを繰り返し乗算する場合，上記のaのn乗ではn - 1回の乗算が必要になります．
つまり，計算量は \(\mathcal{O}(n)\)になります．
べき乗と累乗を高速に計算する方法では，乗算の回数を減らすことを目標とします．
その分，加算，ビット演算やシフト演算の回数は増えますが，一般的には乗算の方が処理時間が長いので，全体として高速に計算できます．
べき乗と累乗をpow関数で実装するコードを知りたいあなたはこちらからどうぞ．
べき乗と累乗を高速に計算する方法は，仮想通貨の技術である暗号理論でよく利用されます．
仮想通貨のイーサリアムを知りたいあなたはこちらからどうぞ．
本記事では，コードを読みやすくするためにunsigned long型を採用していますが，べき乗と累乗の計算では算術オーバーフローが発生しやすいです．
C言語の任意精度算術ライブラリGMPを利用することで，算術オーバーフローを回避できます．
算術オーバーフローと回避方法を知りたいあなたはこちらからどうぞ．
繰り返し2乗法（バイナリ法）
繰り返し2乗法（バイナリ法：binary method）とは，べき乗を高速に計算するアルゴリズムです．
\((a^{x})^2 = a^{2x}\)となる性質を利用します．
つまり，\(a^1\)，\(a^2\)，\(a^4\)，...，の中で適切なものを選べば良いことがわかります．
例えば，\(a^9 = a^{8 + 1} = a^8 * a^1\)になります．
\(a^9\)をバイナリ法で計算すると以下のように4回の乗算になり，aを繰り返し8回乗算するより少ないことがわかります．
- \(a^2 = a^1 * a^1\)
- \(a^4 = a^2 * a^2\)
- \(a^8 = a^4 * a^4\)
- \(a^9 = a^8 * a^1\)
バイナリ法で\(a^n\)を求める計算量は\(\mathcal{O}(\log n)\)になります．
バイナリ法のコードは以下になります．
バイナリ法には，上位ビットから計算する右向きバイナリ法（left-to-right binary method）と下位ビットから計算する左向きバイナリ法（right-to-left binary method）があります．
右向きバイナリ法はleft2right_binary関数，左向きバイナリ法はright2left_binary関数になります．
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 |
/* * Author: Hiroyuki Chishiro * License: 2-Clause BSD */ #include <stdio.h> int get_msb(unsigned long n) { int i = 0; while (n > 0) { n >>= 1; i++; } /* return -1 if a == 0. */ return i - 1; } unsigned long left2right_binary(unsigned long a, unsigned long n) { unsigned long x = 1; int i; for (i = get_msb(n); i >= 0; i--) { x *= x; if ((n >> i) & 1) { x *= a; } } return x; } unsigned long right2left_binary(unsigned long a, unsigned long n) { unsigned long x = 1; int i = get_msb(n); int j = 0; for (j = 0; j <= i; j++) { if ((n >> j) & 1) { x *= a; } a *= a; } return x; } int main(void) { unsigned long a, n; do { printf("Please input two natural numbers: "); scanf("%lu%lu", &a, &n); } while (a == 0 || n == 0); printf("left2right_binary(%lu, %lu) = %lu\n", a, n, left2right_binary(a, n)); printf("right2left_binary(%lu, %lu) = %lu\n", a, n, right2left_binary(a, n)); return 0; } |
実行結果は以下になります．
3行目の行末で「2 5」を入力したら，\(2^5 = 32\)と正しく計算できていることがわかります．
|
1 2 3 4 5 |
$ gcc binary.c $ a.out Please input two natural numbers: 2 5 left2right_binary(2, 5) = 32 right2left_binary(2, 5) = 32 |
バイナリ法はサイドチャネル攻撃に弱いという欠点があります．
この理由として，コード\(n_i\)（nのi番目のビット）が1の時だけ発生する処理（28，45行目）があるので，終了までの処理時間が変わるからです．
nの値によって出力までの処理時間が変わるということは，サイドチャネル攻撃のタイミング攻撃によってnの値が露呈する恐れがあります．
モンゴメリべき乗法
モンゴメリべき乗法（Montgomery's ladder）は，バイナリ法とは異なり，サイドチャネル攻撃に強いアルゴリズムです．
この理由として，モンゴメリべき乗法では，nのi番目のビット\(n_i\)が1か0かで同等の計算を実行するからです．
モンゴメリべき乗法のコードは以下になります．
25～31行目で，ビットが1か0に関わらず同等の計算を行っていることがわかります．
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 |
/* * Author: Hiroyuki Chishiro * License: 2-Clause BSD */ #include <stdio.h> int get_msb(unsigned long n) { int i = 0; while (n > 0) { n >>= 1; i++; } /* return -1 if a == 0. */ return i - 1; } unsigned long montogomery(unsigned long a, unsigned long n) { unsigned long x = a, xx = a * a; int i; for (i = get_msb(n) - 1; i >= 0; i--) { if ((n >> i) & 1) { x *= xx; xx *= xx; } else { xx *= x; x *= x; } } return x; } int main(void) { unsigned long a, n; do { printf("Please input two natural numbers: "); scanf("%lu%lu", &a, &n); } while (a == 0 || n == 0); printf("montogomery(%lu, %lu) = %lu\n", a, n, montogomery(a, n)); return 0; } |
実行結果は以下になります．
|
1 2 3 4 |
$ gcc montgomery.c $ a.out Please input two natural numbers: 2 5 montogomery(2, 5) = 32 |
上記のモンゴメリべき乗法のコードは，サイドチャネル攻撃の中のキャッシュによるタイミング攻撃に対応していません．
キャッシュによりタイミング攻撃に対応したモンゴメリべき乗法の擬似コードはこちらにあります．
\(2^k\)-ary法
\(2^k\)-ary法は，\(a^2\)，\(a^3\)，\(a^4\)，...，\(a^{{2^d} - 1}\)（dは自然数）の計算結果を配列（テーブル）に保持して利用することで，乗算の回数を減らすアルゴリズムです．
d = 3の場合の\(2^k\)-ary法のコードは以下になります．
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 |
/* * Author: Hiroyuki Chishiro * License: 2-Clause BSD */ #include <stdio.h> #define D 3 #define K (1 << D) /* K = 2^D */ #define MASK (K - 1) /* MASK = K - 1 */ int get_msb(unsigned long n) { int i = 0; while (n > 0) { n >>= 1; i++; } /* return -1 if a == 0. */ return i - 1; } unsigned long _2k_ary(unsigned long a, unsigned long n) { unsigned long x = 1; int i, j; unsigned long table[K]; unsigned long remaining; table[0] = 1; table[1] = a; for (i = 2; i < K; i++) { table[i] = a * table[i - 1]; } i = get_msb(n); remaining = i % D; j = (n >> (i - remaining)) & MASK; x = table[j]; i -= remaining; while (i >= D) { for (j = 0; j < D; j++) { x *= x; } j = (n >> (i - D)) & MASK; x *= table[j]; i -= D; } return x; } int main(void) { unsigned long a, n; do { printf("Please input two natural numbers: "); scanf("%lu%lu", &a, &n); } while (a == 0 || n == 0); printf("_2k_ary(%lu, %lu) = %lu\n", a, n, _2k_ary(a, n)); return 0; } |
実行結果は以下になります．
|
1 2 3 4 |
$ gcc 2k_ary.c $ a.out Please input two natural numbers: 2 5 left2right_2k_ary(2, 5) = 32 |
スライディングウィンドウ法
スライディングウィンドウ法（Sliding Window Method）は，\(2^k\)-ary法を拡張したもので，最長dビットの1または0が連続するビット列をウィンドウと定義して，まとめて計算するアルゴリズムです．
\(2^k\)-ary法とは異なり，スライディングウィンドウ法は\(a^3\)，\(a^5\)，\(a^7\)，...，\(a^{{2^d} - 1}\)（dは自然数）の計算結果を配列（テーブル）に保持して利用します．
この理由は，1が連続するウィンドウの最下位ビットは必ず1（つまり奇数）になるためです．
0が連続するウィンドウの処理は，\(a = a^{2^d}\)を計算するのみです．（1が連続するウィンドウでも同様の処理を行います．）
d=3の場合のスライディングウィンドウ法を利用するコードは以下になります．
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 |
/* * Author: Hiroyuki Chishiro * License: 2-Clause BSD */ #include <stdio.h> #include <limits.h> #define D 3 struct windows { struct window { int one; int length; } window[sizeof(unsigned long) * CHAR_BIT]; int length; }; int get_msb(unsigned long n) { int i = 0; while (n > 0) { n >>= 1; i++; } /* return -1 if a == 0. */ return i; } void decompose_bits(unsigned long n, struct windows *windows) { int one, j; int i = get_msb(n) - 1; windows->length = 0; while (i >= 0) { one = (n >> i) & 1; for (j = 1; j < D && i >= j && one == ((n >> (i - j)) & 1); j++) { } windows->window[windows->length].one = one; windows->window[windows->length].length = j; windows->length++; i -= j; } } unsigned long sliding_window(unsigned long a, unsigned long n) { unsigned long x = 1, aa = a * a; int i, j; unsigned long table[D + 1]; struct windows windows; table[0] = a; for (i = 1; i < D + 1; i++) { table[i] = aa * table[i - 1]; } decompose_bits(n, &windows); for (i = 0; i < windows.length; i++) { for (j = 0; j < windows.window[i].length; j++) { x *= x; } if (windows.window[i].one) { j = (1 << (windows.window[i].length - 1)) - 1; x *= table[j]; } } return x; } int main(void) { unsigned long a, n; do { printf("Please input two natural numbers: "); scanf("%lu%lu", &a, &n); } while (a == 0 || n == 0); printf("sliding_window(%lu, %lu) = %lu\n", a, n, sliding_window(a, n)); return 0; } |
実行結果は以下になります．
|
1 2 3 4 |
$ gcc sliding_window.c $ a.out Please input two natural numbers: 2 5 sliding_window(2, 5) = 32 |
まとめ
C言語でべき乗と累乗を高速に計算する方法を紹介しました．
具体的には，以下のアルゴリズムのコードを解説しました．
- 繰り返し2乗法（バイナリ法）
- モンゴメリべき乗法
- \(2^k\)-ary法
- スライディングウィンドウ法
べき乗と累乗における乗算の回数を減らせるので，確かめてみましょう．
C言語を独学で習得することは難しいです．
私にC言語の無料相談をしたいあなたは，公式LINE「ChishiroのC言語」の友だち追加をお願い致します．
私のキャパシティもあり，一定数に達したら終了しますので，今すぐ追加しましょう！
独学が難しいあなたは，元東大教員がおすすめするC言語を学べるオンラインプログラミングスクール5社で自分に合うスクールを見つけましょう．後悔はさせません！
におけるべき乗と累乗の違いは以下になります．){kind=link}