Number theory — GCD, sieve আর modular arithmetic

১ · বিভাজ্যতা আর মৌলিক সংখ্যা — √n পর্যন্ত ভাগই যথেষ্ট

পরীক্ষার হলে ৩৬ জন পরীক্ষার্থীকে সমান সারিতে বসাতে হবে। ৪ সারিতে ৯ জন, ৬ সারিতে ৬ জন, ৯ সারিতে ৪ জন — হয়। কিন্তু ৩৭ জন হলে? ১ সারি বা ৩৭ সারি ছাড়া কোনো উপায় নেই। ৩৭ একটি মৌলিক সংখ্যা (prime). Number theory বা সংখ্যাতত্ত্ব পূর্ণসংখ্যার এমন ধর্ম নিয়ে কাজ করে, আর অ্যালগরিদমে এর প্রয়োগ সর্বত্র: hash-এর mod (পর্ব ৩২), কনটেস্টের “উত্তর 10⁹ + 7 দিয়ে ভাগশেষ”, আর ইন্টারনেটের নিরাপত্তা।

ভাষাটা আগে ঠিক করি। a | b মানে a দিয়ে b নিঃশেষে বিভাজ্য (b = a·k). ১-এর বড় যে সংখ্যার উৎপাদক শুধু ১ আর সে নিজে, সেটি মৌলিক; বাকিরা যৌগিক (composite). Fundamental theorem of arithmetic বলে, প্রতিটি সংখ্যা মৌলিকের গুণফলে ঠিক একভাবে ভাঙে: 360 = 2³·3²·5. মৌলিকরা তাই সংখ্যার “ইট”।

Trial division আর √n-এর যুক্তি

n মৌলিক কি না জানতে 2 থেকে n − 1 পর্যন্ত সব দিয়ে ভাগ করার দরকার নেই। যদি n = a·b আর a ≤ b হয়, তাহলে a·a ≤ a·b = n, অর্থাৎ a ≤ √n. তাই যৌগিক n-এর অন্তত একটি উৎপাদক √n-এর সমান বা ছোট — সেখানে কেউ ভাগ না করলে n মৌলিক। এটিই trial division, সময় O(√n). 36-এর উৎপাদক-জোড়া (1, 36), (2, 18), (3, 12), (4, 9), (6, 6) — প্রতিটি জোড়ার ছোটটি ≤ 6 = √36.

Pseudocode

IS-PRIME(n):
    if n < 2: return false
    for d = 2 while d * d <= n:        // d*d avoids floating-point sqrt
        if n mod d == 0: return false
    return true

মৌলিক কতগুলো? Euclid-এর Elements-এর (আনুমানিক খ্রিস্টপূর্ব ৩০০) প্রমাণ — proof by contradiction: ধরুন সসীম, p₁, …, pₖ. N = p₁·p₂·…·pₖ + 1-কে কোনো pᵢ ভাগ করলে ভাগশেষ 1, অথচ N-এর একটি মৌলিক উৎপাদক থাকতেই হবে — যা তালিকায় নেই। স্ববিরোধ। আর তারা কত ঘন? Prime number theorem বলে n পর্যন্ত মৌলিক প্রায় n / ln n টি: 10⁶ পর্যন্ত আসলে 78498টি, সূত্রে ≈ 72382 — কাছাকাছি, আর n বড় হলে অনুপাত 1-এর দিকে যায়।

২ · Sieve of Eratosthenes — একবারে সব মৌলিক

10⁶ পর্যন্ত প্রতিটি সংখ্যায় আলাদা করে trial division চালালে মোট কাজ মোটামুটি n√n / ln n-এর ক্রমে — অনেক। উল্টো পথে ভাবুন: মৌলিক খোঁজার বদলে যৌগিকগুলো কেটে দিন। চালনিতে আটা চালার মতো — গুণিতকরা ঝরে পড়ে, মৌলিকরা থেকে যায়। এটিই sieve of Eratosthenes: 2 থেকে n লিখুন; সবচেয়ে ছোট না-কাটা সংখ্যা p মৌলিক; তার গুণিতকগুলো p² থেকে কাটুন (2p, 3p, … (p−1)p আগেই ছোট মৌলিকের হাতে কাটা পড়েছে); p² > n হলে থামুন। n = 30:

O(n log log n) কেন — একটি স্কেচ

