Primality test আর factoring — Miller–Rabin থেকে কঠিনতার ধারণা

১ · দুটি আলাদা প্রশ্ন — “মৌলিক কি না” আর “উৎপাদক কী”

বাজারঘরের প্রকৌশলীরা দুটি কারণে প্রায়ই বড় মৌলিক সংখ্যা খোঁজেন। প্রথমত, গ্রাহকের id-কে বালতিতে ভাগ করার hash function h(k) = k mod m-এ m মৌলিক নিলে ছাঁচে-ঢালা id-ও ভালোভাবে ছড়ায় (পর্ব ২৬ “সংখ্যাতত্ত্বের প্রয়োগ — check digit, hashing আর pseudorandom সংখ্যা”-এ এর বিশ্লেষণ)। প্রথম খসড়ায় তাঁরা চান 10⁶-এর ঠিক পরের মৌলিক। দ্বিতীয়ত, পেমেন্টের নিরাপত্তায় যে public-key ব্যবস্থা লাগে, তার চাবি বানাতে শত শত অঙ্কের দুটি মৌলিক সংখ্যা লাগে (গণিত পর্ব ২৮ “RSA-এর গণিত — key, encryption আর কেন কাজ করে”-এ)। দুই ক্ষেত্রেই কাজটি একই: একটি বড় বিজোড় সংখ্যা নেওয়া, মৌলিক কি না দেখা, না হলে পরেরটি।

এখানে দুটি প্রশ্ন আলাদা করা জরুরি। Primality: n কি মৌলিক? উত্তর শুধু “হ্যাঁ” বা “না”। Factoring: n যৌগিক হলে তার উৎপাদক কী? দ্বিতীয়টির উত্তর জানলে প্রথমটি জানা যায়, কিন্তু উল্টোটা নয়। এই পর্বের বড় শিক্ষা হলো, প্রথম প্রশ্নটি একটিও উৎপাদক না খুঁজে নিশ্চিতভাবে “যৌগিক” বলতে পারে। পর্ব ২৪-এ এর প্রথম ঝলক দেখেছেন: 2¹⁴ mod 15 = 4 ≠ 1, তাই 15 যৌগিক। একটি primality test ঠিক এমন একটি পদ্ধতি; আর integer factorization মানে সংখ্যাটিকে মৌলিক উৎপাদকে ভেঙে দেখানো।

সবচেয়ে সরল পথ: trial division আর তার খরচ

পর্ব ২০ “বিভাজ্যতা আর মৌলিক সংখ্যা — division algorithm আর unique factorization”-এ প্রমাণ করেছেন: যৌগিক n-এর অন্তত একটি মৌলিক উৎপাদক ≤ √n. তাই 2, 3, 5, … দিয়ে √n পর্যন্ত ভাগ করে কোনো উৎপাদক না পেলে n মৌলিক। এই trial division নির্ভুল, আর ছোট উৎপাদক থাকলে দ্রুত। নিচের stepper বাজারঘরের একটি অর্ডার id 1147-এ এটি চালায়।

1n = 1147, √n ≈ 33.872ভাজক: 2, 3, 5, …, 3131147 mod 2, 3, 5 = 1, 1, 24mod 7, 11, 13, 17 = 6, 3, 3, 85mod 19, 23, 29 = 7, 20, 166mod 31 = 0 ⇒ 1147 = 31 · 377d অঙ্ক ⇒ প্রায় 10^(d/2) ভাগ

এই খরচকে সঠিক ভাষায় বলা দরকার। একটি algorithm-এর খরচ মাপা হয় ইনপুটের দৈর্ঘ্য দিয়ে, আর n-এর দৈর্ঘ্য তার অঙ্ক বা bit-এর সংখ্যা, প্রায় log n. √n = 2^{(log₂ n)/2}, তাই trial division-এর সময় দৈর্ঘ্যের তুলনায় exponential. বিপরীতে pow(a, b, n) (binary exponentiation, পর্ব ২২ আর ১৩) মাত্র প্রায় 2 log₂ b-টি গুণ করে, যা দৈর্ঘ্যের তুলনায় polynomial. এই পর্বের সব দ্রুত পরীক্ষা এই দ্বিতীয় ধরনের।

নিচের ব্লক ৬, ৮, ১০ আর ১২ অঙ্কের সবচেয়ে বড় মৌলিক সংখ্যায় trial division চালিয়ে ভাগের সংখ্যা আর সময় মাপে, তারপর সেই হার ধরে বড় সংখ্যার খরচ অনুমান করে। প্রথমবার চালালে ব্রাউজারে Python (Pyodide) নামতে কয়েক সেকেন্ড লাগে।

import time

def trial(n):                     # smallest factor of n, and how many divisions
    if n % 2 == 0:
        return 2, 1
    d, tries = 3, 1
    while d * d <= n:             # only up to sqrt(n): see lesson 20
        tries += 1
        if n % d == 0:
            return d, tries
        d += 2
    return n, tries               # no factor up to sqrt(n): n is prime

per_div = []
for n in (999983, 99999989, 9999999967, 999999999989):   # largest primes
    t = time.perf_counter()                               # with 6..12 digits
    f, tries = trial(n)
    sec = time.perf_counter() - t
    per_div.append(sec / tries)
    print(f"{len(str(n)):2d} digits  prime={f == n}  divisions={tries:>7,}"
          f"  time={sec * 1000:7.1f} ms")

c = per_div[-1]                   # seconds per division (rough)
for digits in (20, 30, 40):
    sec = 10 ** (digits / 2) / 2 * c          # sqrt(n)/2 odd divisions
    print(f"{digits} digits, worst case: about {sec:.0e} s"
          f" = {sec / 31557600:.0e} years")

for n in (1147, 2026111001):      # bazarghar order ids
    print(n, "-> smallest factor, divisions:", trial(n))

 6 digits  prime=True  divisions=    500  time=    0.1 ms
 8 digits  prime=True  divisions=  5,000  time=    0.5 ms
10 digits  prime=True  divisions= 50,000  time=    7.7 ms
12 digits  prime=True  divisions=500,000  time=   67.8 ms
20 digits, worst case: about 7e+02 s = 2e-05 years
30 digits, worst case: about 7e+07 s = 2e+00 years
40 digits, worst case: about 7e+12 s = 2e+05 years
1147 -> smallest factor, divisions: (31, 16)
2026111001 -> smallest factor, divisions: (191, 96)

