Diffie-Hellman鍵交換プロトコル：理論とPython実装

import random
import time

random.seed(20230907)  # 再現性のため乱数シードを固定

p = 23  # 大きな素数（教科書的な例のため小さい値を使用）
g = 5   # p を法とする原始根

alice_private_key = random.randint(1, p - 2)
bob_private_key = random.randint(1, p - 2)

alice_public_key = pow(g, alice_private_key, p)  # A = g^a mod p
bob_public_key = pow(g, bob_private_key, p)      # B = g^b mod p

shared_secret_alice = pow(bob_public_key, alice_private_key, p)  # S_A = B^a mod p
shared_secret_bob = pow(alice_public_key, bob_private_key, p)    # S_B = A^b mod p

print(f"公開パラメータ: p = {p}, g = {g}")
print(f"アリスの秘密鍵 a = {alice_private_key}, 公開鍵 A = {alice_public_key}")
print(f"ボブの秘密鍵   b = {bob_private_key}, 公開鍵 B = {bob_public_key}")
print(f"アリス側の共通鍵 S_A = B^a mod p = {shared_secret_alice}")
print(f"ボブ側の共通鍵   S_B = A^b mod p = {shared_secret_bob}")
assert shared_secret_alice == shared_secret_bob
print(f"一致確認: S_A == S_B -> {shared_secret_alice}")


def brute_force_dlog(g, A, p):
    """g, A, p から離散対数 a を総当たりで探す（教育目的。小さい p でのみ現実的）"""
    for x in range(1, p - 1):
        if pow(g, x, p) == A:
            return x
    return None


t0 = time.perf_counter()
recovered_a = brute_force_dlog(g, alice_public_key, p)
t1 = time.perf_counter()
print(f"\nEve が総当たりで復元した a = {recovered_a}（正解: {alice_private_key}）")
print(f"総当たりに要した時間: {(t1 - t0) * 1e6:.2f} us（探索空間: p-2 = {p - 2} 通り）")

# a さえ割り出せば、CDH問題（g^{ab} mod p の計算）も自明に解ける
eve_shared_secret = pow(bob_public_key, recovered_a, p)
print(f"Eve が計算した共通鍵: {eve_shared_secret}（アリス・ボブの鍵と一致: {eve_shared_secret == shared_secret_alice}）")

公開パラメータ: p = 23, g = 5
アリスの秘密鍵 a = 4, 公開鍵 A = 4
ボブの秘密鍵   b = 8, 公開鍵 B = 16
アリス側の共通鍵 S_A = B^a mod p = 9
ボブ側の共通鍵   S_B = A^b mod p = 9
一致確認: S_A == S_B -> 9

Eve が総当たりで復元した a = 4（正解: 4）
総当たりに要した時間: 1.83 us（探索空間: p-2 = 21 通り）
Eve が計算した共通鍵: 9（アリス・ボブの鍵と一致: True）

import random

random.seed(20230907)

p = 23
g = 5

# 正規の当事者
alice_private = random.randint(1, p - 2)
bob_private = random.randint(1, p - 2)
alice_public = pow(g, alice_private, p)  # A = g^a mod p
bob_public = pow(g, bob_private, p)      # B = g^b mod p

# Eve は Alice-Eve 間、Eve-Bob 間で別々の秘密指数を用意する
eve_secret_vs_alice = random.randint(1, p - 2)  # Bob になりすます際に使う e_B
eve_secret_vs_bob = random.randint(1, p - 2)    # Alice になりすます際に使う e_A
fake_B = pow(g, eve_secret_vs_alice, p)  # Alice に「Bobの公開鍵」として送る
fake_A = pow(g, eve_secret_vs_bob, p)    # Bob に「Aliceの公開鍵」として送る

print(f"公開パラメータ: p = {p}, g = {g}")
print(f"Alice: a = {alice_private}, A = {alice_public}")
print(f"Bob  : b = {bob_private}, B = {bob_public}")
print(f"Eve  : e_B(対Alice用) = {eve_secret_vs_alice} -> 偽B' = {fake_B}")
print(f"Eve  : e_A(対Bob用)   = {eve_secret_vs_bob} -> 偽A' = {fake_A}")

# Alice は fake_B を「Bobの公開鍵」だと思って共通鍵を計算する
K_Alice_Eve = pow(fake_B, alice_private, p)                   # (g^{e_B})^a mod p
K_Eve_with_Alice = pow(alice_public, eve_secret_vs_alice, p)  # (g^a)^{e_B} mod p（Eve側の計算）