মৌলিক p-এর জন্য কাটার চেষ্টা প্রায় n/p. মোট কাজ n·(1/2 + 1/3 + 1/5 + 1/7 + …), শুধু মৌলিক p ≤ n-এর ওপর। সব k-এর ওপর 1/k-এর যোগফল harmonic series, ≈ ln n (পর্ব ০৬). মৌলিকরা অনেক বিরল — k-তম মৌলিক প্রায় k ln k — তাই তাদের উল্টোর যোগফল অনেক ধীরে বাড়ে: Mertens-এর (১৮৭৪) ফল অনুযায়ী ∑ 1/p ≈ ln ln n + 0.26. সুতরাং:

ln ln 10⁶ ≈ 2.6, অর্থাৎ বাস্তবে এটি প্রায় রৈখিক। লুপের শর্ত p * p <= n-এ n ≈ 2 × 10⁹ হলে int-এ গুণফল overflow করে — long long দিন বা p <= n / p লিখুন। জায়গা বাঁচাতে C++-এ vector<bool> বা bitset (প্রতি সংখ্যায় এক bit), আর শুধু বিজোড় সংখ্যা রাখলে আরও অর্ধেক।

৩ · Linear sieve আর smallest prime factor — উৎপাদকে ভাঙা O(log n)-এ

stepper-এ 30 কাটা হয়েছে তিনবার (2, 3, 5 — তিনজনই তার উৎপাদক)। Linear sieve (Gries ও Misra, ১৯৭৮) প্রতিটি যৌগিক x-কে কাটে ঠিক একবার — তার smallest prime factor spf[x] দিয়ে। নিয়ম: প্রতিটি i-এর জন্য এ পর্যন্ত পাওয়া মৌলিক p = 2, 3, 5, … নিয়ে spf[i·p] = p লিখুন, কিন্তু p > spf[i] হলেই থামুন।

কেন ঠিক একবার? যেকোনো যৌগিক x-কে একভাবেই x = i·p লেখা যায় যেখানে p = spf[x] আর i = x / p; তখন i-এর সব উৎপাদক ≥ p, তাই p ≤ spf[i] — শর্ত মানে, x লেখা হয় i-এর ধাপে। অন্য কোনো জোড়া (i′, p′) যেখানে p′ > spf[x], সেখানে spf[i′] = spf[x] < p′, তাই লুপ আগেই থেমে যায়। ফলে মোট লেখা = যৌগিকের সংখ্যা, সময় Θ(n), জায়গা Θ(n) int (spf টেবিল)।

spf টেবিলের আসল লাভ: x-কে উৎপাদকে ভাঙতে বারবার p = spf[x] দিয়ে ভাগ করুন। প্রতিটি ভাগে x অন্তত অর্ধেক হয়, তাই ধাপ ≤ log₂ x — 10⁵টি প্রশ্নেও মুহূর্তের কাজ, যেখানে প্রতিটিতে trial division হলে O(√x). ছোট উদাহরণ: 2 থেকে 12-এর spf = 2, 3, 2, 5, 2, 7, 2, 3, 2, 11, 2. তাই 12 → (÷2) 6 → (÷2) 3 → (÷3) 1, অর্থাৎ 12 = 2²·3 — তিনটি ভাগ, কোনো খোঁজাখুঁজি নেই।

#include <initializer_list>
#include <iostream>
#include <vector>
using namespace std;

int main() {
    const int N = 1000000;
    vector<int> spf(N + 1, 0), primes;              // spf[x] = smallest prime factor of x
    long long marks = 0;                             // how many times a composite is written
    for (int i = 2; i <= N; ++i) {
        if (spf[i] == 0) { spf[i] = i; primes.push_back(i); }    // nobody marked i: prime
        for (int p : primes) {                       // x = i * p gets spf[x] = p
            if (p > spf[i] || (long long)i * p > N) break;
            spf[i * p] = p;
            ++marks;
        }
    }
    cout << "primes <= " << N << ": " << primes.size() << ", largest " << primes.back() << '\n';
    cout << "composites marked: " << marks << " (= " << N - 1 - (long long)primes.size() << " composites)\n";
    for (int x : {360, 2026, 997920, 999983}) {
        cout << x << " =";
        for (char sep = ' '; x > 1; sep = '*') {     // O(log x): each step divides by >= 2
            int p = spf[x], e = 0;
            while (x % p == 0) { x /= p; ++e; }
            cout << sep << p;
            if (e > 1) cout << '^' << e;
        }
        cout << '\n';
    }
}

primes <= 1000000: 78498, largest 999983
composites marked: 921501 (= 921501 composites)
360 = 2^3*3^2*5
2026 = 2*1013
997920 = 2^5*3^4*5*7*11
999983 = 999983