২ · Fermat-এর পরীক্ষা কোথায় আটকায়, আর 1-এর বর্গমূল

পর্ব ২৪ “Fermat-এর ছোট উপপাদ্য আর Euler-এর totient”-এর Fermat-এর পরীক্ষা এক লাইনে: aⁿ⁻¹ mod n ≠ 1 হলে n যৌগিক। এটি দ্রুত, কিন্তু উল্টো দিকে খাটে না। সেখানে দেখেছেন, 341 base 2-এ পাস করে, আর 561 = 3 · 11 · 17 একটি Carmichael সংখ্যা: প্রতিটি সহমৌলিক base a-তে a⁵⁶⁰ ≡ 1 (mod 561). Carmichael সংখ্যা অসীম সংখ্যক, তাই base বাড়িয়ে Fermat-এর পরীক্ষাকে নির্ভরযোগ্য করা যায় না। সহমৌলিক নয় এমন base খুঁজে পাওয়া মানে আসলে একটি উৎপাদক খুঁজে পাওয়া, যা বড় সংখ্যায় ঠিক সেই কঠিন কাজ।

তাহলে aⁿ⁻¹ ≡ 1-এর চেয়ে বেশি তথ্য কোথায় পাওয়া যায়? উত্তর: সেই 1-এ পৌঁছানোর পথে। n − 1 জোড় হলে aⁿ⁻¹ হলো a^(n−1)/2-এর বর্গ। তাই a^(n−1)/2 হলো 1-এর একটি বর্গমূল mod n. মৌলিক modulus-এ 1-এর বর্গমূল মাত্র দুটি।

প্রমাণ:

1. (⇐) x ≡ ±1 হলে x² ≡ (±1)² = 1. — গুণ congruence মানে, পর্ব ২২
2. (⇒) x² ≡ 1 মানে p ∣ x² − 1 = (x − 1)(x + 1). — congruence-এর সংজ্ঞা
3. p মৌলিক, তাই Euclid-এর lemma অনুযায়ী p ∣ (x − 1) বা p ∣ (x + 1), অর্থাৎ x ≡ 1 বা x ≡ −1 (mod p). ∎ — Euclid-এর lemma, পর্ব ২১

ধাপ ৩-এ মৌলিকতা অপরিহার্য। যৌগিক modulus-এ 1-এর আরও বর্গমূল থাকতে পারে: mod 15-এ 4² = 16 ≡ 1, অথচ 4 ≢ ±1 (−1 ≡ 14). কারণ 15 ∣ (4 − 1)(4 + 1) = 3 · 5, কিন্তু 15 কোনো একটি উৎপাদককে আলাদাভাবে ভাগ করে না। এমন x-কে বলে 1-এর অদ্ভুত বর্গমূল (nontrivial square root of 1). আর এর একটি দারুণ পার্শ্বফল আছে: n ∣ (x − 1)(x + 1) অথচ n কোনোটিকে ভাগ করে না, তাই gcd(x − 1, n) হলো n-এর একটি প্রকৃত উৎপাদক। mod 15-এ gcd(3, 15) = 3, gcd(5, 15) = 5. অর্থাৎ অদ্ভুত বর্গমূল খুঁজে পাওয়া মানে শুধু “যৌগিক” জানা নয়, উৎপাদকও হাতে পাওয়া।

৩ · Miller–Rabin পরীক্ষা — বিবৃতি, প্রমাণ আর 561

বিজোড় n > 2-এ n − 1 জোড়। তাকে যতবার পারা যায় 2 দিয়ে ভাগ করুন: n − 1 = 2^s · d, যেখানে s ≥ 1 আর d বিজোড়। যেমন 560 = 2⁴ · 35, 96 = 2⁵ · 3. এবার a^d থেকে শুরু করে s বার বর্গ করলে aⁿ⁻¹ পাওয়া যায়:

a^d, a^2d, a^4d, …, a^2s−1d, a^2sd = aⁿ⁻¹ (সব mod n)

প্রমাণ: x_r = a^2rd mod p লিখুন, r = 0, 1, …, s. তাহলে x_r+1 ≡ x_r².

1. x_s = a^p−1 ≡ 1 (mod p). — Fermat-এর ছোট উপপাদ্য, p মৌলিক আর p ∤ a
2. x₀ ≡ 1 হলে প্রথম বিকল্পটি সত্য, প্রমাণ শেষ। — case ১
3. নাহলে ক্রমে অন্তত একটি পদ (x₀) 1 নয়, আর শেষ পদ 1. r নিন সবচেয়ে বড় সূচক যেখানে x_r ≢ 1; তাহলে r ≤ s − 1 আর x_r+1 ≡ 1. — well-ordering, সসীম সেটে সবচেয়ে বড় সদস্য
4. x_r² ≡ x_r+1 ≡ 1, তাই Lemma ২.১-এ x_r ≡ ±1. x_r ≢ 1 বলে x_r ≡ −1, অর্থাৎ দ্বিতীয় বিকল্পটি সত্য। ∎ — Lemma ২.১

উপপাদ্যটিকে contrapositive আকারে পড়লেই পরীক্ষা। বিজোড় n আর 1 < a < n − 1 নিন। যদি a^d ≢ 1 আর প্রতিটি 0 ≤ r ≤ s − 1-এ a^2rd ≢ −1 (mod n), তবে n নিশ্চিতভাবে যৌগিক।

লক্ষ করুন, strong probable prime সবসময় Fermat-এর পরীক্ষাতেও পাস করে: শর্তের দুই ক্ষেত্রেই বর্গ করতে করতে শেষে aⁿ⁻¹ ≡ 1 আসে। তাই নতুন পরীক্ষা পুরোনোটির চেয়ে কখনো দুর্বল নয়। আর এটি কতটা শক্তিশালী, 561-এ দেখুন। Fermat base 2-এ 561 পাস করে, কিন্তু চেইনের মাঝখানে কী হয়?

1n = 561, a = 22560 = 2⁴ · 35 (s = 4, d = 35)3x₀ = 2³⁵ mod 561 = 2634x₁ = 263² mod 561 = 1665x₂ = 166² mod 561 = 676x₃ = 67² mod 561 = 1767 ≢ ±1 ⇒ 561 যৌগিক8gcd(66, 561) = 33, gcd(68, 561) = 17