# Bob は fake_A を「Aliceの公開鍵」だと思って共通鍵を計算する
K_Bob_Eve = pow(fake_A, bob_private, p)                # (g^{e_A})^b mod p
K_Eve_with_Bob = pow(bob_public, eve_secret_vs_bob, p)  # (g^b)^{e_A} mod p（Eve側の計算）

print(f"\nAlice が確立したと思っている共通鍵 K_AE = {K_Alice_Eve}")
print(f"Eve が Alice との間で計算した共通鍵        = {K_Eve_with_Alice}（一致: {K_Alice_Eve == K_Eve_with_Alice}）")
print(f"Bob が確立したと思っている共通鍵   K_EB = {K_Bob_Eve}")
print(f"Eve が Bob との間で計算した共通鍵          = {K_Eve_with_Bob}（一致: {K_Bob_Eve == K_Eve_with_Bob}）")
print(f"\nK_AE == K_EB か? {K_Alice_Eve == K_Bob_Eve}")

公開パラメータ: p = 23, g = 5
Alice: a = 4, A = 4
Bob  : b = 8, B = 16
Eve  : e_B(対Alice用) = 17 -> 偽B' = 15
Eve  : e_A(対Bob用)   = 1 -> 偽A' = 5

Alice が確立したと思っている共通鍵 K_AE = 2
Eve が Alice との間で計算した共通鍵        = 2（一致: True）
Bob が確立したと思っている共通鍵   K_EB = 16
Eve が Bob との間で計算した共通鍵          = 16（一致: True）

K_AE == K_EB か? False

import random
from sympy import isprime, factorint, primitive_root

random.seed(20230907)

# p-1 が複数の小さな素因数を持つ素数を用意する
p = 2003
assert isprime(p)
print(f"p = {p}, p-1 の素因数分解: {factorint(p - 1)}")

g = primitive_root(p)
print(f"原始根 g = {g}")