N = 10**6
spf = [0] * (N + 1)                    # spf[x] = smallest prime factor of x
primes, marks = [], 0                  # marks = how many times a composite is written
for i in range(2, N + 1):
    if spf[i] == 0:                    # nobody marked i: prime
        spf[i] = i
        primes.append(i)
    for p in primes:                   # x = i * p gets spf[x] = p
        if p > spf[i] or i * p > N:
            break
        spf[i * p] = p
        marks += 1
print(f"primes <= {N}: {len(primes)}, largest {primes[-1]}")
print(f"composites marked: {marks} (= {N - 1 - len(primes)} composites)")
for x in (360, 2026, 997920, 999983):
    parts, y = [], x
    while y > 1:                       # O(log x): each step divides by >= 2
        p, e = spf[y], 0
        while y % p == 0:
            y, e = y // p, e + 1
        parts.append(f"{p}^{e}" if e > 1 else str(p))
    print(f"{x} = " + "*".join(parts))

primes <= 1000000: 78498, largest 999983
composites marked: 921501 (= 921501 composites)
360 = 2^3*3^2*5
2026 = 2*1013
997920 = 2^5*3^4*5*7*11
999983 = 999983

৪ · Euclid-এর GCD — প্রমাণ, Lamé-এর সীমা আর LCM

252 সেমি × 198 সেমি একটি ঘরের মেঝেতে সবচেয়ে বড় কোন মাপের বর্গাকার টাইল বসালে কাটতে হবে না? উত্তর দুই দৈর্ঘ্যের gcd — গসাগু। স্কুলে শেখা “দুটিকেই উৎপাদকে ভাঙো” পদ্ধতি বড় সংখ্যায় অচল, কারণ উৎপাদকে ভাঙাই কঠিন। Euclid-এর অ্যালগরিদম উৎপাদক ছাড়াই কাজ করে: মেঝে থেকে যত বড় সম্ভব বর্গ কেটে নিন, বাকি আয়তক্ষেত্রে একই কাজ করুন।

প্রমাণ। a = q·b + r লিখুন, r = a mod b. d যদি a আর b দুটিকেই ভাগ করে, তাহলে r = a − q·b-কেও ভাগ করে। উল্টো দিকে d যদি b আর r-কে ভাগ করে, তাহলে a = q·b + r-কেও ভাগ করে। অর্থাৎ (a, b) আর (b, r)-এর সাধারণ উৎপাদকের সেট হুবহু এক — তাই সবচেয়ে বড়টিও এক। প্রতিটি ধাপে দ্বিতীয় সংখ্যা কঠোরভাবে কমে (r < b), তাই অ্যালগরিদম শেষ হয়।

কত ধাপ? O(log min(a, b)) আর Lamé

মূল পর্যবেক্ষণ: a ≥ b হলে a mod b < a/2. কারণ b ≤ a/2 হলে ভাগশেষ < b ≤ a/2; আর b > a/2 হলে a mod b = a − b < a/2. জোড়া (a, b) → (b, r) → (r, …), তাই দুই ধাপে প্রথম সংখ্যা অন্তত অর্ধেক — মোট ধাপ ≤ 2 log₂ min(a, b) + 1, অর্থাৎ O(log min(a, b)) ভাগ। 10¹⁸-এর মতো সংখ্যাতেও ১০০-এর কম ধাপ।

সবচেয়ে খারাপ ইনপুট কোনটি? প্রতিটি ধাপে q = 1 হলে সংখ্যা সবচেয়ে ধীরে কমে — ঠিক পরপর দুটি Fibonacci সংখ্যা: gcd(89, 55) → (55, 34) → (34, 21) → … → (1, 0), ৯টি ধাপ। Lamé ১৮৪৪ সালে এটি প্রমাণ করেন: ধাপের সংখ্যা ছোট সংখ্যার দশমিক অঙ্কের ৫ গুণের বেশি নয়। অনেকে একে কোনো অ্যালগরিদমের সময়-বিশ্লেষণের প্রথম দিকের উদাহরণ বলেন।

LCM আর overflow — আগে ভাগ, পরে গুণ