এবার পরীক্ষাটি নিজে লিখুন। নিচের strong_test পুরো চেইন ছাপে, অদ্ভুত বর্গমূল পেলে gcd দিয়ে উৎপাদক দেখায়, আর miller_rabin নির্দিষ্ট seed থেকে ২০টি এলোমেলো base নেয়।

from math import gcd
import random

def strong_test(n, a, show=False):    # n odd, n > 3, 1 < a < n - 1
    s, d = 0, n - 1
    while d % 2 == 0:                  # n - 1 = 2^s * d, d odd
        s, d = s + 1, d // 2
    chain = [pow(a, d, n)]
    for _ in range(s):                 # a^d, a^2d, a^4d, ..., a^(n-1)
        chain.append(chain[-1] ** 2 % n)
    passed = chain[0] == 1 or (n - 1) in chain[:s]
    if show:
        print(f"n={n}, a={a}: n-1 = 2^{s}*{d}, chain {chain}")
        for x, y in zip(chain, chain[1:]):
            if y == 1 and x not in (1, n - 1):        # x^2 = 1, x != +-1
                print(f"   {x}^2 = 1 (mod n): gcd({x}-1, n) =", gcd(x - 1, n))
                break
        print("   ->", "probable prime" if passed else "COMPOSITE, a is a witness")
    return passed

def miller_rabin(n, k=20, rng=random.Random(25)):
    if n < 4:
        return n in (2, 3)
    if n % 2 == 0:
        return False
    return all(strong_test(n, rng.randrange(2, n - 1)) for _ in range(k))

strong_test(561, 2, show=True)        # Carmichael: Fermat says 1 for a = 2
strong_test(2047, 2, show=True)       # 23 * 89 fools base 2 ...
strong_test(2047, 3, show=True)       # ... but not base 3
print("Fermat base 2 on 561:", pow(2, 560, 561))
for n in (561, 1105, 1729, 2**61 - 1, 2**67 - 1, 10**18 + 9):
    print(n, "->", "prime" if miller_rabin(n) else "composite")

n=561, a=2: n-1 = 2^4*35, chain [263, 166, 67, 1, 1]
   67^2 = 1 (mod n): gcd(67-1, n) = 33
   -> COMPOSITE, a is a witness
n=2047, a=2: n-1 = 2^1*1023, chain [1, 1]
   -> probable prime
n=2047, a=3: n-1 = 2^1*1023, chain [1565, 1013]
   -> COMPOSITE, a is a witness
Fermat base 2 on 561: 1
561 -> composite
1105 -> composite
1729 -> composite
2305843009213693951 -> prime
147573952589676412927 -> composite
1000000000000000009 -> prime

৪ · ভুলের সম্ভাবনা, নির্দিষ্ট base আর AKS

Miller–Rabin-এর “যৌগিক” উত্তর সবসময় সঠিক, কারণ তার পেছনে একটি witness আছে। ভুল হতে পারে শুধু উল্টো দিকে: যৌগিক n-এর জন্য এমন base বাছা, যে strong liar. তাহলে প্রশ্ন, একটি যৌগিক সংখ্যার কতগুলো strong liar থাকতে পারে?

প্রমাণটি এখানে বাদ দেওয়া হলো: এতে mod n-এর সহমৌলিক ভাগশেষের group-গঠন আর CRT-এর সূক্ষ্ম গণনা লাগে (group-এর ভাষা পর্ব ৭০ “Group — symmetry, cyclic group আর Lagrange-এর উপপাদ্য”-এ আসবে)। তবে মূল ধারণাটি অনুশীলন ১-এ নিজে প্রমাণ করবেন: n-এর দুটি ভিন্ন বিজোড় মৌলিক উৎপাদক থাকলে 1-এর অদ্ভুত বর্গমূল থাকতেই হবে, আর বেশিরভাগ base চেইনের কোথাও সেটিকে বা অন্য কোনো অসংগতিকে সামনে আনে। একই সীমায় 9 বাদ, কারণ সেখানে 8টি base-এর মধ্যে 2টি (1 আর 8) liar, আর φ(9)/4 = 1.5.

সীমাটি পড়ার সময় একটি সূক্ষ্মতা মনে রাখুন। 4^−k হলো একটি নির্দিষ্ট যৌগিক n-এর জন্য, base বাছাইয়ের এলোমেলোতার ওপর সম্ভাবনা। এটি “পাস করা সংখ্যাটি যৌগিক হওয়ার সম্ভাবনা” নয়; সেটি নির্ভর করে কোন সংখ্যা থেকে n এসেছে তার ওপর (Bayes-এর উপপাদ্য, পর্ব ৩৯)। বাস্তবে এলোমেলো বড় বিজোড় সংখ্যায় liar-এর অনুপাত সীমার চেয়ে অনেক কম, তাই কয়েকটি base-ই যথেষ্ট। k = 20-এ সীমা 4⁻²⁰ ≈ 9 · 10⁻¹³. এমন algorithm, যা নির্দিষ্ট সময়ে শেষ হয় কিন্তু খুব ছোট সম্ভাবনায় ভুল বলতে পারে, তাকে Monte Carlo algorithm বলে, আর বারবার চালিয়ে ভুল কমানো হলো amplification (DSA কোর্সের পর্ব ৬১).

নিচের ব্লক 800-এর নিচে প্রতিটি বিজোড় যৌগিক n-এ সব base গুনে Fermat liar আর strong liar তুলনা করে, Monier–Rabin-এর অনুপাত দেখে, তারপর 10⁵-এর নিচে কোন যৌগিক সংখ্যা base 2-কে, আর কোনটি একসঙ্গে 2 আর 3-কে ঠকায়, তা খোঁজে।

from math import gcd

def strong_test(n, a):                 # same test as before, no printing
    s, d = 0, n - 1
    while d % 2 == 0:
        s, d = s + 1, d // 2
    x = pow(a, d, n)
    if x in (1, n - 1):
        return True
    for _ in range(s - 1):
        x = x * x % n
        if x == n - 1:
            return True
    return False

N = 10**5
prime = [False, False] + [True] * (N - 2)     # sieve of Eratosthenes
for i in range(2, int(N ** 0.5) + 1):
    if prime[i]:
        prime[i * i::i] = [False] * len(range(i * i, N, i))

