【非対称鍵暗号】CPU・AVX2・GPUで挑むハッキング実験と量子コンピュータがRSAを一撃で破る理

/* rsa_manual_cpu.c より抜粋：試し割り法 */
static u64 trial_division(u64 n) {
    if (n % 2 == 0) return 2;
    for (u64 d = 3; d * d <= n; d += 2) {
        if (n % d == 0) return d;
    }
    return 0; /* n自体が素数 */
}

/* rsa_manual_cpu.c より抜粋：試し割り法 */
static u64 trial_division(u64 n) {
    if (n % 2 == 0) return 2;
    for (u64 d = 3; d * d <= n; d += 2) {
        if (n % d == 0) return d;
    }
    return 0; /* n自体が素数 */
}

/* rsa_manual_avx2.c より抜粋 */
static u64 trial_division_avx2(u64 n) {
    if (n % 2 == 0) return 2;
    double dn = (double)n;
    __m256d vn = _mm256_set1_pd(dn);
    __m256d vzero = _mm256_set1_pd(0.0);

    for (u64 d = 3; (double)d * (double)d <= dn; d += 8) {
        double arr[4] = { (double)d, (double)(d + 2), (double)(d + 4), (double)(d + 6) };
        __m256d vd   = _mm256_loadu_pd(arr);
        __m256d q    = _mm256_div_pd(vn, vd);      // 除算
        __m256d qf   = _mm256_floor_pd(q);        // 床関数
        __m256d prod = _mm256_mul_pd(qf, vd);
        __m256d rem  = _mm256_sub_pd(vn, prod);     // 剰余の算出
        __m256d cmp  = _mm256_cmp_pd(rem, vzero, _CMP_EQ_OQ);
        int mask = _mm256_movemask_pd(cmp);
        if (mask) {
            for (int i = 0; i < 4; i++) {
                if (mask & (1 << i)) return (u64)arr[i];
            }
        }
    }
    return 0;
}

/* rsa_manual_avx2.c より抜粋 */
static u64 trial_division_avx2(u64 n) {
    if (n % 2 == 0) return 2;
    double dn = (double)n;
    __m256d vn = _mm256_set1_pd(dn);
    __m256d vzero = _mm256_set1_pd(0.0);

    for (u64 d = 3; (double)d * (double)d <= dn; d += 8) {
        double arr[4] = { (double)d, (double)(d + 2), (double)(d + 4), (double)(d + 6) };
        __m256d vd   = _mm256_loadu_pd(arr);
        __m256d q    = _mm256_div_pd(vn, vd);      // 除算
        __m256d qf   = _mm256_floor_pd(q);        // 床関数
        __m256d prod = _mm256_mul_pd(qf, vd);
        __m256d rem  = _mm256_sub_pd(vn, prod);     // 剰余の算出
        __m256d cmp  = _mm256_cmp_pd(rem, vzero, _CMP_EQ_OQ);
        int mask = _mm256_movemask_pd(cmp);
        if (mask) {
            for (int i = 0; i < 4; i++) {
                if (mask & (1 << i)) return (u64)arr[i];
            }
        }
    }
    return 0;
}

/* rsa_manual_avx512.c より抜粋（高精度ハイブリッド判定） */
static u64 trial_division_avx512(u64 n) {
    if (n % 2 == 0) return 2;
    double dn = (double)n;
    __m512d vn = _mm512_set1_pd(dn);
    __m512i vn_64 = _mm512_set1_epi64((long long)n);
    __m512i vzero_i = _mm512_setzero_si512();

    for (u64 d = 3; (double)d * (double)d <= dn; d += 16) {
        double arr[8] = {
            (double)d, (double)(d + 2), (double)(d + 4), (double)(d + 6),
            (double)(d + 8), (double)(d + 10), (double)(d + 12), (double)(d + 14)
        };
        __m512d vd   = _mm512_loadu_pd(arr);
        __m512d q    = _mm512_div_pd(vn, vd);
        __m512d qf   = _mm512_floor_pd(q);

        // 32ビット整数に一度変換し、64ビットへ符号拡張（商qfと除数vdは共に32ビット素数以下なので完全に正確）
        __m256i vqf_32 = _mm512_cvttpd_epi32(qf);
        __m512i vqf = _mm512_cvtepi32_epi64(vqf_32);

        __m256i vd_32 = _mm512_cvttpd_epi32(vd);
        __m512i vd_64 = _mm512_cvtepi32_epi64(vd_32);

        // 正確な64ビット乗算
        __m512i vprod = _mm512_mul_epu32(vqf, vd_64);

        // vrem = n - qf * d
        __m512i vrem = _mm512_sub_epi64(vn_64, vprod);

        // 剰余が0であるか比較（マスク判定）
        __mmask8 cmp = _mm512_cmpeq_epi64_mask(vrem, vzero_i);
        if (cmp) {
            for (int i = 0; i < 8; i++) {
                if (cmp & (1 << i)) return (u64)arr[i];
            }
        }
    }
    return 0;
}

/* rsa_manual_avx512.c より抜粋（高精度ハイブリッド判定） */
static u64 trial_division_avx512(u64 n) {
    if (n % 2 == 0) return 2;
    double dn = (double)n;
    __m512d vn = _mm512_set1_pd(dn);
    __m512i vn_64 = _mm512_set1_epi64((long long)n);
    __m512i vzero_i = _mm512_setzero_si512();

    for (u64 d = 3; (double)d * (double)d <= dn; d += 16) {
        double arr[8] = {
            (double)d, (double)(d + 2), (double)(d + 4), (double)(d + 6),
            (double)(d + 8), (double)(d + 10), (double)(d + 12), (double)(d + 14)
        };
        __m512d vd   = _mm512_loadu_pd(arr);
        __m512d q    = _mm512_div_pd(vn, vd);
        __m512d qf   = _mm512_floor_pd(q);

        // 32ビット整数に一度変換し、64ビットへ符号拡張（商qfと除数vdは共に32ビット素数以下なので完全に正確）
        __m256i vqf_32 = _mm512_cvttpd_epi32(qf);
        __m512i vqf = _mm512_cvtepi32_epi64(vqf_32);

        __m256i vd_32 = _mm512_cvttpd_epi32(vd);
        __m512i vd_64 = _mm512_cvtepi32_epi64(vd_32);

        // 正確な64ビット乗算
        __m512i vprod = _mm512_mul_epu32(vqf, vd_64);

        // vrem = n - qf * d
        __m512i vrem = _mm512_sub_epi64(vn_64, vprod);

        // 剰余が0であるか比較（マスク判定）
        __mmask8 cmp = _mm512_cmpeq_epi64_mask(vrem, vzero_i);
        if (cmp) {
            for (int i = 0; i < 8; i++) {
                if (cmp & (1 << i)) return (u64)arr[i];
            }
        }
    }
    return 0;
}

/* rsa_manual_opencl.c：GPU側カーネル */
__kernel void trial_division(__global const ulong *n_ptr,
                              __global const ulong *base_ptr,
                              __global ulong *result) {
    ulong n = *n_ptr;
    ulong base = *base_ptr;
    size_t gid = get_global_id(0);
    ulong d = base + gid * 2; /* 奇数候補のみを並列走査 */
    if (d >= 3 && n % d == 0) {
        *result = d; /* 発見した因数を直接上書き */
    }
}

/* rsa_manual_opencl.c：GPU側カーネル */
__kernel void trial_division(__global const ulong *n_ptr,
                              __global const ulong *base_ptr,
                              __global ulong *result) {
    ulong n = *n_ptr;
    ulong base = *base_ptr;
    size_t gid = get_global_id(0);
    ulong d = base + gid * 2; /* 奇数候補のみを並列走査 */
    if (d >= 3 && n % d == 0) {
        *result = d; /* 発見した因数を直接上書き */
    }
}

/*
 * rsa_manual_cpu.c
 * 汎用CPU（スカラー演算）版：
 *   1. 試し割り法による素因数分解
 *   2. ミラー・ラビン素数判定
 *   3. 超小型RSAの鍵生成・暗号化・復号デモ
 *   4. 桁数を増やしていったときの分解時間の計測
 *
 * コンパイル (Windows / MinGW-w64, GCC):
 *   gcc -O2 -o rsa_manual_cpu.exe rsa_manual_cpu.c
 *
 * コンパイル (Linux):
 *   gcc -O2 -o rsa_manual_cpu rsa_manual_cpu.c
 */