lcm(a, b) · gcd(a, b) = a · b, তাই lcm = a·b / gcd. কিন্তু C++-এ a * b / g লিখলে a, b ≈ 10¹² হলে গুণফল 10²⁴ — long long (≈ 9.2 × 10¹⁸) ছাড়িয়ে যায়, যদিও উত্তর হয়তো ছোট। সঠিক ক্রম a / g * b: g দিয়ে a নিঃশেষে বিভাজ্য, তাই ভাগে কিছু হারায় না। C++17-এর std::gcd, std::lcm (<numeric>) আর Python-এর math.gcd, math.lcm তৈরি আছে; অনেক সংখ্যার gcd = gcd(gcd(a, b), c) … দৈনন্দিন একটি ব্যবহার: 1920 × 1080 পর্দার অনুপাত — gcd = 120, তাই 16 : 9.

৫ · Extended Euclid — Bézout-এর x আর y

Euclid-এর ধাপগুলো উল্টো দিকে জুড়ে দিলে আরও একটি তথ্য পাওয়া যায়। Bézout's identity: যেকোনো a, b-এর জন্য এমন পূর্ণসংখ্যা x, y আছে যাতে

Extended Euclidean algorithm এই x, y বের করে। recursive রূপের যুক্তি: gcd(b, a mod b)-এর জন্য b·x′ + (a mod b)·y′ = g জানা থাকলে, a mod b = a − ⌊a/b⌋·b বসিয়ে সাজালে a·y′ + b·(x′ − ⌊a/b⌋·y′) = g. তাই x = y′, y = x′ − ⌊a/b⌋·y′; ভিত্তি b = 0-এ (x, y) = (1, 0). Iterative রূপে প্রতিটি সারিতে একটি সংখ্যা r আর তার জোড়া (x, y) রাখি, এই invariant ধরে: r = a·x + b·y. প্রতিটি নতুন সারি = আগের-আগের সারি − q × আগের সারি, তাই invariant টিকে থাকে। নিচে a = 3120, b = 17 (এই সংখ্যাগুলো অংশ ৮-এর খেলনা RSA-তে আবার আসবে):

প্রয়োগ — linear Diophantine equation: a·x + b·y = c-এর পূর্ণসংখ্যা সমাধান থাকে ঠিক তখন, যখন g = gcd(a, b) দিয়ে c বিভাজ্য। একটি সমাধান (x·c/g, y·c/g), আর সব সমাধান x + k·b/g, y − k·a/g. ধরুন ৭ টাকা আর ৫ টাকার স্ট্যাম্প দিয়ে ঠিক ৩১ টাকা বানাতে হবে: gcd = 1 | 31, তাই সমাধান আছে — যেমন 3·7 + 2·5 = 31.

৬ · Modular arithmetic — ভাগশেষের হিসাব, দ্রুত ঘাত আর inverse

আজ সোমবার। ১০০ দিন পর কী বার? ১০০ দিন গোনার দরকার নেই: 100 mod 7 = 2, তাই সোমবার + ২ = বুধবার। সপ্তাহের বার, ঘড়ির ঘণ্টা — এগুলো modular arithmetic: একটি সীমা m-এ পৌঁছে আবার শূন্য থেকে। লেখা হয় a ≡ b (mod m), মানে m | (a − b).

কেন মাঝপথে mod নেওয়া নিরাপদ? a = q₁·m + r₁ আর b = q₂·m + r₂ লিখুন। তাহলে a + b = (q₁ + q₂)·m + (r₁ + r₂) আর a·b = (q₁q₂m + q₁r₂ + q₂r₁)·m + r₁·r₂ — m-এর গুণিতক অংশ ভাগশেষে কোনো প্রভাব ফেলে না, তাই শুধু r₁, r₂ নিয়ে হিসাব করলেই চলে। ভাগে এটি ভাঙে: 12 / 4 = 3, অথচ mod 8-এ 12 ≡ 4, আর 4 / 4 = 1 ≠ 3. এ কারণেই ভাগের জন্য inverse লাগে।

কনটেস্টে উত্তর কেন “mod 10⁹ + 7”

“কতভাবে সাজানো যায়” ধরনের উত্তর প্রায়ই শত শত অঙ্কের। পুরো সংখ্যা ছাপাতে big integer লাগত, আর সময় বাড়ত। তাই প্রশ্নে বলা থাকে “উত্তর 10⁹ + 7 দিয়ে ভাগশেষ দাও”। এই সংখ্যাটি বাছার কারণ: (১) এটি মৌলিক, তাই 1…p − 1-এর প্রতিটির inverse আছে — ভাগ করা যায়; (২) এটি 2³⁰-এর চেয়ে ছোট, তাই দুটি ভাগশেষের যোগফল 32-bit signed int-এ আঁটে; (৩) দুটির গুণফল < 2⁶³, তাই long long-এ আঁটে। আরেকটি প্রচলিত মৌলিক 998244353 = 119·2²³ + 1, যা number-theoretic transform-এর (দ্রুত বহুপদী গুণ) উপযোগী।