worst = (0, None)
for n in (n for n in range(11, 800, 2) if not prime[n]):   # odd composites, n > 9
    strong = sum(strong_test(n, a) for a in range(1, n))
    phi = sum(gcd(a, n) == 1 for a in range(1, n))
    worst = max(worst, (strong / phi, n))
    if n in (91, 341, 561):
        fermat = sum(pow(a, n - 1, n) == 1 for a in range(1, n))
        print(f"n={n}: phi={phi}, Fermat liars {fermat}, strong liars {strong}")
print(f"largest strong liars / phi(n) for 9 < n < 800: {worst[0]:.3f} at n = {worst[1]}")

fool2 = [n for n in range(9, N, 2) if not prime[n] and strong_test(n, 2)]
print("odd composites below 10^5 passing base 2:", len(fool2), fool2[:5], "...")
print("passing bases 2 and 3:", [n for n in fool2 if strong_test(n, 3)])
print("1373653 = 829*1657, bases 2, 3, 5:", [strong_test(1373653, a) for a in (2, 3, 5)])

n=91: phi=72, Fermat liars 36, strong liars 18
n=341: phi=300, Fermat liars 100, strong liars 50
n=561: phi=320, Fermat liars 320, strong liars 10
largest strong liars / phi(n) for 9 < n < 800: 0.250 at n = 703
odd composites below 10^5 passing base 2: 16 [2047, 3277, 4033, 4681, 8321] ...
passing bases 2 and 3: []
1373653 = 829*1657, bases 2, 3, 5: [True, True, False]

এলোমেলো base ছাড়া: নির্দিষ্ট base-এর তালিকা

শেষ লাইনগুলো একটি কৌশলের ইঙ্গিত দেয়। যদি জানা থাকে যে একটি সীমার নিচে কোনো যৌগিক সংখ্যা একসঙ্গে কয়েকটি নির্দিষ্ট base-কে ঠকায় না, তাহলে সেই সীমার নিচে পরীক্ষাটি পুরোপুরি deterministic, কোনো ভুলের সম্ভাবনা নেই। এমন সীমা বের করা হয়েছে বিশাল কম্পিউটার-অনুসন্ধানে। নিচের তালিকার প্রতিটি “সবচেয়ে ছোট যৌগিক” সংখ্যা প্রকাশিত গবেষণা থেকে নেওয়া; এই পর্বের জন্য ডেস্কটপে (§৬) প্রতিটির উৎপাদক, তালিকার সব base-এ পাস করা আর পরের মৌলিক base-এ ধরা পড়া যাচাই করা হয়েছে, আর 2 · 10⁶ পর্যন্ত সম্পূর্ণ খুঁজে দ্বিতীয় সারির দাবিটি মেলানো হয়েছে।

অর্থাৎ n < 3,215,031,751 হলে base 2, 3, 5, 7-এ পাস করা মানেই মৌলিক, আর 64-bit পূর্ণসংখ্যার (n < 2⁶⁴ ≈ 1.8 · 10¹⁹) পুরো পরিসর প্রথম ১২টি মৌলিক base-এ নিশ্চিতভাবে ঢাকা পড়ে। প্রথম চারটি সারির সীমা Carl Pomerance, John Selfridge আর Samuel Wagstaff-এর ১৯৮০ সালের হিসাব বলে পরিচিত; শেষ দুটি Jonathan Sorenson আর Jonathan Webster-এর ২০১৫ সালের দিকের কাজ। সংখ্যাগুলো উদ্ধৃতি; এই পর্বের script সীমার নিচের সব সংখ্যা পরীক্ষা করেনি, শুধু উপরের দাবিগুলো।

AKS — মৌলিকতা polynomial সময়ে

Miller–Rabin দ্রুত কিন্তু randomized; নির্দিষ্ট base-এর তালিকা শুধু একটি সীমা পর্যন্ত। সব n-এর জন্য deterministic আর polynomial-সময়ের পরীক্ষা কি আছে? ২০০২ সালে কানপুরের IIT-এর Manindra Agrawal আর তাঁর দুই ছাত্র Neeraj Kayal ও Nitin Saxena হ্যাঁ বলে প্রমাণ দেন; শিরোনাম ছিল “PRIMES is in P”, আর প্রকাশ ২০০৪ সালে। ভিত্তিটি Fermat-এরই একটি polynomial রূপ: gcd(a, n) = 1 হলে

n মৌলিক ⇔ (x + a)ⁿ ≡ xⁿ + a (mod n), polynomial হিসেবে

(⇒) দিকটি আপনি এখনই দেখতে পারেন: p মৌলিক হলে 0 < k < p-এ binomial coefficient C(p, k) = p!/(k!(p − k)!)-এর লবে p আছে কিন্তু হরে নেই, তাই p ∣ C(p, k); মাঝের সব পদ মুছে যায় আর a^p ≡ a (Fermat). পুরো polynomial যাচাই করতে n-টি সহগ লাগে, যা আবার exponential; AKS-এর আসল কৌশল হলো এটিকে একটি ছোট polynomial x^r − 1-এর modulus-এ কয়েকটি a-তে পরীক্ষা করা। মূল সংস্করণের সময় মোটামুটি (log n)¹²-এর মতো, পরে Hendrik Lenstra আর Carl Pomerance-এর উন্নত রূপে প্রায় (log n)⁶. তত্ত্বে এটি যুগান্তকারী, কিন্তু বাস্তবে Miller–Rabin আর তার আত্মীয়রা অনেক দ্রুত, তাই সফটওয়্যারে প্রায় সবসময় সেগুলোই চলে।

৫ · Factoring কেন কঠিন — Pollard-এর rho

Miller–Rabin বলল 561 যৌগিক, আর ভাগ্যক্রমে উৎপাদকও দিল। কিন্তু সাধারণভাবে Miller–Rabin উৎপাদক দেয় না: 2047 base 3-এ ধরা পড়েছে, অথচ চেইনে কোনো অদ্ভুত বর্গমূল আসেনি। “যৌগিক” জানা আর “কীভাবে ভাঙে” জানা দুটি আলাদা কাজ। Trial division-এর চেয়ে অনেক ভালো কিন্তু তবু সহজ একটি পদ্ধতি হলো John Pollard-এর ১৯৭৫ সালের rho algorithm.