# 攻撃者は「公開鍵」と称して、位数 t=7 の元を送りつける
t = 7
assert (p - 1) % t == 0
h = pow(g, (p - 1) // t, p)  # 位数ちょうど t の元
print(f"位数 {t} の元 h = g^((p-1)/{t}) mod p = {h}")
order_check = [pow(h, i, p) for i in range(1, t + 1)]
print(f"h, h^2, ..., h^{t} mod p = {order_check}（h^{t} = 1 のはず）")

# Bob は h を「相手の公開鍵」として受け取り、部分群チェックをせずに共通鍵を計算する
distinct_secrets = set()
trials = 2000
for _ in range(trials):
    b = random.randint(1, p - 2)  # Bob の本物の秘密鍵（毎回ランダム）
    S_B = pow(h, b, p)
    distinct_secrets.add(S_B)

print(f"\nBob の秘密鍵 b を {trials} 回変えて共通鍵 S_B = h^b mod p を計算")
print(f"実際に出現した S_B の相異なる値: {sorted(distinct_secrets)}")
print(f"個数: {len(distinct_secrets)}（高々 t={t} 通りに収まる）")

p = 2003, p-1 の素因数分解: {2: 1, 7: 1, 11: 1, 13: 1}
原始根 g = 5
位数 7 の元 h = g^((p-1)/7) mod p = 874
h, h^2, ..., h^7 mod p = [874, 733, 1685, 485, 1257, 974, 1]（h^7 = 1 のはず）

Bob の秘密鍵 b を 2000 回変えて共通鍵 S_B = h^b mod p を計算
実際に出現した S_B の相異なる値: [1, 485, 733, 874, 974, 1257, 1685]
個数: 7（高々 t=7 通りに収まる）

from sympy import factorint
from sympy.ntheory.modular import crt

# p-1 = 40 = 2^3 * 5 (smooth: 小さな素因数のみ)
p = 41
g = 6  # 原始根
x_true = 15  # 秘密の離散対数（Pohlig-Hellmanで復元する）
h = pow(g, x_true, p)
print(f"p = {p}, g = {g}, p-1 = {p - 1} の素因数分解 = {factorint(p - 1)}")
print(f"求めたい x: g^x = h mod p, h = {h}")

factors = factorint(p - 1)
remainders = []
moduli = []
for q, e in factors.items():
    qe = q**e
    g_sub = pow(g, (p - 1) // qe, p)  # 位数 qe の部分群に落とし込んだ g
    h_sub = pow(h, (p - 1) // qe, p)  # 位数 qe の部分群に落とし込んだ h
    x_i = None
    for candidate in range(qe):  # 小さな群なので総当たりで解ける
        if pow(g_sub, candidate, p) == h_sub:
            x_i = candidate
            break
    remainders.append(x_i)
    moduli.append(qe)
    print(f"  素因数冪 {qe:>3d}: 部分群での離散対数 x ≡ {x_i} (mod {qe})")

x_recovered, mod_total = crt(moduli, remainders)
print(f"\nCRT で復元した x = {x_recovered} (mod {mod_total})")
print(f"検算: g^x_recovered mod p = {pow(g, int(x_recovered), p)}（h = {h}）")
assert pow(g, int(x_recovered), p) == h
print(f"真の秘密指数と一致: {int(x_recovered) == x_true}")

p = 41, g = 6, p-1 = 40 の素因数分解 = {2: 3, 5: 1}
求めたい x: g^x = h mod p, h = 3

  素因数冪   8: 部分群での離散対数 x ≡ 7 (mod 8)
  素因数冪   5: 部分群での離散対数 x ≡ 0 (mod 5)

CRT で復元した x = 15 (mod 40)
検算: g^x_recovered mod p = 3（h = 3）
真の秘密指数と一致: True

import time
import random
from sympy import isprime, factorint, primitive_root, discrete_log, primerange

random.seed(20230907)

small_primes = list(primerange(2, 100))  # smooth性を保つための小さな素因数プール


def find_smooth_prime(target_bits, tries=20000):
    """p-1 の素因数がすべて100未満になるsmoothな素数を探す"""
    for _ in range(tries):
        random.shuffle(small_primes)
        prod = 1
        for q in small_primes:
            if prod.bit_length() >= target_bits:
                break
            prod *= q
        p = prod + 1
        if prod.bit_length() >= target_bits - 2 and isprime(p):
            return p
    raise RuntimeError("not found")


def find_safe_prime(bits):
    while True:
        q = random.getrandbits(bits - 1) | 1 | (1 << (bits - 2))
        if not isprime(q):
            continue
        p = 2 * q + 1
        if isprime(p):
            return p


for bits in [24, 32, 40, 48]:
    p_smooth = find_smooth_prime(bits)
    p_safe = find_safe_prime(bits)

    g_smooth = primitive_root(p_smooth)
    g_safe = primitive_root(p_safe)

    x_smooth = random.randint(2, p_smooth - 2)
    h_smooth = pow(g_smooth, x_smooth, p_smooth)
    x_safe = random.randint(2, p_safe - 2)
    h_safe = pow(g_safe, x_safe, p_safe)

    t0 = time.perf_counter()
    xr_smooth = discrete_log(p_smooth, h_smooth, g_smooth)
    t_smooth = time.perf_counter() - t0
    assert xr_smooth == x_smooth

    t0 = time.perf_counter()
    xr_safe = discrete_log(p_safe, h_safe, g_safe)
    t_safe = time.perf_counter() - t0
    assert xr_safe == x_safe

    max_factor_smooth = max(factorint(p_smooth - 1).keys())
    max_factor_safe_bits = max(factorint(p_safe - 1).keys()).bit_length()

    print(
        f"bit~{bits}: smooth p={p_smooth} (最大素因数 {max_factor_smooth}) -> {t_smooth * 1000:.3f} ms | "
        f"safe p={p_safe} (最大素因数 {max_factor_safe_bits}bit) -> {t_safe * 1000:.3f} ms | "
        f"倍率 {t_safe / t_smooth:.1f}x"
    )

bit~24: smooth p=16297859 (最大素因数 89) -> 0.107 ms | safe p=11600867 (最大素因数 23bit) -> 0.365 ms | 倍率 3.4x
bit~32: smooth p=10747150679 (最大素因数 89) -> 0.194 ms | safe p=3978555587 (最大素因数 31bit) -> 9.649 ms | 倍率 49.8x
bit~40: smooth p=1401201149063 (最大素因数 67) -> 0.262 ms | safe p=837601172687 (最大素因数 39bit) -> 173.684 ms | 倍率 664.1x
bit~48: smooth p=212996770326583 (最大素因数 71) -> 0.353 ms | safe p=253901762228183 (最大素因数 47bit) -> 16301.000 ms | 倍率 46200.3x