দ্রুত modular ঘাত

aᵉ mod m-এ e-বার গুণ নয়: পর্ব ৪৯-এর (Divide and conquer) binary exponentiation — e-কে binary-তে ভেঙে বারবার বর্গ, প্রতিটি ধাপে mod. গুণ লাগে O(log e), আর m ≤ 2⁶³ হলে প্রতিটি গুণ 128-bit-এ (__int128) করে mod নিতে হয়। Python-এ তৈরি আছে: pow(a, e, m). ছোট উদাহরণ, 3¹³ mod 7: 13 = 1101₂; বারবার বর্গে 3¹ = 3, 3² ≡ 2, 3⁴ ≡ 4, 3⁸ ≡ 2; যে bit-গুলো 1 (8, 4, 1), তাদের গুণ: 2·4·3 = 24 ≡ 3. ১২টি গুণের বদলে ৫টি।

Modular inverse — দুই পথ

a-এর modular inverse x, যাতে a·x ≡ 1 (mod m). যেমন mod 7-এ 3⁻¹ = 5, কারণ 15 = 2·7 + 1. inverse থাকে ঠিক তখন, যখন gcd(a, m) = 1 (gcd = g > 1 হলে a·x − k·m সবসময় g-এর গুণিতক, কখনো 1 নয়)।

• পথ ১ — extended Euclid: a·x + m·y = 1 থেকে mod m নিলে a·x ≡ 1. যেকোনো m-এ চলে (শুধু gcd = 1 লাগে), O(log m).
• পথ ২ — Fermat: Fermat's little theorem বলে, p মৌলিক আর p ∤ a হলে aᵖ⁻¹ ≡ 1 (mod p). তাই a · aᵖ⁻² ≡ 1, অর্থাৎ a⁻¹ ≡ aᵖ⁻² (mod p) — একটি দ্রুত ঘাত, O(log p). শুধু মৌলিক mod-এ।

Fermat-এর উপপাদ্যের ছোট প্রমাণ: 1, 2, …, p − 1-কে a দিয়ে গুণ করে mod p নিন। কোনো দুটি সমান হতে পারে না (a·i ≡ a·j হলে p | a(i − j), p ∤ a বলে p | (i − j), তাই i = j), আর কেউ 0 নয়। তাই a, 2a, …, (p−1)a আসলে 1, …, p − 1-এরই অন্য ক্রম। দুই তালিকার গুণফল সমান: aᵖ⁻¹·(p−1)! ≡ (p−1)!, আর (p−1)!-এর inverse আছে বলে কেটে দিলে aᵖ⁻¹ ≡ 1. একটি উদাহরণে কাজে লাগাই: 3²⁰⁰ mod 13? 3¹² ≡ 1, আর 200 = 16·12 + 8, তাই 3²⁰⁰ ≡ 3⁸ = 6561 ≡ 9.

৭ · nCr mod p, Euler-এর φ আর Chinese remainder theorem

Factorial টেবিল দিয়ে nCr mod p

“n জনের ক্রিকেট স্কোয়াড থেকে r জনের দল কতভাবে?” — C(n, r) = n! / (r!·(n−r)!). mod p-এ ভাগ মানে inverse দিয়ে গুণ। অনেক প্রশ্নের জন্য একবার টেবিল বানান: fact[i] = fact[i−1]·i; তারপর শুধু inv_fact[N] = fact[N]ᵖ⁻² একটি দ্রুত ঘাতে, আর নিচে নামতে inv_fact[i−1] = inv_fact[i]·i (কারণ 1/(i−1)! = i/i!). এরপর প্রতিটি প্রশ্ন:

প্রস্তুতি O(N + log p), প্রতিটি প্রশ্ন O(1), জায়গা O(N). শর্ত: n < p, নইলে n!-এ p নিজেই ঢুকে 0 করে দেয়; n ≥ p হলে লাগে Lucas's theorem. কনটেস্টে n ≤ 10⁶ আর p = 10⁹ + 7 — শর্ত নিশ্চিন্তে মানে।

Euler-এর totient φ(n) আর Euler-এর উপপাদ্য