ধারণাটি birthday paradox থেকে (পর্ব ৪২ “Birthday paradox, hash collision আর balls-into-bins”-এ গণিত): p-টি ঘরে এলোমেলো বল ফেললে মাত্র প্রায় √p বলেই দুটি একই ঘরে পড়ার ভালো সম্ভাবনা। ধরুন n-এর একটি অজানা মৌলিক উৎপাদক p. একটি ক্রম চালান x₀ = 2, x_i+1 = f(x_i) = x_i² + 1 mod n. এই ক্রমকে mod p দেখলে (যা আমরা সরাসরি দেখতে পাই না) সেটিও একই নিয়মে চলে, কারণ যোগ-গুণ congruence মানে: x_i+1 mod p = (x_i mod p)² + 1 mod p. মাত্র p-টি সম্ভাব্য মান, তাই pigeonhole নীতিতে কোনো এক সময় মান পুনরাবৃত্ত হয়, আর তারপর থেকে ক্রমটি চক্রে ঘোরে। লেজ আর চক্র মিলে আকারটি গ্রিক অক্ষর ρ (rho)-এর মতো।

ছবিতে x₃ = 677 আর x₆ = 871, mod 8051-এ আলাদা কিন্তু mod 97-এ দুটিই 95. তাই 97 ∣ (871 − 677) = 194, অথচ 8051 ∤ 194. ফলে gcd(194, 8051) = 97, একটি প্রকৃত উৎপাদক, আর আমরা 97 আগে থেকে না জেনেই সেটি পেলাম। বাকি প্রশ্ন: দুই পদ যে mod p-এ মিলেছে, সেটা জানব কীভাবে, যখন p-ই অজানা? উত্তর: প্রতিটি সম্ভাব্য জোড়ার gcd দেখা ব্যয়বহুল, তাই Floyd-এর cycle detection (কচ্ছপ আর খরগোশ) দিয়ে শুধু x_i আর x_2i জোড়া দেখা হয় (DSA কোর্সের পর্ব ১২-এ linked list-এ একই কৌশল)।

প্রমাণ:

1. চক্রে ঢোকার পর ক্রম প্রতি λ ধাপে নিজেকে পুনরাবৃত্তি করে, তাই i, j ≥ μ হলে x_i = x_j ⇔ λ ∣ (j − i). — চক্রের সংজ্ঞা, ভাগশেষ mod λ
2. j = 2i নিলে শর্তটি হয় λ ∣ i. তাই এমন i চাই যা λ-এর গুণিতক আর ≥ μ (এবং ≥ 1). — ধাপ ১
3. λ, 2λ, 3λ, … এর মধ্যে প্রথম যে গুণিতকটি ≥ max(μ, 1), সেটি তার আগের গুণিতক থেকে ঠিক λ বড়, আর আগেরটি μ-এর ছোট (বা শূন্য)। তাই সেটি ≤ μ + λ. ∎ — division algorithm, পর্ব ২০

Lemma-টি mod p-এর ক্রমে লাগালে Pollard-এর algorithm পাওয়া যায়: x = y = 2 থেকে শুরু, প্রতি ধাপে x এক পা (x ← f(x)), y দুই পা (y ← f(f(y))), আর d = gcd(|x − y|, n) দেখা। d = 1 হলে চলতে থাকুন; 1 < d < n হলে একটি উৎপাদক পাওয়া গেল; d = n হলে এই f-এ ব্যর্থ, অন্য c নিয়ে f(x) = x² + c দিয়ে আবার চেষ্টা। mod p-এর চক্রে μ + λ ধাপের মধ্যে মিল হবেই। যদি f এলোমেলো ফাংশনের মতো আচরণ করে, তাহলে birthday-যুক্তিতে μ + λ গড়ে প্রায় √(πp/2) ≈ 1.25√p. এটি একটি heuristic, প্রমাণিত উপপাদ্য নয়: x² + c সত্যিই এলোমেলো নয়, তবে অভিজ্ঞতায় অনুমানটি ভালো খাটে। ফলে সবচেয়ে ছোট মৌলিক উৎপাদক p হলে খরচ প্রায় √p ≤ n^1/4 ধাপ, trial division-এর p ধাপের তুলনায় নাটকীয়ভাবে কম।

1n = 8051, x = y = 22i=1: x = 5, y = 263    gcd(21, 8051) = 14i=2: x = 26, y = 74745    gcd(7448, 8051) = 16i=3: x = 677, y = 8717    gcd(194, 8051) = 9788051 = 97 · 83

এবার একটি বড় পরীক্ষা। নিচের ব্লক ছোট উদাহরণগুলো চালায়, তারপর ১০,০০০তম মৌলিক 104729 আর ১,০০,০০০তম মৌলিক 1299709-এর গুণফল (১২ অঙ্কের একটি semiprime) ভাঙে, ৪০টি এমন সংখ্যায় গড় ধাপ মাপে, আর একটি ব্যর্থতার উদাহরণ দেখায়।

from math import gcd, isqrt

def rho(n, c=1, x0=2):
    f = lambda v: (v * v + c) % n     # x -> x^2 + c (mod n)
    x = y = x0
    steps = 0
    while True:
        x, y = f(x), f(f(y))          # tortoise: 1 step, hare: 2 steps
        steps += 1
        d = gcd(abs(x - y), n)
        if d != 1:
            return d, steps           # d == n means this c failed

for n in (8051, 10403, 2026111001):
    d, steps = rho(n)
    print(f"n={n}: factor {d} = n/{n // d}, steps {steps}")