#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <time.h>
/* 128bit中間値が欲しいので __int128 を使う（GCC拡張） */
typedef unsigned long long u64;
typedef unsigned __int128 u128;
/* --- モジュラー演算 --- */
static u64 mulmod(u64 a, u64 b, u64 m) {
    return (u64)(((u128)a * (u128)b) % m);
}
static u64 powmod(u64 base, u64 exp, u64 mod) {
    u64 result = 1;
    base %= mod;
    while (exp > 0) {
        if (exp & 1) result = mulmod(result, base, mod);
        base = mulmod(base, base, mod);
        exp >>= 1;
    }
    return result;
}
/* --- ミラー・ラビン素数判定（決定的：64bit範囲で十分な基底を使用） --- */
static int miller_rabin(u64 n) {
    if (n < 2) return 0;
    if (n % 2 == 0) return n == 2;
    u64 d = n - 1;
    int r = 0;
    while (d % 2 == 0) { d /= 2; r++; }
    /* 64bit全域で決定的に判定できることが知られている基底集合 */
    u64 witnesses[] = {2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37};
    for (size_t i = 0; i < sizeof(witnesses)/sizeof(witnesses[0]); i++) {
        u64 a = witnesses[i];
        if (a >= n) continue;
        u64 x = powmod(a, d, n);
        if (x == 1 || x == n - 1) continue;
        int composite = 1;
        for (int j = 0; j < r - 1; j++) {
            x = mulmod(x, x, n);
            if (x == n - 1) { composite = 0; break; }
        }
        if (composite) return 0;
    }
    return 1;
}
/* --- 拡張ユークリッドの互除法（逆元計算：秘密鍵d導出用） --- */
static long long ext_gcd(long long a, long long b, long long *x, long long *y) {
    if (b == 0) { *x = 1; *y = 0; return a; }
    long long x1, y1;
    long long g = ext_gcd(b, a % b, &x1, &y1);
    *x = y1;
    *y = x1 - (a / b) * y1;
    return g;
}
static u64 modinv(u64 e, u64 phi) {
    long long x, y;
    ext_gcd((long long)e, (long long)phi, &x, &y);
    long long r = x % (long long)phi;
    if (r < 0) r += phi;
    return (u64)r;
}
/* --- 試し割り法：素因数（の一つ)を返す。見つからなければ0 --- */
static u64 trial_division(u64 n) {
    if (n % 2 == 0) return 2;
    for (u64 d = 3; d * d <= n; d += 2) {
        if (n % d == 0) return d;
    }
    return 0; /* n自体が素数 */
}
/* --- 手作業で追いやすい超小型RSAデモ --- */
static void tiny_rsa_demo(void) {
    printf("=== 超小型RSAデモ（暗算追跡用） ===\n");
    u64 p = 61, q = 53;             /* テキスト定番の例 */
    u64 n = p * q;                  /* 3233 */
    u64 phi = (p - 1) * (q - 1);    /* 3120 */
    u64 e = 17;                     /* gcd(e, phi) = 1 を満たす小さい値 */
    u64 d = modinv(e, phi);
    printf("p=%llu, q=%llu, n=p*q=%llu\n", p, q, n);
    printf("phi(n)=(p-1)(q-1)=%llu\n", phi);
    printf("公開鍵 (e, n) = (%llu, %llu)\n", e, n);
    printf("秘密鍵 (d, n) = (%llu, %llu)\n", d, n);
    u64 m = 65; /* 平文（0〜n-1の任意の値） */
    u64 c = powmod(m, e, n);
    u64 m2 = powmod(c, d, n);
    printf("平文 m=%llu -> 暗号文 c=m^e mod n=%llu -> 復号 m'=c^d mod n=%llu\n", m, c, m2);
    printf("復号結果が平文と一致: %s\n\n", (m == m2) ? "成功" : "失敗");
}
/* --- 桁数を増やしたときの素因数分解時間の計測 --- */
static void benchmark_factoring(void) {
    printf("=== 試し割り法による素因数分解ベンチマーク（CPUスカラー版） ===\n");
    printf("%-20s %-20s %-15s\n", "合成数 n (2つの素数の積)", "見つかった因数", "所要時間(秒)");
    /* 桁数を段階的に増やした半素数を用意（p, qは実際に素数判定済み） */
    u64 primes_small[] = {97, 997, 9973, 99991, 999983, 9999991ULL, 99999989ULL};
    int cnt = sizeof(primes_small)/sizeof(primes_small[0]);
    for (int i = 0; i < cnt; i++) {
        u64 p = primes_small[i];
        u64 q = primes_small[i]; /* 近い大きさの2素数の積で半素数を作る（q側は少しずらす） */
        /* qを次に大きい素数に置き換える簡易処理 */
        u64 qq = p + 2;
        while (!miller_rabin(qq)) qq += 2;
        u64 n = p * qq;
        clock_t t0 = clock();
        u64 factor = trial_division(n);
        clock_t t1 = clock();
        double sec = (double)(t1 - t0) / CLOCKS_PER_SEC;
        printf("n=%-18llu factor=%-18llu %.6f\n", n, factor, sec);
    }
    printf("\n");
}
/* 巨大合成数（1回約5秒・2回で10秒負荷）ベンチマーク */
static void benchmark_factoring_20sec(void) {
    u64 p = 4294967291ULL;
    u64 q = 4294967291ULL;
    u64 n = p * q;
    int iterations = 2;
    printf("=== CPUスカラー版 10秒負荷（巨大合成数・2イテレーション） ===\n");
    printf("n = %llu (factor = %llu)\n", n, p);
    printf("イテレーション数: %d\n", iterations);
    printf("実行中...\n");
    clock_t t0 = clock();
    u64 factor = 0;
    for (int i = 0; i < iterations; i++) {
        factor = trial_division(n);
    }
    clock_t t1 = clock();
    double sec = (double)(t1 - t0) / CLOCKS_PER_SEC;
    printf("完了。見つかった因数: %llu\n", factor);
    printf("総所要時間: %.6f 秒 (1イテレーション平均: %.6f 秒)\n\n", sec, sec / iterations);
}
int main(void) {
    tiny_rsa_demo();
    benchmark_factoring();
    benchmark_factoring_20sec();
    return 0;
}

/*
 * rsa_manual_cpu.c
 * 汎用CPU（スカラー演算）版：
 *   1. 試し割り法による素因数分解
 *   2. ミラー・ラビン素数判定
 *   3. 超小型RSAの鍵生成・暗号化・復号デモ
 *   4. 桁数を増やしていったときの分解時間の計測
 *
 * コンパイル (Windows / MinGW-w64, GCC):
 *   gcc -O2 -o rsa_manual_cpu.exe rsa_manual_cpu.c
 *
 * コンパイル (Linux):
 *   gcc -O2 -o rsa_manual_cpu rsa_manual_cpu.c
 */
