RSA暗号の仕組みを追ってみた
最近では量子コンピュータの技術が発展してきて2030年頃には2048bitのRSA暗号が突破されると言われているらしいです。RSA暗号は馴染みが深いところでいうとhttpsでのssl・tls通信でも使われている技術になりますので改めて勉強してみようかと思います。 確認の実装にはRustを使っています。
RSAの歴史
公開鍵暗号の概念は1976年にデルフィ、ヘルマンにより発表されたがその時点では具体的な暗号の実現には至らなかった。翌年の1977年にリベスト、シャミア、エーデルマンが具体的な公開鍵暗号が提唱され3人の頭文字をとってRSA暗号と命名された。ところが英国政府通信本部(GCHQ)が1973年に同等の暗号を考案していたが守秘義務があったためそのことを公表することが出来なかったと言われています。英国政府通信本部は映画イミテーションゲームでも題材になったアラン・チューリングのナチス暗号機エニグマの解読の功績が有名です。
RSAのアルゴリズム
鍵の生成(公開鍵と秘密鍵)
鍵生成処理の流れ
まず素数pとqをランダムに生成する。ssh-keygenコマンドで鍵を生成するときに指定するbit数はここで決める素数p,qの値のものになります。それからN=pqを求め、GCD(φ(N), e) = 1となるeを選択します。GCDは最大公約数のことで1ということは φ(N)とeが互いに素であることを意味しています。φ(N)はオイラーの関数というものでφ(N)は1からNの範囲でNに対して素な数の合計となります。オイラーの関数についてm、nが互いに素であるときφ(mn)=φ(m)φ(n)が成り立ち、これを乗法性と言います。ここで使うφ(N)はφ(pq)でありp,qは互いに素なのでφ(N)=φ(p)φ(q)=(p-1)
(q-1)となります。(1~pの範囲で素数pに対して素な数はp-1になります。)よってGCD(φ(N), e) = 1の条件を満たすeを見つけることはGCD((p-1)(q-1), e) = 1を満たすeを見つけるのと同じになります。それからed = 1 mod φ(N)となるdを見つけます。modは剰余演算のことでed % φ(N) = ed % (p-1)(q-1) = 1と同様です。条件を満たすdを見つけるには拡張ユークリッドの互除法を使います。
鍵の生成はこれで完了で公開鍵はNとeで秘密鍵はdとなります。
簡潔にすると以下のようになる
1. 素数p,qをランダムに生成
2. N = pq を計算する
3. GCD(φ(N), e) = 1となるeを選択する
4. ed = 1 mod φ(N) を満たすdを計算する
5. 公開鍵(N, e)、秘密鍵(d)を返す
では、これらについて確認のため一つずつ実装してみたいと思います。
素数p,qをランダムに生成
素数の一覧を返す関数として以下を実装しました。
fn primes(a: i32) -> Vec<i32> {
let mut vec = Vec::new();
if a >= 2 { vec.push(2); }
if a >= 3 { vec.push(3); }
fn check(a: i32, vec: &Vec<i32>) -> bool {
let mut result = true;
for prime in vec {
if a % prime == 0 {
result = false;
break;
}
}
result
}
for i in 1..((a+1)/6 + 1) {
let num1 = i * 6 - 1;
if check(num1, &vec) { vec.push(num1); }
let num2 = i * 6 + 1;
if num2 <= a && check(num2, &vec) { vec.push(num2); }
}
vec
}
6の倍数の前後を求めていますが、これは2,3以外の素数は6の倍数の前後に現れる(6n±1)というものでmod 6したときの結果が0, 2, 3, 4は素数でなく残りが1, 5(6n±1)であることから説明がつくと思います。 実際に確認してみたところ大丈夫そうでした。
#[test]
fn primes_check(){
println!("{:?}", primes(59));
}
[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59]
素数の一覧を取得できたのでこれを元に素数p,qを以下のように求めておきます。
use rand::Rng;
let prime_vec = &primes(100)[5..];
let vec_len = prime_vec.len();
let mut rng = rand::thread_rng();
let mut p;
let mut q;
loop {
p = prime_vec[rng.gen::<usize>() % vec_len];
q = prime_vec[rng.gen::<usize>() % vec_len];
if p != 1 && q != 1 && p != q {
break;
}
}
実際はもっと大きな数値じゃないと使い物にならないはずですが、今回は確認ということで小さい数値で確認してみたいと思います。
N = pq を計算する
これはただ掛け合わせるだけです。実際は2048bitの数値とかを扱うので大きな数値の計算を考慮する必要がありそうです。
let N = p * q;
GCD(φ(N), e) = 1となるeを選択する
ここでは以下を確認してみたいと思います。 - φ(N) = φ(pq) = φ(p)φ(q) = (p-1)(q-1)となるか確認(オイラーの関数の乗法性) - GCD(φ(N), e) = 1 -> (p-1)(q-1) に対して素なeを拡張ユークリッドの互除法で求める
φ(p)=p-1となることはpが素数であることより自明ですが今回は確認のために地道に求めてみたいと思います。まず最小公倍数の確認として以下のユークリッド関数を作成しました。(拡張ユークリッドとは別です)
// ユークリッド互除法で最大公約数を計算
fn euclid(a: i32, b: i32) -> i32 {
let (small, big) = match a > b {
true => (b, a),
false => (a, b)
};
match big % small {
0 => small,
x => euclid(small, x)
}
}
テストコードを実行してみたところ、正しく動いてそうです。
#[test]
fn euclid_spec() {
assert_eq!(euclid(24, 40), 8);
assert_eq!(euclid(24, 15), 3);
assert_eq!(euclid(32, 40), 8);
assert_eq!(euclid(16, 21), 1);
}
次にオイラー関数を以下のように実装しました。やっていることとしては1 ~ aの範囲の数値でaとの最小公倍数が1(互いに素)なものをVectorに追加していって長さとともに返すだけです。
fn euler_func(a: i32) -> (usize, Vec<i32>) {
let mut vec: Vec<i32> = Vec::new();
for i in 1..a {
if euclid(i, a) == 1 {
vec.push(i);
}
}
(vec.len(), vec)
}
φ(p)=p-1となることについては以下で確認ができます。
#[test]
fn euler_func_spec() {
let prime_vec = &primes(100)[5..];
for prime in prime_vec {
let (len, _) = euler_func(*prime);
assert_eq!(len as i32, *prime-1);
}
}
次に乗法性の確認にとして以下を実行します。
#[test]
fn euler_func_multi_spec() {
let prime_vec = &primes(100)[5..];
for i in 0..prime_vec.len()-1 {
for j in i+1..prime_vec.len() {
let (len, _) = euler_func(prime_vec[i]*prime_vec[j]);
assert_eq!(len as i32, (prime_vec[i]-1)*(prime_vec[j]-1));
}
}
}
それから GCD(φ(N), e) = 1 -> (p-1)(q-1) に対して素なeを見つけるについては拡張ユークリッドの互除法を使うのですが以下のように実装しました。
fn ext_euclid(r1: i32, r2: i32) -> (i32, i32) {
fn ext_euclid_temp(r1: i32, r2: i32, u1: i32, u2: i32, v1: i32, v2: i32) -> (i32, i32) {
match r1 % r2 {
0 => (v1, v2),
n => {
let m = r1 / r2;
ext_euclid_temp(r2, n, v1, v2, u1 - m * v1, u2 - m * v2)
}
}
}
ext_euclid_temp(r1, r2, 1, 0, 0, 1)
}
拡張ユークリッドの互除法では2つの数x,yを渡したときにax + byが最大公約数となるようなa,bを出力するので以下の関数で確認しました。
#[test]
fn ext_euclid_spec() {
fn check(x: i32, y: i32) -> bool {
let (a, b) = ext_euclid(x, y);
a * x + b * y == euclid(x, y)
}
for x in 1..1000 {
for y in 1..1000 {
assert!(check(x, y));
}
}
}
GCD(φ(N), e) = 1となるeを見つけるについてe=3や65537を固定で使うといった手法もあるようですが、今回は以下のようにランダムに求めてみたいと思います。
let euler = (p-1) * (q-1);
let mut e:i32 = rng.gen::<i32>() % euler;
if e < 0 { e += euler;}
loop {
if e > 1 && euclid(euler, e) == 1 {
break;
} else {
e = rng.gen::<i32>() % euler;
if e < 0 { e += euler;}
}
};
ed = 1 mod φ(N) を満たすdを計算する
先に求めた拡張ユークリッドの関数に対して互いに素なeとφ(N)を渡すとa * d + b * φ(N) = 1となるa,bを返すのでこれを使います。ただこのとき返してくる数値が負の数の場合があるので、その時はφ(N)を足して0よりも大きくなるようにします。
let euler = (p-1) * (q-1);
let (mut d, _) = ext_euclid(e, euler);
if d < 0 {
d += euler;
}
公開鍵(N, e)、秘密鍵(d)を返す
最後に公開鍵(N, e)、秘密鍵(d)を返します。関数全体は以下のようになります。
fn rsa_keygen(k: i32) -> ((i32,i32), i32, (i32, i32)) {
let prime_vec = &primes(k)[5..];
println!("{:?}", prime_vec);
let vec_len = prime_vec.len();
let mut rng = rand::thread_rng();
let mut p;
let mut q;
loop {
p = prime_vec[rng.gen::<usize>() % vec_len];
q = prime_vec[rng.gen::<usize>() % vec_len];
if p != 1 && q != 1 && p != q {
break;
}
}
let mut rng = rand::thread_rng();
let N = p * q;
let euler = (p-1) * (q-1);
let mut e:i32 = rng.gen::<i32>() % euler;
if e < 0 { e += euler;}
loop {
if e > 1 && euclid(euler, e) == 1 {
break;
} else {
e = rng.gen::<i32>() % euler;
if e < 0 { e += euler;}
}
};
let (mut d, _) = ext_euclid(e, euler);
if d < 0 {
d += euler;
}
((N,e), d, (p, q))
}
今回はp, qの上限値を指定するようにしていますが、ここは実用的にするにはbit数を指定できるようにした方が良いのかと思います。
公開鍵を使って暗号化
公開鍵、秘密鍵が生成されたので次は公開鍵で暗号化を試してみたいと思います。やることは単純で以下のように入力mをe乗して mod Nを求めたものになります。
c ≡ m^e mod N
以下のように実装します。ここでは最後にmod Nするのではなく毎回mod Nするようにしていますが、mを aN * b と表すと aN ≡ 0 mod N よりm2 ≡ (aN * b)m ≡ と(aN)m * bm ≡ 0 + bm ≡ mod Nとなるので毎回mod Nしても結果が変わらないことがわかると思います。
fn rsa_enc(N: i32, e: i32, m: i32) -> i32 {
let mut result = 1;
for _ in 0..e {
result = result * m % N;
}
result
}
秘密鍵を使って復号化
最後に秘密鍵を使っての復号化ですが、暗号化の出力結果cをd乗してmod Nしたものが暗号化前の入力mと同一になります。
m' ≡ c^d mod N
関数を分けていますが暗号化で使ったものと同じ内容です。
fn rsa_dec(N: i32, d: i32, c: i32) -> i32 {
let mut result = 1;
for _ in 0..d {
result = result * c % N;
}
result
}
全体としては以下の通りになり、テストコードもパスすることが確認できます。
#[test]
fn check_rsa() {
use rand::Rng;
fn rsa_keygen(k: i32) -> ((i32,i32), i32, (i32, i32)) {
let prime_vec = &primes(k)[5..];
let vec_len = prime_vec.len();
let mut rng = rand::thread_rng();
let mut p;
let mut q;
loop {
p = prime_vec[rng.gen::<usize>() % vec_len];
q = prime_vec[rng.gen::<usize>() % vec_len];
if p != 1 && q != 1 && p != q {
break;
}
}
let mut rng = rand::thread_rng();
let N = p * q;
let euler = (p-1) * (q-1);
let mut e:i32 = rng.gen::<i32>() % euler;
if e < 0 { e += euler;}
loop {
if e > 1 && euclid(euler, e) == 1 {
break;
} else {
e = rng.gen::<i32>() % euler;
if e < 0 { e += euler;}
}
};
let (mut d, _) = ext_euclid(e, euler);
if d < 0 {
d += euler;
}
((N,e), d, (p, q))
}
fn rsa_enc(N: i32, e: i32, m: i32) -> i32 {
let mut result = 1;
for _ in 0..e {
result = result * m % N;
}
result
}
fn rsa_dec(N: i32, d: i32, c: i32) -> i32 {
let mut result = 1;
for _ in 0..d {
result = result * c % N;
}
result
}
let ((N, e), d, (p, q)) = rsa_keygen(100);
println!("(p,q) = {:?}", (p,q));
println!("((N,e), d) = {:?}", ((N,e), d));
for i in 0..(if N < 254 {N} else {254}) {
let original = i;
let enc = rsa_enc(N, e, original);
let dec = rsa_dec(N, d, enc);
println!("(original, enc, dec) = {:?}", (original, enc, dec));
assert_eq!(original, dec);
}
}
実際に試してみて大丈夫そうな感じではありましたが、念の為証明もしてみたいと思います。確認したい内容は復号化結果m'が入力mと同一になることなので以下の証明になります。
m' ≡ c^d ≡ (m^e)^d ≡ m^(ed) ≡ m^(k*φ(N) + 1) ≡ m^(k*φ(N)) * m
ed = 1 mod φ(N)より ed = k*φ(N)+1としています。
これについてGCD(m,N) = 1とGCD(m,N) = pまたはqの場合で証明できれば十分になります。
GCD(m,N) = 1について
オイラーの定理ではaと互いなnについてmφ(n) = 1 mod nが成り立ちます。これより以下が成り立ちます。
m' ≡ m^(k*φ(N)) * m ≡ 1^k * m ≡ m mod N
GCD(m,N) = pについて
pとqは互いに素なので拡張ユークリッドの互除法より以下の条件を満たすx, yが存在します。
p*x + q*y = 1
それからmedの法をNでなくp, qで評価すると
m^(ed) ≡ 0^(ed) ≡ 0 mod p m^(ed) ≡ m^(k*φ(N) + 1) ≡ m^(k(p-1)(q-1) + 1) ≡ m * (m^(q-1))^(k(p-1)) ≡ m * 1^(k(p-1)) ≡ m mod q
となります。ここで m * (mq-1)k(p-1) ≡ m * 1k(p-1)の変換については、qが素数でaがqの倍数でないときにaq-1 ≡ 1 mod qが成り立つというフェルマーの小定理を利用しています。
次にp,qが互いに素でa ≡ 0 mod p, b ≡ m mod qが成り立つので中国の剰余定理より以下のようになります。
m^(ed) ≡ a*q*y + b*p*x mod pq
これに対してGCD(m,N) = pを利用して式を進めると
m^(ed) ≡ a*q*y + b*p*x mod pq ≡ 0*q*y + m*p*x mod N ≡ m*p*x mod N ≡ m(1 - q*y) mod N ≡ m - m*q*y mod N ≡ m mod N
となりGCD(m,N) = pの場合もmed ≡ m mod Nが成り立つことがわかります。このときmqy mod N≡0 mod Nの変換はmがpの倍数であるためm*qがNの倍数になることを利用しています。
冪乗余の計算効率化
今回の実装では暗号化、復号化のための冪乗余の計算で以下のように素直に掛け合わせていましたが大きい数を扱う場合時間がかかってしまいます。
fn rsa_enc(N: i32, e: i32, m: i32) -> i32 {
let mut result = 1;
for _ in 0..e {
result = result * m % N;
}
result
}
冪乗余の計算高速化の手法としてバイナリ法というものがあります。