n = 104729 * 1299709                  # 10,000th prime * 100,000th prime
d, steps = rho(n)
p = min(d, n // d)
print(f"n={n} ({len(str(n))} digits): {d} * {n // d}")
print(f"   rho steps {steps}, sqrt(p) = {isqrt(p)}, trial division ~ {p // 2} divisions")

total = 0                             # average over many 12-digit semiprimes
ps = [p for p in range(100003, 101000, 2) if all(p % q for q in range(3, 320, 2))]
qs = [q for q in range(10**6 + 3, 10**6 + 2000, 2) if all(q % r for r in range(3, 1001, 2))]
pairs = list(zip(ps[:40], qs[:40]))
for p, q in pairs:
    d, steps = rho(p * q)
    total += steps
print(f"average rho steps over {len(pairs)} semiprimes (p ~ 10^5): {total / len(pairs):.0f}",
      f"| sqrt(p) ~ {isqrt(100000)}")

print("rho on 25 = 5*5:", rho(25), "  with c = 2:", rho(25, c=2))

n=8051: factor 97 = n/83, steps 3
n=10403: factor 101 = n/103, steps 9
n=2026111001: factor 191 = n/10607911, steps 15
n=136117223861 (12 digits): 104729 * 1299709
   rho steps 278, sqrt(p) = 323, trial division ~ 52364 divisions
average rho steps over 40 semiprimes (p ~ 10^5): 332 | sqrt(p) ~ 316
rho on 25 = 5*5: (25, 3)   with c = 2: (5, 2)

৬ · ডেস্কটপে সত্যিকারের run — sympy, gmpy2 আর বড় সংখ্যা

ব্রাউজারের ব্লকগুলো ছোট সংখ্যায় সীমিত। বড় সংখ্যায় আমাদের হাতে লেখা পরীক্ষা পেশাদার লাইব্রেরির সঙ্গে কেমন পাল্লা দেয়, তা দেখতে এই কোর্সের লেখার সময় একটি সাধারণ Windows ডেস্কটপে সত্যিকারের run করা হয়েছে: Python 3.12.13, sympy 1.14.0 (BSD লাইসেন্স) আর gmpy2 2.3.2 (GMP লাইব্রেরির Python সংযোগ, LGPL লাইসেন্স), uv দিয়ে বানানো একটি ছোট প্রজেক্টে। নিজে চালাতে চাইলে:

cd ~/bazarghar-dm
uv add sympy gmpy2
uv run python bench.py
uv run python verify.py

bench.py-এর মূল অংশ নিচে। is_sprp (strong test), mr আর rho ফাংশন তিনটি §৩ আর §৫-এর ব্রাউজার-ব্লকের ফাংশনের হুবহু সমতুল্য, তাই এখানে বাদ দেওয়া হলো। প্রতিটি পরীক্ষা ২০০ বার চালিয়ে এক বারের গড় সময় মাপা হয়েছে (us মানে মাইক্রোসেকেন্ড)।

import sys, time, random
import sympy, gmpy2
# is_sprp, mr and rho: same as the in-page blocks

def per_call(fn, arg, reps):
    t = time.perf_counter()
    for _ in range(reps):
        r = fn(arg)
    return r, (time.perf_counter() - t) / reps * 1e6

print("Python", sys.version.split()[0], "| sympy", sympy.__version__, "| gmpy2", gmpy2.version())
p20, p40, p60 = (sympy.nextprime(10**k) for k in (19, 39, 59))
# Chernick form (6k+1)(12k+1)(18k+1): a Carmichael number when all three are prime
k = next(k for k in range(10**19, 10**19 + 10**6)
         if all(sympy.isprime(m * k + 1) for m in (6, 12, 18)))
car = (6*k + 1) * (12*k + 1) * (18*k + 1)
tests = [("prime 20d", p20), ("prime 40d", p40), ("prime 60d", p60),
         ("p40*q20 (60d)", p40 * sympy.nextprime(10**20)),
         ("Carmichael %dd" % len(str(car)), car)]
print("%-16s %-6s %-6s %-6s %10s %10s %10s" % ("number", "MR20", "sympy", "gmpy2", "MR us", "sympy us", "gmpy2 us"))
for name, n in tests:
    a, t1 = per_call(mr, n, 200)
    b, t2 = per_call(sympy.isprime, n, 200)
    c, t3 = per_call(lambda v: gmpy2.is_prime(v, 25), n, 200)
    print("%-16s %-6s %-6s %-6s %10.1f %10.1f %10.1f" % (name, a, b, c, t1, t2, t3))
print("Carmichael: 2^(n-1) mod n == 1 ?", pow(2, car - 1, car) == 1, "| strong base 2 ?", is_sprp(car, 2))

rng = random.Random(2026)
p = sympy.nextprime(rng.randrange(10**10, 10**11))
q = sympy.nextprime(rng.randrange(10**12, 10**13))
n = p * q
t = time.perf_counter(); d, steps = rho(n); t_rho = time.perf_counter() - t
t = time.perf_counter(); fi = sympy.factorint(n); t_f = time.perf_counter() - t
# ... four print lines for these results (see the output)

$ uv run python bench.py
Python 3.12.13 | sympy 1.14.0 | gmpy2 2.3.2
number           MR20   sympy  gmpy2       MR us   sympy us   gmpy2 us
prime 20d        True   True   True        168.5        9.5        4.7
prime 40d        True   True   True        433.7       33.9       20.1
prime 60d        True   True   True        880.6       68.8       51.9
p40*q20 (60d)    False  False  False        47.0        8.7        7.7
Carmichael 61d   False  False  False        55.2        9.7        7.8
Carmichael: 2^(n-1) mod n == 1 ? True | strong base 2 ? False
semiprime n = 448985231080606748777401 (24 digits)
our rho : factor 53461288997 steps 122067 sqrt(p) ~ 231216 time 0.09 s
factorint: {53461288997: 1, 8398324086533: 1} time 0.12 s
trial division would need about 2.7e+10 odd divisions

ফলগুলো পড়ুন। প্রথমত, তিনটি পদ্ধতির উত্তর প্রতিটি সারিতে এক। ৬০ অঙ্কের একটি মৌলিক সংখ্যায় আমাদের ২০-base Miller–Rabin এক মিলিসেকেন্ডেরও কম সময় নেয়; trial division-এ একই কাজে মহাবিশ্বের বয়সের চেয়েও বেশি লাগত। মৌলিক সংখ্যায় সময় বেশি, কারণ ২০টি base-ই পুরো চলে; যৌগিক সংখ্যায় প্রথম base-ই সাধারণত ধরে ফেলে, তাই প্রায় এক base-এর সময়। দ্বিতীয়ত, sympy আর gmpy2 আমাদের চেয়ে অনেক দ্রুত: মৌলিক সংখ্যায় প্রায় ১৩ থেকে ৩৬ গুণ, যৌগিকে প্রায় ৫ থেকে ৭ গুণ। sympy কম রাউন্ড চালায়, আর gmpy2 ২৫ রাউন্ড চালিয়েও দ্রুত, কারণ বড় সংখ্যার গুণ হয় GMP-র অপ্টিমাইজ করা C কোডে। sympy-র নিজের বিবরণ অনুযায়ী 2⁶⁴-এর নিচে এটি নির্দিষ্ট base-এ deterministic Miller–Rabin চালায়, আর তার ওপরে strong BPSW পরীক্ষা, যার কোনো counterexample এখনো জানা নেই। তৃতীয়ত, Chernick-এর সূত্রে বানানো ৬১ অঙ্কের একটি Carmichael সংখ্যা Fermat base 2-এ পাস করে, কিন্তু strong test-এ base 2-এই ধরা পড়ে।

শেষ অংশ factoring. ২৪ অঙ্কের semiprime-এর ছোট উৎপাদক p ≈ 5.3 · 10¹⁰, আর আমাদের rho ১,২২,০৬৭ ধাপে, ০.০৯ সেকেন্ডে সেটি পেয়েছে; √p ≈ ২,৩১,২১৬, তাই ধাপের সংখ্যা √p-এর একই মাপের। একই সংখ্যায় trial division-এ প্রায় 2.7 · 10¹⁰টি ভাগ লাগত। sympy.factorint একই উত্তর দিয়েছে; তার ভেতরে trial division, Pollard rho, Pollard-এর p − 1 আর elliptic curve পদ্ধতি পরপর চেষ্টা করা হয়।

$ uv run python verify.py
odd composites below 2*10^6 that pass base 2: 73 [2047, 3277, 4033, 4681, 8321, 15841] ...
... that pass bases 2 and 3: [1373653, 1530787, 1987021]
search time 0.7 s
1373653 {829: 1, 1657: 1} passes [2, 3, 17, 19, 29, 31, 43, 47] first witness 5
25326001 {2251: 1, 11251: 1} passes 5 first witness 7
3215031751 {151: 1, 751: 1, 28351: 1} passes [2, 3, 5, 7, 19, 37, 47] first witness 11
318665857834031151167461 {399165290221: 1, 798330580441: 1} passes 37 first witness 41
3317044064679887385961981 {1287836182261: 1, 2575672364521: 1} passes 41 first witness 43

verify.py §৪-এর তালিকা যাচাই করে। প্রথম অংশ 2 · 10⁶-এর নিচের প্রতিটি বিজোড় যৌগিক সংখ্যা দেখে: base 2-কে ৭৩টি ঠকায়, আর 2 ও 3 দুটোকেই মাত্র তিনটি, যার সবচেয়ে ছোটটি 1373653. তারপর তালিকার প্রতিটি সংখ্যার জন্য: sympy-র উৎপাদক, কোন মৌলিক base-এ (50-এর নিচে) পাস করে, আর প্রথম কোন base তাকে ধরে। একটি তালিকা ছাপা হলে তা সব পাস করা base; শুধু একটি সংখ্যা ছাপা হলে 2 থেকে সেই মৌলিক পর্যন্ত সব base-এ পাস। লক্ষ করুন, 3215031751 base 2, 3, 5, 7-এ পাস করে কিন্তু 11-এ ধরা পড়ে, আর শেষ সংখ্যাটি 2 থেকে 41 পর্যন্ত ১৩টি base-এ পাস করে, 43-এ ধরা পড়ে।

৭ · বাস্তবের আকার, রেকর্ড আর quantum প্রশ্ন

Pollard-এর rho ছোট উৎপাদক থাকলে দারুণ, কিন্তু দুটি সমান-আকারের বড় মৌলিকের গুণফলে √p নিজেই বিশাল: ৬০০ অঙ্কের n-এ p প্রায় ৩০০ অঙ্কের, √p প্রায় 10¹⁵⁰. এমন সংখ্যার জন্য আরও উন্নত পদ্ধতি আছে: quadratic sieve, elliptic curve method, আর সবচেয়ে শক্তিশালী number field sieve. এদের সময় exponential-এর চেয়ে কম কিন্তু polynomial-এর চেয়ে বেশি (“sub-exponential”); আজ পর্যন্ত সাধারণ কম্পিউটারে polynomial-সময়ের কোনো factoring algorithm জানা নেই।

রেকর্ডগুলো এই ধীর অগ্রগতি দেখায়। RSA Laboratories-এর factoring challenge-এর সংখ্যা RSA-250 (২৫০ দশমিক অঙ্ক, ৮২৯ bit) ২০২০ সালের ফেব্রুয়ারিতে Fabrice Boudot, Pierrick Gaudry, Aurore Guillevic, Nadia Heninger, Emmanuel Thomé আর Paul Zimmermann-এর দল number field sieve দিয়ে ভাঙে; প্রকাশিত হিসাবে মোট কাজ ছিল কয়েক হাজার CPU core-বছরের সমান। তার আগে RSA-129 (১২৯ অঙ্ক) ১৯৯৪ সালে ইন্টারনেটে শত শত স্বেচ্ছাসেবকের কম্পিউটারে quadratic sieve দিয়ে ভাঙা হয়েছিল। ছোট অঙ্কের এই রেকর্ডগুলো থেকে কোনো নির্দিষ্ট key-আকার “এখন ভাঙা যায়” বলে ধরে নেবেন না; নিরাপদ key-আকারের সিদ্ধান্ত মানদণ্ড-সংস্থার সুপারিশ আর যাচাই-করা লাইব্রেরির ওপর ছেড়ে দিন (Cybersecurity কোর্সের পর্ব ১০).

Quantum কম্পিউটার গল্পটি বদলাতে পারে। Peter Shor ১৯৯৪ সালে দেখান, একটি যথেষ্ট বড় ও নির্ভুল quantum কম্পিউটারে Shor-এর algorithm polynomial সময়ে factoring করতে পারে। মূল ধারণা এই কোর্সেরই গণিত: একটি এলোমেলো a-এর order mod n (পর্ব ২৪) বের করা, আর order জোড় হলে a^order/2 প্রায়ই 1-এর একটি অদ্ভুত বর্গমূল, যার gcd থেকে উৎপাদক বেরোয়, ঠিক §২-এর মতো। Quantum অংশটি শুধু order দ্রুত খুঁজে দেয়। এখন পর্যন্ত প্রকাশিত পরীক্ষাগুলোতে quantum যন্ত্রে এই পদ্ধতিতে খুব ছোট সংখ্যাই (যেমন 15) ভাঙা গেছে বলে জানা যায়, তবে সতর্কতা হিসেবে নতুন “post-quantum” মানদণ্ড ইতিমধ্যে তৈরি হচ্ছে।

৮ · পরীক্ষার অঙ্ক আর একটি ভুল প্রমাণ

বিশ্ববিদ্যালয়, চাকরি বা ইন্টারভিউর লিখিত পরীক্ষায় এই বিষয়ের যে ধাঁচের অঙ্ক আসে, তার তিনটি পূর্ণ সমাধান। প্রতিটি সংখ্যা scratch-এ Python দিয়ে মেলানো।

অঙ্ক ১: Miller–Rabin দিয়ে base 2-এ 341 পরীক্ষা করুন

1. 340 = 2² · 85, তাই s = 2, d = 85. দেখতে হবে 2⁸⁵ আর 2¹⁷⁰ mod 341.
2. 2¹⁰ = 1024 = 3 · 341 + 1 ≡ 1, তাই 2⁸⁵ = (2¹⁰)⁸ · 2⁵ ≡ 32. এটি 1 বা 340 নয়।
3. 32² = 1024 ≡ 1. অর্থাৎ 32 হলো 1-এর অদ্ভুত বর্গমূল, আর চেইনে −1 আসেনি। 341 যৌগিক; 2 একটি witness.
4. উৎপাদক: gcd(31, 341) = 31 আর gcd(33, 341) = 11, তাই 341 = 11 · 31. (Fermat base 2-এ 341 পাস করত, কারণ 2³⁴⁰ ≡ 1.)

অঙ্ক ২: base 5-এ 97 পরীক্ষা করুন

1. 96 = 2⁵ · 3, তাই s = 5, d = 3. 5³ = 125 ≡ 28 (mod 97).
2. বর্গ করতে থাকুন: 28² = 784 ≡ 8, 8² = 64, 64² = 4096 ≡ 22, 22² = 484 ≡ 96 ≡ −1.
3. r = 4-এ −1 এসেছে, তাই 97 base 5-এ strong probable prime. এটি একা প্রমাণ নয়; তবে 97 < 2047, আর §৪-এর তালিকা অনুযায়ী 2047-এর নিচে base 2-ই যথেষ্ট, তাই base 2 দিয়ে মিলিয়ে নিলে (96 = 2⁵ · 3, 2³ = 8, 64, 22, 96: আবার −1) 97 নিশ্চিতভাবে মৌলিক।

অঙ্ক ৩: ভুলের সম্ভাবনা 10⁻¹²-এর নিচে নামাতে কয়টি এলোমেলো base লাগবে?

1. Monier–Rabin সীমায় k base-এ ভুলের সম্ভাবনা ≤ 4^−k. চাই 4^−k ≤ 10⁻¹², অর্থাৎ k ≥ 12 · log 10 / log 4 ≈ 19.93.
2. তাই k = 20 যথেষ্ট (4⁻²⁰ ≈ 9.1 · 10⁻¹³). যদি কেউ Fermat-এর মতো 1/2 সীমা ভুল করে ধরেন, তাহলে 40টি লাগত বলে ভাববেন।

প্রমাণের ভুল ধরুন: “Miller–Rabin-এ পাস মানেই মৌলিক”

একজন লিখলেন: “দাবি: বিজোড় n আর কোনো base a-তে n যদি strong probable prime হয়, তবে n মৌলিক। প্রমাণ: চেইনটি a^d, a^2d, …, aⁿ⁻¹. শর্ত অনুযায়ী 1-এর ঠিক আগে −1 এসেছে, অর্থাৎ চেইনে 1-এর কোনো অদ্ভুত বর্গমূল নেই। Lemma ২.১ বলে 1-এর বর্গমূল শুধু ±1, তাই সব শর্ত মেলে, আর n মৌলিক। ∎” পড়ে ভুলটি নিজে খুঁজুন।

ভুল কোথায়: Lemma ২.১-এর শর্তই মৌলিক modulus. প্রমাণটি যা প্রমাণ করতে চায় (n মৌলিক), সেটিকেই Lemma লাগানোর সময় ধরে নিয়েছে: এটি চক্রাকার যুক্তি। তার ওপর, “চেইনে অদ্ভুত বর্গমূল দেখা যায়নি” মানে “অদ্ভুত বর্গমূল নেই” নয়; যৌগিক n-এ এরা থাকে, শুধু এই a-এর চেইন তাদের পার করেনি। একটি counterexample-ই যথেষ্ট: 2047 = 23 · 89, আর 2046 = 2 · 1023, 2¹⁰²³ ≡ 1 (mod 2047), তাই 2047 base 2-এ strong probable prime অথচ যৌগিক। সঠিক বক্তব্য এক দিকেই: strong test-এ ব্যর্থ হলে n যৌগিক।

৯ · AI-এর যুগে — মৌলিকতা আর factoring-এর দাবি যাচাই

“এই সংখ্যাটি কি মৌলিক?” প্রশ্নে AI সহকারী প্রায়ই আত্মবিশ্বাসী উত্তর দেয়, কিন্তু বড় সংখ্যায় মাথার হিসাবে আর পরীক্ষার শর্তে নিয়মিত ভুল করে। কয়েকটি ধরন, আর কীভাবে ধরবেন:

AI-এর উত্তর হাতে পেলে — যাচাইয়ের চেকলিস্ট

• n − 1 = 2^s · d নিজে মিলিয়েছি, d সত্যিই বিজোড়
• চেইনের প্রতিটি মান pow দিয়ে মিলিয়েছি
• “মৌলিক” দাবিটি probable নাকি প্রমাণিত, আর কোন সীমায়, তা স্পষ্ট
• দেওয়া উৎপাদক গুণ করে মূল সংখ্যা পেয়েছি, আর প্রতিটি উৎপাদক মৌলিক কি না দেখেছি
• কোডটি ছোট সব n-এ brute force-এর সঙ্গে মিলিয়েছি, আর পরিচিত counterexample (561, 2047) চালিয়েছি
• ইতিহাস, রেকর্ড আর সালগুলো বিশ্বস্ত উৎসে মিলিয়েছি

এই বিষয়ে যাচাই প্রায় বিনা খরচে হয়: একটি উৎপাদক গুণ করা বা একটি pow চালানো এক সেকেন্ডের কাজ। খোঁজা কঠিন, যাচাই সহজ, factoring-এর এই অসমতাই AI-এর উত্তর পরীক্ষা করার সবচেয়ে বড় সুবিধা।