#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <time.h>
/* 128bit中間値が欲しいので __int128 を使う（GCC拡張） */
typedef unsigned long long u64;
typedef unsigned __int128 u128;
/* --- モジュラー演算 --- */
static u64 mulmod(u64 a, u64 b, u64 m) {
    return (u64)(((u128)a * (u128)b) % m);
}
static u64 powmod(u64 base, u64 exp, u64 mod) {
    u64 result = 1;
    base %= mod;
    while (exp > 0) {
        if (exp & 1) result = mulmod(result, base, mod);
        base = mulmod(base, base, mod);
        exp >>= 1;
    }
    return result;
}
/* --- ミラー・ラビン素数判定（決定的：64bit範囲で十分な基底を使用） --- */
static int miller_rabin(u64 n) {
    if (n < 2) return 0;
    if (n % 2 == 0) return n == 2;
    u64 d = n - 1;
    int r = 0;
    while (d % 2 == 0) { d /= 2; r++; }
    /* 64bit全域で決定的に判定できることが知られている基底集合 */
    u64 witnesses[] = {2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37};
    for (size_t i = 0; i < sizeof(witnesses)/sizeof(witnesses[0]); i++) {
        u64 a = witnesses[i];
        if (a >= n) continue;
        u64 x = powmod(a, d, n);
        if (x == 1 || x == n - 1) continue;
        int composite = 1;
        for (int j = 0; j < r - 1; j++) {
            x = mulmod(x, x, n);
            if (x == n - 1) { composite = 0; break; }
        }
        if (composite) return 0;
    }
    return 1;
}
/* --- 拡張ユークリッドの互除法（逆元計算：秘密鍵d導出用） --- */
static long long ext_gcd(long long a, long long b, long long *x, long long *y) {
    if (b == 0) { *x = 1; *y = 0; return a; }
    long long x1, y1;
    long long g = ext_gcd(b, a % b, &x1, &y1);
    *x = y1;
    *y = x1 - (a / b) * y1;
    return g;
}
static u64 modinv(u64 e, u64 phi) {
    long long x, y;
    ext_gcd((long long)e, (long long)phi, &x, &y);
    long long r = x % (long long)phi;
    if (r < 0) r += phi;
    return (u64)r;
}
/* --- 試し割り法：素因数（の一つ)を返す。見つからなければ0 --- */
static u64 trial_division(u64 n) {
    if (n % 2 == 0) return 2;
    for (u64 d = 3; d * d <= n; d += 2) {
        if (n % d == 0) return d;
    }
    return 0; /* n自体が素数 */
}
/* --- 手作業で追いやすい超小型RSAデモ --- */
static void tiny_rsa_demo(void) {
    printf("=== 超小型RSAデモ（暗算追跡用） ===\n");
    u64 p = 61, q = 53;             /* テキスト定番の例 */
    u64 n = p * q;                  /* 3233 */
    u64 phi = (p - 1) * (q - 1);    /* 3120 */
    u64 e = 17;                     /* gcd(e, phi) = 1 を満たす小さい値 */
    u64 d = modinv(e, phi);
    printf("p=%llu, q=%llu, n=p*q=%llu\n", p, q, n);
    printf("phi(n)=(p-1)(q-1)=%llu\n", phi);
    printf("公開鍵 (e, n) = (%llu, %llu)\n", e, n);
    printf("秘密鍵 (d, n) = (%llu, %llu)\n", d, n);
    u64 m = 65; /* 平文（0〜n-1の任意の値） */
    u64 c = powmod(m, e, n);
    u64 m2 = powmod(c, d, n);
    printf("平文 m=%llu -> 暗号文 c=m^e mod n=%llu -> 復号 m'=c^d mod n=%llu\n", m, c, m2);
    printf("復号結果が平文と一致: %s\n\n", (m == m2) ? "成功" : "失敗");
}
/* --- 桁数を増やしたときの素因数分解時間の計測 --- */
static void benchmark_factoring(void) {
    printf("=== 試し割り法による素因数分解ベンチマーク（CPUスカラー版） ===\n");
    printf("%-20s %-20s %-15s\n", "合成数 n (2つの素数の積)", "見つかった因数", "所要時間(秒)");
    /* 桁数を段階的に増やした半素数を用意（p, qは実際に素数判定済み） */
    u64 primes_small[] = {97, 997, 9973, 99991, 999983, 9999991ULL, 99999989ULL};
    int cnt = sizeof(primes_small)/sizeof(primes_small[0]);
    for (int i = 0; i < cnt; i++) {
        u64 p = primes_small[i];
        u64 q = primes_small[i]; /* 近い大きさの2素数の積で半素数を作る（q側は少しずらす） */
        /* qを次に大きい素数に置き換える簡易処理 */
        u64 qq = p + 2;
        while (!miller_rabin(qq)) qq += 2;
        u64 n = p * qq;
        clock_t t0 = clock();
        u64 factor = trial_division(n);
        clock_t t1 = clock();
        double sec = (double)(t1 - t0) / CLOCKS_PER_SEC;
        printf("n=%-18llu factor=%-18llu %.6f\n", n, factor, sec);
    }
    printf("\n");
}
/* 巨大合成数（1回約5秒・2回で10秒負荷）ベンチマーク */
static void benchmark_factoring_20sec(void) {
    u64 p = 4294967291ULL;
    u64 q = 4294967291ULL;
    u64 n = p * q;
    int iterations = 2;
    printf("=== CPUスカラー版 10秒負荷（巨大合成数・2イテレーション） ===\n");
    printf("n = %llu (factor = %llu)\n", n, p);
    printf("イテレーション数: %d\n", iterations);
    printf("実行中...\n");
    clock_t t0 = clock();
    u64 factor = 0;
    for (int i = 0; i < iterations; i++) {
        factor = trial_division(n);
    }
    clock_t t1 = clock();
    double sec = (double)(t1 - t0) / CLOCKS_PER_SEC;
    printf("完了。見つかった因数: %llu\n", factor);
    printf("総所要時間: %.6f 秒 (1イテレーション平均: %.6f 秒)\n\n", sec, sec / iterations);
}
int main(void) {
    tiny_rsa_demo();
    benchmark_factoring();
    benchmark_factoring_20sec();
    return 0;
}

/*
 * rsa_manual_avx2.c
 * AVX2版：倍精度浮動小数点(__m256d)で商を算出し、64bit整数ベクトル演算で剰余を
 * 正確に求めるハイブリッド試し割り法。
 * これにより倍精度の精度限界(2^53)を超えた64bit範囲(2^64-1)までの合成数に対して
 * 正確に動作する並列演算を実現します。
 *
 * コンパイル (Windows / MinGW-w64, GCC):
 *   gcc -O2 -mavx2 -o rsa_manual_avx2.exe rsa_manual_avx2.c
 *
 * コンパイル (Linux):
 *   gcc -O2 -mavx2 -o rsa_manual_avx2 rsa_manual_avx2.c
 */
