C++ Baillie-PSW Primality Test: การทดสอบเฉพาะจำนวนแบบ Probabilistic

อัลกอริทึม (High Level)

1. Trial Division ด้วย Primes เล็กๆ (Optional Optimization)
2. Miller-Rabin Test Base 2 (Strong Pseudoprime Test)
3. Lucas Test (หา D, P, Q ที่ Jacobi Symbol (D/n) = -1)
4. ทั้งคู่ Pass → "Probably Prime"

โค้ดตัวอย่าง: Miller-Rabin ส่วนประกอบ (Simple Version)

#include <iostream>
#include <cstdint>
#include <random>

// Modular Multiplication (หลีกเลี่ยง Overflow)
uint64_t mul_mod(uint64_t a, uint64_t b, uint64_t mod) {
    return (__int128)a * b % mod; // GCC/Clang Extension
    // หรือใช้ Russian Peasant Algorithm สำหรับ Portable
}

// Modular Exponentiation
uint64_t pow_mod(uint64_t base, uint64_t exp, uint64_t mod) {
    uint64_t result = 1;
    while (exp > 0) {
        if (exp & 1) result = mul_mod(result, base, mod);
        base = mul_mod(base, base, mod);
        exp >>= 1;
    }
    return result;
}

// Miller-Rabin Single Base
bool millerRabin(uint64_t n, uint64_t a) {
    if (n < 2) return false;
    if (n % 2 == 0) return n == 2;
    
    uint64_t d = n - 1;
    int s = 0;
    while ((d & 1) == 0) { d >>= 1; ++s; }
    
    uint64_t x = pow_mod(a, d, n);
    if (x == 1 || x == n - 1) return true;
    
    for (int r = 1; r < s; ++r) {
        x = mul_mod(x, x, n);
        if (x == n - 1) return true;
    }
    return false;
}

// Baillie-PSW (Simplified: Miller-Rabin Base 2 + Lucas - Full Implementation ยาว)
bool isPrimeBailliePSW(uint64_t n) {
    // 1. Small Prime Trial Division
    static const uint32_t smallPrimes[] = {2,3,5,7,11,13,17,19,23,29,31,37};
    for (uint32_t p : smallPrimes) {
        if (n == p) return true;
        if (n % p == 0) return false;
    }
    
    // 2. Miller-Rabin Base 2
    if (!millerRabin(n, 2)) return false;
    
    // 3. Lucas Test (Omitted for Brevity - ต้องการ Jacobi Symbol, Lucas Sequence)
    // Full Implementation ต้องการ ~100 Lines
    
    return true; // Placeholder
}

int main() {
    uint64_t nums[] = {23, 45, 1000003, 18446744073709551557ULL}; // Large Prime
    for (auto n : nums) {
        std::cout << n << ": " << (isPrimeBailliePSW(n) ? "Probable Prime" : "Composite") << '\n';
    }
}

หมายเหตุสำคัญ

• โค้ดด้านบนเป็น Simplified Miller-Rabin ไม่ใช่ Full Baillie-PSW (ข้าม Lucas Test)
• Full Baillie-PSW Implementation ซับซ้อน (ต้องหา Lucas Parameters, Jacobi Symbol)
• สำหรับ Production ใช้ Library: Boost.Multiprecision::miller_rabin_test, GMP, OpenSSL
• Deterministic Miller-Rabin สำหรับ 64-bit: Test Bases {2, 325, 9375, 28178, 450775, 9780504, 1795265022} แล้วได้คำตอบแน่ชัด

Deterministic Miller-Rabin for 64-bit (แนะนำแทน Baillie-PSW สำหรับ 64-bit)

bool isPrime64(uint64_t n) {
    if (n < 2) return false;
    for (uint64_t p : {2,3,5,7,11,13,17,19,23,29,31,37}) {
        if (n % p == 0) return n == p;
    }
    // Deterministic Bases for 64-bit
    static const uint64_t bases[] = {2, 325, 9375, 28178, 450775, 9780504, 1795265022};
    for (uint64_t a : bases) {
        if (!millerRabin(n, a)) return false;
    }
    return true;
}

สรุป

Baillie-PSW เป็น Strong Probabilistic Test แต่ Implementation ซับซ้อน สำหรับ 64-bit Integers ใช้ Deterministic Miller-Rabin (7 Bases) ให้คำตอบแน่ชัด เร็วกว่า และง่ายกว่า สำหรับ Big Integers (>64-bit) ใช้ Boost.Multiprecision + Miller-Rabin Multiple Bases