φ(n) = 1…n-এর মধ্যে n-এর সাথে সহমৌলিক সংখ্যা কয়টি। φ(12) = 4 (1, 5, 7, 11); মৌলিক p-এ φ(p) = p − 1. সূত্র: n-এর প্রতিটি আলাদা মৌলিক উৎপাদক p-এর জন্য “p-এর গুণিতক বাদ” — φ(n) = n·∏(1 − 1/p). 36 = 2²·3², তাই φ(36) = 36 · ½ · ⅔ = 12. একটি n-এর জন্য trial division-এ O(√n); 1…N সবার জন্য sieve-এর মতো করে O(N log log N).

Euler-এর উপপাদ্য: gcd(a, n) = 1 হলে a^φ(n) ≡ 1 (mod n) — Fermat-এর সাধারণ রূপ, প্রমাণও একই (সহমৌলিক φ(n)টি সংখ্যাকে a দিয়ে গুণ করলে আবার তারাই)। উদাহরণ: 5⁴ = 625 = 52·12 + 1 ≡ 1 (mod 12). ফলে যেকোনো mod n-এ inverse = a^(φ(n)−1), আর খুব বড় ঘাতকে ছোট করা যায়: aᵉ ≡ a^(e mod φ(n)).

Chinese remainder theorem — Sunzi-র ধাঁধা

প্রাচীন চীনা গণিতগ্রন্থ Sunzi Suanjing-এর (খ্রিস্টীয় তৃতীয় থেকে পঞ্চম শতকের মধ্যে কোনো সময়) একটি প্রশ্ন: কিছু জিনিস আছে, সংখ্যা অজানা। ৩টি করে গুনলে ২টি বাকি থাকে, ৫টি করে গুনলে ৩টি, ৭টি করে গুনলে ২টি। জিনিস কয়টি?

Chinese remainder theorem (CRT) বলে, mod-গুলো জোড়ায় জোড়ায় সহমৌলিক হলে M = 3·5·7 = 105-এর ভেতরে উত্তর ঠিক একটি। গঠন: প্রতিটি mᵢ-এর জন্য Mᵢ = M / mᵢ, আর yᵢ = Mᵢ⁻¹ mod mᵢ; তাহলে x = ∑ rᵢ·Mᵢ·yᵢ mod M. এখানে M₁ = 35 ≡ 2 (mod 3), y₁ = 2; M₂ = 21 ≡ 1 (mod 5), y₂ = 1; M₃ = 15 ≡ 1 (mod 7), y₃ = 1. x = 2·35·2 + 3·21·1 + 2·15·1 = 140 + 63 + 30 = 233 ≡ 23 (mod 105). যাচাই: 23 = 7·3 + 2 = 4·5 + 3 = 3·7 + 2. কেন কাজ করে: rᵢ·Mᵢ·yᵢ পদটি mod mᵢ-এ rᵢ, আর বাকি সব mᵢ-এর গুণিতক বলে 0. প্রতিটি inverse extended Euclid-এ আসে, তাই মোট সময় O(k log M). mod-গুলো সহমৌলিক না হলে সাবধান: x ≡ r₁ (mod m₁) আর x ≡ r₂ (mod m₂)-এর সমাধান থাকে ঠিক তখন, যখন gcd(m₁, m₂) দিয়ে r₂ − r₁ বিভাজ্য; তখন extended Euclid দিয়ে দুটি মিলে একটি শর্ত হয়, mod lcm(m₁, m₂)-এ। যেমন x ≡ 1 (mod 4) আর x ≡ 2 (mod 6)-এর কোনো সমাধান নেই — প্রথমটি বলে x বিজোড়, দ্বিতীয়টি বলে জোড়। বাস্তবে CRT দিয়ে বড় mod-এর হিসাব কয়েকটি ছোট মৌলিক mod-এ ভাগ করে সমান্তরালে চালানো যায় (যেমন RSA-র ডিক্রিপশন দ্রুত করতে)।

৮ · Miller–Rabin — বড় সংখ্যা মৌলিক কি না, দ্রুত

n ≈ 10¹⁸ হলে trial division-এ 10⁹ ভাগ — কনটেস্টে অসম্ভব। Fermat-এর উপপাদ্য উল্টে একটি পরীক্ষা বানানো যায়: কোনো a-তে aⁿ⁻¹ ≢ 1 (mod n) হলে n নিশ্চিতভাবে যৌগিক। কিন্তু উল্টোটা সত্য নয়। 561 = 3·11·17 যৌগিক, তবু 561-এর সাথে সহমৌলিক প্রতিটি a-তে a⁵⁶⁰ ≡ 1 (mod 561). এমন সংখ্যাকে বলে Carmichael number, আর এরা অসীম সংখ্যক — Fermat-এর পরীক্ষা এদের কাছে অন্ধ।