#include <stdio.h>
#include <stdint.h>
#include <time.h>
#include <immintrin.h>
typedef unsigned long long u64;
/* --- スカラー版（比較用） --- */
static u64 trial_division_scalar(u64 n) {
    if (n % 2 == 0) return 2;
    for (u64 d = 3; d * d <= n; d += 2) {
        if (n % d == 0) return d;
    }
    return 0;
}
/* --- AVX2版：4候補同時チェック（高精度64bit整数剰余判定） --- */
static u64 trial_division_avx2(u64 n) {
    if (n % 2 == 0) return 2;
    double dn = (double)n;
    __m256d vn = _mm256_set1_pd(dn);
    __m256i vn_64 = _mm256_set1_epi64x((long long)n);
    __m256i vzero_i = _mm256_setzero_si256();
    for (u64 d = 3; (double)d * (double)d <= dn; d += 8) {
        double arr[4] = { (double)d, (double)(d + 2), (double)(d + 4), (double)(d + 6) };
        __m256d vd   = _mm256_loadu_pd(arr);
        __m256d q    = _mm256_div_pd(vn, vd);
        __m256d qf   = _mm256_floor_pd(q);

        // 32ビット整数に一度変換し、64ビットへ符号拡張（商qfと除数vdは共に32ビット素数以下なので完全に正確）
        __m128i vqf_32 = _mm256_cvttpd_epi32(qf);
        __m256i vqf = _mm256_cvtepi32_epi64(vqf_32);

        __m128i vd_32 = _mm256_cvttpd_epi32(vd);
        __m256i vd_64 = _mm256_cvtepi32_epi64(vd_32);
        // 正確な64ビット乗算（32ビット同士の積のため、_mm256_mul_epu32 で安全かつ操作）
        __m256i vprod = _mm256_mul_epu32(vqf, vd_64);

        // vrem = n - qf * d
        __m256i vrem = _mm256_sub_epi64(vn_64, vprod);

        // 剰余が0であるか比較
        __m256i vcmp = _mm256_cmpeq_epi64(vrem, vzero_i);
        int mask = _mm256_movemask_pd(_mm256_castsi256_pd(vcmp));
        if (mask) {
            for (int i = 0; i < 4; i++) {
                if (mask & (1 << i)) return (u64)arr[i];
            }
        }
    }
    return 0;
}
static void benchmark_20sec(void) {
    u64 p = 4294967291ULL;
    u64 q = p;
    u64 n = p * q;
    int iterations = 2;
    printf("\n=== スカラー版 vs AVX2版 10秒負荷（巨大合成数・2イテレーション） ===\n");
    printf("n = %llu (factor = %llu)\n", n, p);
    printf("実行中...\n");
    clock_t t0 = clock();
    u64 f1 = 0;
    for (int i = 0; i < iterations; i++) {
        f1 = trial_division_scalar(n);
    }
    clock_t t1 = clock();
    u64 f2 = 0;
    for (int i = 0; i < iterations; i++) {
        f2 = trial_division_avx2(n);
    }
    clock_t t2 = clock();
    double sec_scalar = (double)(t1 - t0) / CLOCKS_PER_SEC;
    double sec_avx2   = (double)(t2 - t1) / CLOCKS_PER_SEC;
    double ratio = (sec_avx2 > 0) ? (sec_scalar / sec_avx2) : 0.0;
    printf("スカラー総時間 : %.6f 秒 (平均: %.6f 秒)\n", sec_scalar, sec_scalar / iterations);
    printf("AVX2総時間    : %.6f 秒 (平均: %.6f 秒)\n", sec_avx2, sec_avx2 / iterations);
    printf("速度倍率       : %.2fx\n\n", ratio);
}
int main(void) {
    printf("=== スカラー版 vs AVX2版（4並列・高精度剰余）試し割り法ベンチマーク ===\n");
    printf("（64bit整数の最大範囲で正確に剰余判定可能）\n\n");
    printf("%-20s %-14s %-14s %-10s\n", "n", "scalar(秒)", "AVX2(秒)", "倍率");
    u64 primes[] = {1000003ULL, 10000019ULL, 30000001ULL, 60000011ULL, 90000049ULL};
    int cnt = sizeof(primes)/sizeof(primes[0]);
    for (int i = 0; i < cnt; i++) {
        u64 p = primes[i];
        u64 q = p;
        u64 n = p * q;
        clock_t t0 = clock();
        u64 f1 = trial_division_scalar(n);
        clock_t t1 = clock();
        u64 f2 = trial_division_avx2(n);
        clock_t t2 = clock();
        double sec_scalar = (double)(t1 - t0) / CLOCKS_PER_SEC;
        double sec_avx2   = (double)(t2 - t1) / CLOCKS_PER_SEC;
        double ratio = (sec_avx2 > 0) ? (sec_scalar / sec_avx2) : 0.0;
        printf("n=%-18llu %-14.6f %-14.6f %.2fx  (factor scalar=%llu, avx2=%llu)\n",
               n, sec_scalar, sec_avx2, ratio, f1, f2);
    }
    benchmark_20sec();
    return 0;
}

/*
 * rsa_manual_avx2.c
 * AVX2版：倍精度浮動小数点(__m256d)で商を算出し、64bit整数ベクトル演算で剰余を
 * 正確に求めるハイブリッド試し割り法。
 * これにより倍精度の精度限界(2^53)を超えた64bit範囲(2^64-1)までの合成数に対して
 * 正確に動作する並列演算を実現します。
 *
 * コンパイル (Windows / MinGW-w64, GCC):
 *   gcc -O2 -mavx2 -o rsa_manual_avx2.exe rsa_manual_avx2.c
 *
 * コンパイル (Linux):
 *   gcc -O2 -mavx2 -o rsa_manual_avx2 rsa_manual_avx2.c
 */
#include <stdio.h>
#include <stdint.h>
#include <time.h>
#include <immintrin.h>
typedef unsigned long long u64;
/* --- スカラー版（比較用） --- */
static u64 trial_division_scalar(u64 n) {
    if (n % 2 == 0) return 2;
    for (u64 d = 3; d * d <= n; d += 2) {
        if (n % d == 0) return d;
    }
    return 0;
}
/* --- AVX2版：4候補同時チェック（高精度64bit整数剰余判定） --- */
static u64 trial_division_avx2(u64 n) {
    if (n % 2 == 0) return 2;
    double dn = (double)n;
    __m256d vn = _mm256_set1_pd(dn);
    __m256i vn_64 = _mm256_set1_epi64x((long long)n);
    __m256i vzero_i = _mm256_setzero_si256();
    for (u64 d = 3; (double)d * (double)d <= dn; d += 8) {
        double arr[4] = { (double)d, (double)(d + 2), (double)(d + 4), (double)(d + 6) };
        __m256d vd   = _mm256_loadu_pd(arr);
        __m256d q    = _mm256_div_pd(vn, vd);
        __m256d qf   = _mm256_floor_pd(q);

        // 32ビット整数に一度変換し、64ビットへ符号拡張（商qfと除数vdは共に32ビット素数以下なので完全に正確）
        __m128i vqf_32 = _mm256_cvttpd_epi32(qf);
        __m256i vqf = _mm256_cvtepi32_epi64(vqf_32);

        __m128i vd_32 = _mm256_cvttpd_epi32(vd);
        __m256i vd_64 = _mm256_cvtepi32_epi64(vd_32);
        // 正確な64ビット乗算（32ビット同士の積のため、_mm256_mul_epu32 で安全かつ操作）
        __m256i vprod = _mm256_mul_epu32(vqf, vd_64);

        // vrem = n - qf * d
        __m256i vrem = _mm256_sub_epi64(vn_64, vprod);

        // 剰余が0であるか比較
        __m256i vcmp = _mm256_cmpeq_epi64(vrem, vzero_i);
        int mask = _mm256_movemask_pd(_mm256_castsi256_pd(vcmp));
        if (mask) {
            for (int i = 0; i < 4; i++) {
                if (mask & (1 << i)) return (u64)arr[i];
            }
        }
    }
    return 0;
}
static void benchmark_20sec(void) {
    u64 p = 4294967291ULL;
    u64 q = p;
    u64 n = p * q;
    int iterations = 2;
    printf("\n=== スカラー版 vs AVX2版 10秒負荷（巨大合成数・2イテレーション） ===\n");
    printf("n = %llu (factor = %llu)\n", n, p);
    printf("実行中...\n");
    clock_t t0 = clock();
    u64 f1 = 0;
    for (int i = 0; i < iterations; i++) {
        f1 = trial_division_scalar(n);
    }
    clock_t t1 = clock();
    u64 f2 = 0;
    for (int i = 0; i < iterations; i++) {
        f2 = trial_division_avx2(n);
    }
    clock_t t2 = clock();
    double sec_scalar = (double)(t1 - t0) / CLOCKS_PER_SEC;
    double sec_avx2   = (double)(t2 - t1) / CLOCKS_PER_SEC;
    double ratio = (sec_avx2 > 0) ? (sec_scalar / sec_avx2) : 0.0;
    printf("スカラー総時間 : %.6f 秒 (平均: %.6f 秒)\n", sec_scalar, sec_scalar / iterations);
    printf("AVX2総時間    : %.6f 秒 (平均: %.6f 秒)\n", sec_avx2, sec_avx2 / iterations);
    printf("速度倍率       : %.2fx\n\n", ratio);
}
int main(void) {
    printf("=== スカラー版 vs AVX2版（4並列・高精度剰余）試し割り法ベンチマーク ===\n");
    printf("（64bit整数の最大範囲で正確に剰余判定可能）\n\n");
    printf("%-20s %-14s %-14s %-10s\n", "n", "scalar(秒)", "AVX2(秒)", "倍率");
    u64 primes[] = {1000003ULL, 10000019ULL, 30000001ULL, 60000011ULL, 90000049ULL};
    int cnt = sizeof(primes)/sizeof(primes[0]);
    for (int i = 0; i < cnt; i++) {
        u64 p = primes[i];
        u64 q = p;
        u64 n = p * q;
        clock_t t0 = clock();
        u64 f1 = trial_division_scalar(n);
        clock_t t1 = clock();
        u64 f2 = trial_division_avx2(n);
        clock_t t2 = clock();
        double sec_scalar = (double)(t1 - t0) / CLOCKS_PER_SEC;
        double sec_avx2   = (double)(t2 - t1) / CLOCKS_PER_SEC;
        double ratio = (sec_avx2 > 0) ? (sec_scalar / sec_avx2) : 0.0;
        printf("n=%-18llu %-14.6f %-14.6f %.2fx  (factor scalar=%llu, avx2=%llu)\n",
               n, sec_scalar, sec_avx2, ratio, f1, f2);
    }
    benchmark_20sec();
    return 0;
}

/*
 * rsa_manual_avx512.c
 * AVX-512版：倍精度浮動小数点(__m512d)で商を算出し、64bit整数ベクトル演算で剰余を
 * 正確に求めるハイブリッド試し割り法。
 * これにより倍精度の精度限界(2^53)を超えた64bit範囲(2^64-1)までの合成数に対して
 * 正確に動作する並列演算を実現します。
 *
 * コンパイル (Windows / MinGW-w64, GCC):
 *   gcc -O2 -mavx512f -o rsa_manual_avx512.exe rsa_manual_avx512.c
 */
#include <stdio.h>
#include <stdint.h>
#include <time.h>
#include <immintrin.h>
typedef unsigned long long u64;
/* --- スカラー版（比較用） --- */
static u64 trial_division_scalar(u64 n) {
    if (n % 2 == 0) return 2;
    for (u64 d = 3; d * d <= n; d += 2) {
        if (n % d == 0) return d;
    }
    return 0;
}
/* --- AVX-512版：8候補同時チェック（高精度64bit整数剰余判定） --- */
static u64 trial_division_avx512(u64 n) {
    if (n % 2 == 0) return 2;
    double dn = (double)n;
    __m512d vn = _mm512_set1_pd(dn);
    __m512i vn_64 = _mm512_set1_epi64((long long)n);
    __m512i vzero_i = _mm512_setzero_si512();
    for (u64 d = 3; (double)d * (double)d <= dn; d += 16) {
        double arr[8] = {
            (double)d, (double)(d + 2), (double)(d + 4), (double)(d + 6),
            (double)(d + 8), (double)(d + 10), (double)(d + 12), (double)(d + 14)
        };
        __m512d vd   = _mm512_loadu_pd(arr);
        __m512d q    = _mm512_div_pd(vn, vd);
        __m512d qf   = _mm512_floor_pd(q);

        // 32ビット整数に一度変換し、64ビットへ符号拡張（商qfと除数vdは共に32ビット素数以下なので完全に正確）
        __m256i vqf_32 = _mm512_cvttpd_epi32(qf);
        __m512i vqf = _mm512_cvtepi32_epi64(vqf_32);

        __m256i vd_32 = _mm512_cvttpd_epi32(vd);
        __m512i vd_64 = _mm512_cvtepi32_epi64(vd_32);
        // 正確な64ビット乗算（32ビット同士の積のため、_mm512_mul_epu32 で安全かつ正確）
        __m512i vprod = _mm512_mul_epu32(vqf, vd_64);

        // vrem = n - qf * d
        __m512i vrem = _mm512_sub_epi64(vn_64, vprod);

        // 剰余が0であるか比較（マスク判定）
        __mmask8 cmp = _mm512_cmpeq_epi64_mask(vrem, vzero_i);
        if (cmp) {
            for (int i = 0; i < 8; i++) {
                if (cmp & (1 << i)) return (u64)arr[i];
            }
        }
    }
    return 0;
}
static void benchmark_20sec(void) {
    u64 p = 4294967291ULL;
    u64 q = p;
    u64 n = p * q;
    int iterations = 2;
    printf("\n=== スカラー版 vs AVX-512版 10秒負荷（巨大合成数・2イテレーション） ===\n");
    printf("n = %llu (factor = %llu)\n", n, p);
    printf("実行中...\n");
    clock_t t0 = clock();
    u64 f1 = 0;
    for (int i = 0; i < iterations; i++) {
        f1 = trial_division_scalar(n);
    }
    clock_t t1 = clock();
    u64 f2 = 0;
    for (int i = 0; i < iterations; i++) {
        f2 = trial_division_avx512(n);
    }
    clock_t t2 = clock();
    double sec_scalar = (double)(t1 - t0) / CLOCKS_PER_SEC;
    double sec_avx512 = (double)(t2 - t1) / CLOCKS_PER_SEC;
    double ratio = (sec_avx512 > 0) ? (sec_scalar / sec_avx512) : 0.0;
    printf("スカラー総時間 : %.6f 秒 (平均: %.6f 秒)\n", sec_scalar, sec_scalar / iterations);
    printf("AVX-512総時間  : %.6f 秒 (平均: %.6f 秒)\n", sec_avx512, sec_avx512 / iterations);
    printf("速度倍率       : %.2fx\n\n", ratio);
}
int main(void) {
    printf("=== スカラー版 vs AVX-512版（8並列・高精度剰余）試し割り法ベンチマーク ===\n");
    printf("（64bit整数の最大範囲で正確に剰余判定可能）\n\n");
    printf("%-20s %-14s %-14s %-10s\n", "n", "scalar(秒)", "AVX-512(秒)", "倍率");
    u64 primes[] = {1000003ULL, 10000019ULL, 30000001ULL, 60000011ULL, 90000049ULL};
    int cnt = sizeof(primes)/sizeof(primes[0]);
    for (int i = 0; i < cnt; i++) {
        u64 p = primes[i];
        u64 q = p;
        u64 n = p * q;
        clock_t t0 = clock();
        u64 f1 = trial_division_scalar(n);
        clock_t t1 = clock();
        u64 f2 = trial_division_avx512(n);
        clock_t t2 = clock();
        double sec_scalar = (double)(t1 - t0) / CLOCKS_PER_SEC;
        double sec_avx512 = (double)(t2 - t1) / CLOCKS_PER_SEC;
        double ratio = (sec_avx512 > 0) ? (sec_scalar / sec_avx512) : 0.0;
        printf("n=%-18llu %-14.6f %-14.6f %.2fx  (factor scalar=%llu, avx512=%llu)\n",
               n, sec_scalar, sec_avx512, ratio, f1, f2);
    }
    benchmark_20sec();
    return 0;
}

/*
 * rsa_manual_avx512.c
 * AVX-512版：倍精度浮動小数点(__m512d)で商を算出し、64bit整数ベクトル演算で剰余を
 * 正確に求めるハイブリッド試し割り法。
 * これにより倍精度の精度限界(2^53)を超えた64bit範囲(2^64-1)までの合成数に対して
 * 正確に動作する並列演算を実現します。
 *
 * コンパイル (Windows / MinGW-w64, GCC):
 *   gcc -O2 -mavx512f -o rsa_manual_avx512.exe rsa_manual_avx512.c
 */