Miller–Rabin আরেকটি সত্য যোগ করে: মৌলিক p-এ x² ≡ 1 মানে p | (x − 1)(x + 1), তাই x ≡ ±1 — 1-এর বর্গমূল শুধু 1 আর −1. n − 1 = d·2ˢ লিখুন (d বিজোড়)। x = a^d থেকে s − 1 বার বর্গ করে দেখুন: n মৌলিক হলে ক্রমটি হয় 1 দিয়ে শুরু হয়, নয়তো কোথাও n − 1 আসে। n − 1 না এসে 1-এ পৌঁছানো মানে 1-এর একটি “অস্বাভাবিক” বর্গমূল পাওয়া গেল — n নিশ্চিতভাবে যৌগিক, আর a হলো তার witness (সাক্ষী)।

কতটা নির্ভরযোগ্য? যৌগিক n-এর জন্য 1…n − 1-এর অন্তত ¾ ভাগ a witness (Rabin, ১৯৮০) — এলোমেলো k টি base-এ ভুলের সম্ভাবনা ≤ 4⁻ᵏ: একটি Monte Carlo অ্যালগরিদম, যার ভুল একদিকে (“মৌলিক” উত্তর ভুল হতে পারে, “যৌগিক” কখনো নয়)। আর 64-bit-এর জন্য এলোমেলোতাও লাগে না: প্রথম ১২টি মৌলিক {2, 3, 5, …, 37} base নিলে প্রতিটি n < 2⁶⁴-এ ফল নির্ভুল — পরীক্ষিত সীমাটি 2⁶⁴-এর চেয়ে অনেক বড়, তাই তখন এটি deterministic. প্রতিটি base-এ একটি দ্রুত ঘাত আর ≤ s বর্গ: সময় O(k log n) modular গুণ। n < 2⁶⁴ হলে গুণফল 2¹²⁸ পর্যন্ত যায়, তাই গুণ unsigned __int128-এ।

#include <initializer_list>
#include <iostream>
#include <vector>
using namespace std;
using ll = long long; using u64 = unsigned long long; using u128 = unsigned __int128;
const ll MOD = 1'000'000'007;
ll extgcd(ll a, ll b, ll& x, ll& y) {               // returns g and sets a*x + b*y = g
    if (b == 0) { x = 1; y = 0; return a; }
    ll g = extgcd(b, a % b, y, x);                  // b*y + (a%b)*x = g
    y -= (a / b) * x;
    return g;
}
u64 power(u64 b, u64 e, u64 m) {                     // b^e mod m; products in 128 bits
    u64 r = 1; b %= m;
    for (; e; e >>= 1, b = (u128)b * b % m) if (e & 1) r = (u128)r * b % m;
    return r;
}
bool isPrime(u64 n) {                                // deterministic Miller-Rabin, 64-bit n
    if (n < 2) return false;
    u64 d = n - 1; int s = 0;
    while (d % 2 == 0) { d /= 2; ++s; }              // n - 1 = d * 2^s, d odd
    for (u64 a : {2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37}) {
        if (n % a == 0) return n == a;
        u64 x = power(a, d, n); bool witness = x != 1 && x != n - 1;
        for (int r = 1; r < s && witness; ++r) witness = (x = (u128)x * x % n) != n - 1;
        if (witness) return false;                   // a proves n composite
    }
    return true;
}
int main() {
    ll x, y, g = extgcd(3120, 17, x, y);
    cout << "gcd=" << g << " x=" << x << " y=" << y << ", 17^-1 mod 3120 = " << (y % 3120 + 3120) % 3120 << '\n';
    extgcd(123456789, MOD, x, y);
    cout << "inverse of 123456789: extgcd " << (x % MOD + MOD) % MOD << ", fermat " << power(123456789, MOD - 2, MOD) << '\n';
    const int N = 1000000; vector<ll> f(N + 1, 1), inv(N + 1);
    for (int i = 1; i <= N; ++i) f[i] = f[i - 1] * i % MOD;
    inv[N] = power(f[N], MOD - 2, MOD);              // one Fermat inverse, then walk down
    for (int i = N; i > 0; --i) inv[i - 1] = inv[i] * i % MOD;
    auto C = [&](int n, int r) { return r < 0 || r > n ? 0 : f[n] * inv[r] % MOD * inv[n - r] % MOD; };
    cout << "C(10,3)=" << C(10, 3) << ", C(1000000,500000) mod p = " << C(N, N / 2) << '\n';
    for (u64 n : {561ULL, 3215031751ULL, 1000000007ULL, 999999999989ULL, 18446744073709551557ULL})
        cout << n << (isPrime(n) ? " prime" : " composite") << '\n';
}