#include <stdio.h>
#include <stdint.h>
#include <time.h>
#include <immintrin.h>
typedef unsigned long long u64;
/* --- スカラー版（比較用） --- */
static u64 trial_division_scalar(u64 n) {
    if (n % 2 == 0) return 2;
    for (u64 d = 3; d * d <= n; d += 2) {
        if (n % d == 0) return d;
    }
    return 0;
}
/* --- AVX-512版：8候補同時チェック（高精度64bit整数剰余判定） --- */
static u64 trial_division_avx512(u64 n) {
    if (n % 2 == 0) return 2;
    double dn = (double)n;
    __m512d vn = _mm512_set1_pd(dn);
    __m512i vn_64 = _mm512_set1_epi64((long long)n);
    __m512i vzero_i = _mm512_setzero_si512();
    for (u64 d = 3; (double)d * (double)d <= dn; d += 16) {
        double arr[8] = {
            (double)d, (double)(d + 2), (double)(d + 4), (double)(d + 6),
            (double)(d + 8), (double)(d + 10), (double)(d + 12), (double)(d + 14)
        };
        __m512d vd   = _mm512_loadu_pd(arr);
        __m512d q    = _mm512_div_pd(vn, vd);
        __m512d qf   = _mm512_floor_pd(q);

        // 32ビット整数に一度変換し、64ビットへ符号拡張（商qfと除数vdは共に32ビット素数以下なので完全に正確）
        __m256i vqf_32 = _mm512_cvttpd_epi32(qf);
        __m512i vqf = _mm512_cvtepi32_epi64(vqf_32);

        __m256i vd_32 = _mm512_cvttpd_epi32(vd);
        __m512i vd_64 = _mm512_cvtepi32_epi64(vd_32);
        // 正確な64ビット乗算（32ビット同士の積のため、_mm512_mul_epu32 で安全かつ正確）
        __m512i vprod = _mm512_mul_epu32(vqf, vd_64);

        // vrem = n - qf * d
        __m512i vrem = _mm512_sub_epi64(vn_64, vprod);

        // 剰余が0であるか比較（マスク判定）
        __mmask8 cmp = _mm512_cmpeq_epi64_mask(vrem, vzero_i);
        if (cmp) {
            for (int i = 0; i < 8; i++) {
                if (cmp & (1 << i)) return (u64)arr[i];
            }
        }
    }
    return 0;
}
static void benchmark_20sec(void) {
    u64 p = 4294967291ULL;
    u64 q = p;
    u64 n = p * q;
    int iterations = 2;
    printf("\n=== スカラー版 vs AVX-512版 10秒負荷（巨大合成数・2イテレーション） ===\n");
    printf("n = %llu (factor = %llu)\n", n, p);
    printf("実行中...\n");
    clock_t t0 = clock();
    u64 f1 = 0;
    for (int i = 0; i < iterations; i++) {
        f1 = trial_division_scalar(n);
    }
    clock_t t1 = clock();
    u64 f2 = 0;
    for (int i = 0; i < iterations; i++) {
        f2 = trial_division_avx512(n);
    }
    clock_t t2 = clock();
    double sec_scalar = (double)(t1 - t0) / CLOCKS_PER_SEC;
    double sec_avx512 = (double)(t2 - t1) / CLOCKS_PER_SEC;
    double ratio = (sec_avx512 > 0) ? (sec_scalar / sec_avx512) : 0.0;
    printf("スカラー総時間 : %.6f 秒 (平均: %.6f 秒)\n", sec_scalar, sec_scalar / iterations);
    printf("AVX-512総時間  : %.6f 秒 (平均: %.6f 秒)\n", sec_avx512, sec_avx512 / iterations);
    printf("速度倍率       : %.2fx\n\n", ratio);
}
int main(void) {
    printf("=== スカラー版 vs AVX-512版（8並列・高精度剰余）試し割り法ベンチマーク ===\n");
    printf("（64bit整数の最大範囲で正確に剰余判定可能）\n\n");
    printf("%-20s %-14s %-14s %-10s\n", "n", "scalar(秒)", "AVX-512(秒)", "倍率");
    u64 primes[] = {1000003ULL, 10000019ULL, 30000001ULL, 60000011ULL, 90000049ULL};
    int cnt = sizeof(primes)/sizeof(primes[0]);
    for (int i = 0; i < cnt; i++) {
        u64 p = primes[i];
        u64 q = p;
        u64 n = p * q;
        clock_t t0 = clock();
        u64 f1 = trial_division_scalar(n);
        clock_t t1 = clock();
        u64 f2 = trial_division_avx512(n);
        clock_t t2 = clock();
        double sec_scalar = (double)(t1 - t0) / CLOCKS_PER_SEC;
        double sec_avx512 = (double)(t2 - t1) / CLOCKS_PER_SEC;
        double ratio = (sec_avx512 > 0) ? (sec_scalar / sec_avx512) : 0.0;
        printf("n=%-18llu %-14.6f %-14.6f %.2fx  (factor scalar=%llu, avx512=%llu)\n",
               n, sec_scalar, sec_avx512, ratio, f1, f2);
    }
    benchmark_20sec();
    return 0;
}

/*
 * rsa_manual_opencl.c
 * OpenCL版：多数のGPUスレッド（ワークアイテム）に候補除数を1つずつ割り当て、
 * 全候補を並列でチェックする試し割り法。
 *
 * コンパイル (Windows / MinGW-w64, GCC。OpenCL SDKのヘッダ・ライブラリが必要):
 *   gcc -O2 -I"<OpenCL_include>" -L"<OpenCL_lib>" -o rsa_manual_opencl.exe rsa_manual_opencl.c -lOpenCL
 *
 * または、カレントディレクトリにCLヘッダーを配置し、システムのDLLを直接リンクする場合:
 *   gcc -O2 -I. -o rsa_manual_opencl.exe rsa_manual_opencl.c C:\Windows\System32\OpenCL.dll
 *
 * コンパイル (Linux):
 *   gcc -O2 -o rsa_manual_opencl rsa_manual_opencl.c -lOpenCL
 *
 * 実行にはGPU（またはCPU上で動くOpenCLランタイム）が必要です。
 * 手元にNVIDIA/AMD/Intelいずれかの環境があれば、そのまま動作します。
 */
#define CL_TARGET_OPENCL_VERSION 120
#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <string.h>
#include <time.h>
#include <CL/cl.h>
typedef unsigned long long u64;
/* GPU側カーネル：各ワークアイテムが1つの候補除数dを担当し、n % d == 0 かを判定する
 * メモリ使用量を削減するため、発見した因数を結果バッファのアドレスに直接書き込む */
static const char *kernel_src =
"__kernel void trial_division(__global const ulong *n_ptr,\n"
"                              __global const ulong *base_ptr,\n"
"                              __global ulong *result) {\n"
"    ulong n = *n_ptr;\n"
"    ulong base = *base_ptr;\n"
"    size_t gid = get_global_id(0);\n"
"    ulong d = base + gid * 2; /* 奇数候補のみを走査 */\n"
"    if (d >= 3 && n % d == 0) {\n"
"        *result = d;\n;\n"
"    }\n"
"}\n";
int main(void) {
    cl_int err;
    cl_uint num_platforms = 0;
    clGetPlatformIDs(0, NULL, &num_platforms);
    if (num_platforms == 0) {
        fprintf(stderr, "OpenCLプラットフォームが見つかりません。"
                        "GPUドライバ（またはCPU用OpenCLランタイム）を導入してください。\n");
        return 1;
    }
    cl_platform_id platform;
    clGetPlatformIDs(1, &platform, NULL);
    cl_device_id device;
    err = clGetDeviceIDs(platform, CL_DEVICE_TYPE_GPU, 1, &device, NULL);
    if (err != CL_SUCCESS) {
        /* GPUが無ければCPUデバイスにフォールバック（動作確認用） */
        err = clGetDeviceIDs(platform, CL_DEVICE_TYPE_CPU, 1, &device, NULL);
        if (err != CL_SUCCESS) {
            fprintf(stderr, "使用可能なOpenCLデバイスが見つかりません。\n");
            return 1;
        }
        fprintf(stderr, "[注意] GPUが見つからないため、OpenCL経由のCPUデバイスで実行します。\n");
    }
    cl_context context = clCreateContext(NULL, 1, &device, NULL, NULL, &err);
    cl_command_queue queue = clCreateCommandQueue(context, device, 0, &err);
    cl_program program = clCreateProgramWithSource(context, 1, &kernel_src, NULL, &err);
    err = clBuildProgram(program, 1, &device, NULL, NULL, NULL);
    if (err != CL_SUCCESS) {
        char log[4096];
        clGetProgramBuildInfo(program, device, CL_PROGRAM_BUILD_LOG, sizeof(log), log, NULL);
        fprintf(stderr, "カーネルビルドエラー:\n%s\n", log);
        return 1;
    }
    cl_kernel kernel = clCreateKernel(program, "trial_division", &err);
    /* 対象n（4294967291 * 4294967291 = 18446744030759878681） */
    u64 n = 4294967291ULL * 4294967291ULL;
    u64 base = 3;
    size_t work_items = (size_t)(4294967291ULL / 2); /* 奇数候補数 = 約21億4700万 */
    int iterations = 2;
    cl_mem n_buf = clCreateBuffer(context, CL_MEM_READ_ONLY | CL_MEM_COPY_HOST_PTR, sizeof(u64), &n, &err);
    cl_mem base_buf = clCreateBuffer(context, CL_MEM_READ_ONLY | CL_MEM_COPY_HOST_PTR, sizeof(u64), &base, &err);
    cl_mem result_buf = clCreateBuffer(context, CL_MEM_READ_WRITE, sizeof(u64), NULL, &err);
    clSetKernelArg(kernel, 0, sizeof(cl_mem), &n_buf);
    clSetKernelArg(kernel, 1, sizeof(cl_mem), &base_buf);
    clSetKernelArg(kernel, 2, sizeof(cl_mem), &result_buf);
    printf("=== GPU (OpenCL) 10秒負荷（巨大合成数・2イテレーション） ===\n");
    printf("n = %llu (factor = 4294967291)\n", n);
    printf("ワークアイテム数: %zu\n", work_items);
    printf("実行中...\n");
    struct timespec t0, t1;
    clock_gettime(CLOCK_MONOTONIC, &t0);
    u64 found = 0;
    for (int iter = 0; iter < iterations; iter++) {
        // 結果バッファを 0 にリセット
        u64 zero = 0;
        clEnqueueWriteBuffer(queue, result_buf, CL_TRUE, 0, sizeof(u64), &zero, 0, NULL, NULL);
        err = clEnqueueNDRangeKernel(queue, kernel, 1, NULL, &work_items, NULL, 0, NULL, NULL);
        clFinish(queue);
        // 最後のループのみ結果読み取り
        if (iter == iterations - 1) {
            clEnqueueReadBuffer(queue, result_buf, CL_TRUE, 0, sizeof(u64), &found, 0, NULL, NULL);
        }
    }
    clock_gettime(CLOCK_MONOTONIC, &t1);
    double sec = (t1.tv_sec - t0.tv_sec) + (t1.tv_nsec - t0.tv_nsec) / 1e9;
    printf("完了。見つかった因数: %llu\n", found);
    printf("GPU並列総実行時間: %.6f 秒 (1イテレーション平均: %.6f 秒)\n\n", sec, sec / iterations);
    clReleaseMemObject(n_buf);
    clReleaseMemObject(base_buf);
    clReleaseMemObject(result_buf);
    clReleaseKernel(kernel);
    clReleaseProgram(program);
    clReleaseCommandQueue(queue);
    clReleaseContext(context);
    return 0;
}

/*
 * rsa_manual_opencl.c
 * OpenCL版：多数のGPUスレッド（ワークアイテム）に候補除数を1つずつ割り当て、
 * 全候補を並列でチェックする試し割り法。
 *
 * コンパイル (Windows / MinGW-w64, GCC。OpenCL SDKのヘッダ・ライブラリが必要):
 *   gcc -O2 -I"<OpenCL_include>" -L"<OpenCL_lib>" -o rsa_manual_opencl.exe rsa_manual_opencl.c -lOpenCL
 *
 * または、カレントディレクトリにCLヘッダーを配置し、システムのDLLを直接リンクする場合:
 *   gcc -O2 -I. -o rsa_manual_opencl.exe rsa_manual_opencl.c C:\Windows\System32\OpenCL.dll
 *
 * コンパイル (Linux):
 *   gcc -O2 -o rsa_manual_opencl rsa_manual_opencl.c -lOpenCL
 *
 * 実行にはGPU（またはCPU上で動くOpenCLランタイム）が必要です。
 * 手元にNVIDIA/AMD/Intelいずれかの環境があれば、そのまま動作します。
 */
#define CL_TARGET_OPENCL_VERSION 120
#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <string.h>
#include <time.h>
#include <CL/cl.h>
typedef unsigned long long u64;
/* GPU側カーネル：各ワークアイテムが1つの候補除数dを担当し、n % d == 0 かを判定する
 * メモリ使用量を削減するため、発見した因数を結果バッファのアドレスに直接書き込む */
static const char *kernel_src =
"__kernel void trial_division(__global const ulong *n_ptr,\n"
"                              __global const ulong *base_ptr,\n"
"                              __global ulong *result) {\n"
"    ulong n = *n_ptr;\n"
"    ulong base = *base_ptr;\n"
"    size_t gid = get_global_id(0);\n"
"    ulong d = base + gid * 2; /* 奇数候補のみを走査 */\n"
"    if (d >= 3 && n % d == 0) {\n"
"        *result = d;\n;\n"
"    }\n"
"}\n";
int main(void) {
    cl_int err;
    cl_uint num_platforms = 0;
    clGetPlatformIDs(0, NULL, &num_platforms);
    if (num_platforms == 0) {
        fprintf(stderr, "OpenCLプラットフォームが見つかりません。"
                        "GPUドライバ（またはCPU用OpenCLランタイム）を導入してください。\n");
        return 1;
    }
    cl_platform_id platform;
    clGetPlatformIDs(1, &platform, NULL);
    cl_device_id device;
    err = clGetDeviceIDs(platform, CL_DEVICE_TYPE_GPU, 1, &device, NULL);
    if (err != CL_SUCCESS) {
        /* GPUが無ければCPUデバイスにフォールバック（動作確認用） */
        err = clGetDeviceIDs(platform, CL_DEVICE_TYPE_CPU, 1, &device, NULL);
        if (err != CL_SUCCESS) {
            fprintf(stderr, "使用可能なOpenCLデバイスが見つかりません。\n");
            return 1;
        }
        fprintf(stderr, "[注意] GPUが見つからないため、OpenCL経由のCPUデバイスで実行します。\n");
    }
    cl_context context = clCreateContext(NULL, 1, &device, NULL, NULL, &err);
    cl_command_queue queue = clCreateCommandQueue(context, device, 0, &err);
    cl_program program = clCreateProgramWithSource(context, 1, &kernel_src, NULL, &err);
    err = clBuildProgram(program, 1, &device, NULL, NULL, NULL);
    if (err != CL_SUCCESS) {
        char log[4096];
        clGetProgramBuildInfo(program, device, CL_PROGRAM_BUILD_LOG, sizeof(log), log, NULL);
        fprintf(stderr, "カーネルビルドエラー:\n%s\n", log);
        return 1;
    }
    cl_kernel kernel = clCreateKernel(program, "trial_division", &err);
    /* 対象n（4294967291 * 4294967291 = 18446744030759878681） */
    u64 n = 4294967291ULL * 4294967291ULL;
    u64 base = 3;
    size_t work_items = (size_t)(4294967291ULL / 2); /* 奇数候補数 = 約21億4700万 */
    int iterations = 2;
    cl_mem n_buf = clCreateBuffer(context, CL_MEM_READ_ONLY | CL_MEM_COPY_HOST_PTR, sizeof(u64), &n, &err);
    cl_mem base_buf = clCreateBuffer(context, CL_MEM_READ_ONLY | CL_MEM_COPY_HOST_PTR, sizeof(u64), &base, &err);
    cl_mem result_buf = clCreateBuffer(context, CL_MEM_READ_WRITE, sizeof(u64), NULL, &err);
    clSetKernelArg(kernel, 0, sizeof(cl_mem), &n_buf);
    clSetKernelArg(kernel, 1, sizeof(cl_mem), &base_buf);
    clSetKernelArg(kernel, 2, sizeof(cl_mem), &result_buf);
    printf("=== GPU (OpenCL) 10秒負荷（巨大合成数・2イテレーション） ===\n");
    printf("n = %llu (factor = 4294967291)\n", n);
    printf("ワークアイテム数: %zu\n", work_items);
    printf("実行中...\n");
    struct timespec t0, t1;
    clock_gettime(CLOCK_MONOTONIC, &t0);
    u64 found = 0;
    for (int iter = 0; iter < iterations; iter++) {
        // 結果バッファを 0 にリセット
        u64 zero = 0;
        clEnqueueWriteBuffer(queue, result_buf, CL_TRUE, 0, sizeof(u64), &zero, 0, NULL, NULL);
        err = clEnqueueNDRangeKernel(queue, kernel, 1, NULL, &work_items, NULL, 0, NULL, NULL);
        clFinish(queue);
        // 最後のループのみ結果読み取り
        if (iter == iterations - 1) {
            clEnqueueReadBuffer(queue, result_buf, CL_TRUE, 0, sizeof(u64), &found, 0, NULL, NULL);
        }
    }
    clock_gettime(CLOCK_MONOTONIC, &t1);
    double sec = (t1.tv_sec - t0.tv_sec) + (t1.tv_nsec - t0.tv_nsec) / 1e9;
    printf("完了。見つかった因数: %llu\n", found);
    printf("GPU並列総実行時間: %.6f 秒 (1イテレーション平均: %.6f 秒)\n\n", sec, sec / iterations);
    clReleaseMemObject(n_buf);
    clReleaseMemObject(base_buf);
    clReleaseMemObject(result_buf);
    clReleaseKernel(kernel);
    clReleaseProgram(program);
    clReleaseCommandQueue(queue);
    clReleaseContext(context);
    return 0;
}