gcd=1 x=2 y=-367, 17^-1 mod 3120 = 2753
inverse of 123456789: extgcd 18633540, fermat 18633540
C(10,3)=120, C(1000000,500000) mod p = 996692777
561 composite
3215031751 composite
1000000007 prime
999999999989 prime
18446744073709551557 prime

MOD = 10**9 + 7

def extgcd(a, b):                      # returns (g, x, y) with a*x + b*y = g
    if b == 0:
        return a, 1, 0
    g, x1, y1 = extgcd(b, a % b)
    return g, y1, x1 - (a // b) * y1

def is_prime(n):                       # deterministic Miller-Rabin, 64-bit n
    if n < 2:
        return False
    d, s = n - 1, 0
    while d % 2 == 0:                  # n - 1 = d * 2^s, d odd
        d, s = d // 2, s + 1
    for a in (2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37):
        if n % a == 0:
            return n == a
        x = pow(a, d, n)
        if x in (1, n - 1):            # this base says "probably prime"
            continue
        for _ in range(s - 1):
            x = x * x % n
            if x == n - 1:
                break
        else:                          # never reached n - 1: a is a witness
            return False
    return True

g, x, y = extgcd(3120, 17)
print(f"gcd={g} x={x} y={y}, 17^-1 mod 3120 = {y % 3120}")
_, x, _ = extgcd(123456789, MOD)
print(f"inverse of 123456789: extgcd {x % MOD}, fermat {pow(123456789, MOD - 2, MOD)}")
N = 10**6
f, inv = [1] * (N + 1), [1] * (N + 1)
for i in range(1, N + 1):
    f[i] = f[i - 1] * i % MOD
inv[N] = pow(f[N], MOD - 2, MOD)       # one Fermat inverse, then walk down
for i in range(N, 0, -1):
    inv[i - 1] = inv[i] * i % MOD
C = lambda n, r: 0 if r < 0 or r > n else f[n] * inv[r] % MOD * inv[n - r] % MOD
print(f"C(10,3)={C(10, 3)}, C(1000000,500000) mod p = {C(N, N // 2)}")
for n in (561, 3215031751, 1000000007, 999999999989, 18446744073709551557):
    print(n, "prime" if is_prime(n) else "composite")

gcd=1 x=2 y=-367, 17^-1 mod 3120 = 2753
inverse of 123456789: extgcd 18633540, fermat 18633540
C(10,3)=120, C(1000000,500000) mod p = 996692777
561 composite
3215031751 composite
1000000007 prime
999999999989 prime
18446744073709551557 prime

৯ · দৈনন্দিন জীবনে সংখ্যাতত্ত্ব আর ইতিহাস

১০ · AI-এর যুগে number theory — overflow, ঋণাত্মক mod আর অর্ধেক-সত্য

sieve, gcd, modpow — এসব AI মুহূর্তে লিখে দেয়, আর ছোট উদাহরণে প্রায় সবসময় ঠিক চলে। ভুল লুকায় বড় সংখ্যায় আর কিনারায়: int-এ গুণ, C++-এর ঋণাত্মক %, p যৌগিক হলেও Fermat-এর inverse, আর “Fermat test দিয়ে মৌলিক যাচাই”-এর মতো অর্ধেক-সত্য।

Number theory-তে AI সাধারণত যেসব ভুল করে

AI-এর (বা নিজের) number theory কোড যাচাইয়ের চেকলিস্ট

• প্রতিটি গুণের আগে দেখেছি দুই পাশ কত বড় হতে পারে, আর দরকারে long long বা __int128 দিয়েছি
• ছোট একটি উদাহরণ (gcd(252, 198), 17⁻¹ mod 3120) হাতে করে প্রোগ্রামের সাথে মিলিয়েছি
• sieve-এর ফল যাচাই করেছি: 10⁶ পর্যন্ত 78498টি মৌলিক, 100 পর্যন্ত 25টি
• isPrime-কে 1…10⁵ পর্যন্ত trial division-এর সাথে stress test করেছি, আর 561, 3215031751 আলাদা করে দিয়েছি
• inverse-এর প্রতিটি ফল a·x mod m = 1 কি না পরীক্ষা করেছি, আর mod মৌলিক কি না নিশ্চিত হয়েছি
• nCr-এ n < p শর্ত আর r > n বা r < 0-তে 0 ফেরত দেওয়া — দুটিই আছে