"""
shor_period_finding_simulation.py
Shorのアルゴリズムの核心（量子フーリエ変換による周期発見）を、
古典コンピュータ上でシミュレーションし、「量子重ね合わせ→測定」で
なぜ一撃で周期rが見えるのかを体感するためのスクリプト。
これは量子コンピュータの代わりに、小さな系（N=15）に限って
状態ベクトルを配列として愚直に計算しているだけであり、
本物の量子コンピュータの指数的な資源優位性そのものを再現するものではない
（配列サイズが2^n bitで指数的に増えるため、大きなNではこの古典シミュレーション
自体が破綻する。それこそが量子コンピュータが必要な理由）。
実行:
    python3 shor_period_finding_simulation.py
Windows環境での実行（日本語文字化け防止）:
    $env:PYTHONIOENCODING="utf-8"; python shor_period_finding_simulation.py
    または
    set PYTHONIOENCODING=utf-8 && python shor_period_finding_simulation.py
"""
import numpy as np
def f(a, N, x):
    return pow(a, x, N)
def simulate_period_finding(a, N, Q):
    """
    Q個の位相状態にわたって f(x) = a^x mod N を計算し、
    「周期r ごとに同じ値が現れる」という構造を、
    離散フーリエ変換（QFTの古典版）で可視化する。
    """
    xs = list(range(Q))
    values = np.array([f(a, N, x) for x in xs])
    # f(x)が周期r（余りが最初に1に戻る周期）を持つことを直接確認
    r = None
    for candidate_r in range(1, Q):
        if values[candidate_r] == values[0]:
            r = candidate_r
            break
    # 古典DFT（量子コンピュータならQFTを1回の演算で実行できる部分）
    # values自体ではなく、「f(x)の値が特定の値になっているか」を
    # 示す指標関数（例えばf(x)==1の位置）にDFTをかけると、
    # 周期rに対応する周波数にピークが立つ。
    indicator = (values == values[0]).astype(float)
    spectrum = np.abs(np.fft.fft(indicator))
    peak_freqs = np.argsort(spectrum)[::-1][:6]
    return values, r, spectrum, sorted(peak_freqs.tolist())
if __name__ == "__main__":
    a, N, Q = 7, 15, 32
    values, r, spectrum, peaks = simulate_period_finding(a, N, Q)
    print(f"N = {N} を分解したい。ランダムに選んだ底 a = {a}")
    print(f"f(x) = {a}^x mod {N} を x=0..{Q-1} まで計算:")
    print(values.tolist())
    print(f"\n=> 周期 r = {r} （f(x+r) = f(x) となる最小の正整数）")
    print(f"\n古典DFTスペクトルのピーク周波数（上位6件, Q={Q}）: {peaks}")
    print(f"Q/r = {Q}/{r} = {Q/r:.3f}")
    print("=> ピークがQ/rの整数倍の位置に立っていることが確認できる。")
    print("   本物の量子コンピュータでは、この「ピークの位置」を")
    print("   重ね合わせ状態への1回のQFTと測定だけで得られる。")
    print("   古典コンピュータは、これをxごとに逐次計算してからでないと")
    print("   スペクトルを作れない（=全体でO(Q)の逐次計算が必要）。")
    # 古典後処理：周期rからNの素因数を復元する（Shorのアルゴリズムの仕上げ）
    if r % 2 == 0:
        x = pow(a, r // 2, N)
        import math
        p = math.gcd(x - 1, N)
        q = math.gcd(x + 1, N)
        print(f"\n[古典後処理] a^(r/2) mod N = {x}")
        print(f"gcd({x}-1, {N}) = {p}")
        print(f"gcd({x}+1, {N}) = {q}")
        print(f"=> N = {N} = {p} x {q}")

"""
shor_period_finding_simulation.py
Shorのアルゴリズムの核心（量子フーリエ変換による周期発見）を、
古典コンピュータ上でシミュレーションし、「量子重ね合わせ→測定」で
なぜ一撃で周期rが見えるのかを体感するためのスクリプト。
これは量子コンピュータの代わりに、小さな系（N=15）に限って
状態ベクトルを配列として愚直に計算しているだけであり、
本物の量子コンピュータの指数的な資源優位性そのものを再現するものではない
（配列サイズが2^n bitで指数的に増えるため、大きなNではこの古典シミュレーション
自体が破綻する。それこそが量子コンピュータが必要な理由）。
実行:
    python3 shor_period_finding_simulation.py
Windows環境での実行（日本語文字化け防止）:
    $env:PYTHONIOENCODING="utf-8"; python shor_period_finding_simulation.py
    または
    set PYTHONIOENCODING=utf-8 && python shor_period_finding_simulation.py
"""
import numpy as np
def f(a, N, x):
    return pow(a, x, N)
def simulate_period_finding(a, N, Q):
    """
    Q個の位相状態にわたって f(x) = a^x mod N を計算し、
    「周期r ごとに同じ値が現れる」という構造を、
    離散フーリエ変換（QFTの古典版）で可視化する。
    """
    xs = list(range(Q))
    values = np.array([f(a, N, x) for x in xs])
    # f(x)が周期r（余りが最初に1に戻る周期）を持つことを直接確認
    r = None
    for candidate_r in range(1, Q):
        if values[candidate_r] == values[0]:
            r = candidate_r
            break
    # 古典DFT（量子コンピュータならQFTを1回の演算で実行できる部分）
    # values自体ではなく、「f(x)の値が特定の値になっているか」を
    # 示す指標関数（例えばf(x)==1の位置）にDFTをかけると、
    # 周期rに対応する周波数にピークが立つ。
    indicator = (values == values[0]).astype(float)
    spectrum = np.abs(np.fft.fft(indicator))
    peak_freqs = np.argsort(spectrum)[::-1][:6]
    return values, r, spectrum, sorted(peak_freqs.tolist())
if __name__ == "__main__":
    a, N, Q = 7, 15, 32
    values, r, spectrum, peaks = simulate_period_finding(a, N, Q)
    print(f"N = {N} を分解したい。ランダムに選んだ底 a = {a}")
    print(f"f(x) = {a}^x mod {N} を x=0..{Q-1} まで計算:")
    print(values.tolist())
    print(f"\n=> 周期 r = {r} （f(x+r) = f(x) となる最小の正整数）")
    print(f"\n古典DFTスペクトルのピーク周波数（上位6件, Q={Q}）: {peaks}")
    print(f"Q/r = {Q}/{r} = {Q/r:.3f}")
    print("=> ピークがQ/rの整数倍の位置に立っていることが確認できる。")
    print("   本物の量子コンピュータでは、この「ピークの位置」を")
    print("   重ね合わせ状態への1回のQFTと測定だけで得られる。")
    print("   古典コンピュータは、これをxごとに逐次計算してからでないと")
    print("   スペクトルを作れない（=全体でO(Q)の逐次計算が必要）。")
    # 古典後処理：周期rからNの素因数を復元する（Shorのアルゴリズムの仕上げ）
    if r % 2 == 0:
        x = pow(a, r // 2, N)
        import math
        p = math.gcd(x - 1, N)
        q = math.gcd(x + 1, N)
        print(f"\n[古典後処理] a^(r/2) mod N = {x}")
        print(f"gcd({x}-1, {N}) = {p}")
        print(f"gcd({x}+1, {N}) = {q}")
        print(f"=> N = {N} = {p} x {q}")