Concept Notes — Advanced Number Theory (gcd-এর ওপারে)¶
Level 1-2-তে যে hatiyar গুলো বানিয়েছিলে, এই note-এ সেগুলোকে শান দেব। মূল সুর একটাই: প্রতিটা advanced technique আসলে একটা পুরোনো idea-র চালাক পুনর্ব্যবহার — extended Euclid হলো Euclid-এর পথ মনে রাখা, matrix exponentiation হলো binary exponentiation-এর সংখ্যার বদলে matrix-এ চালানো, Mobius হলো inclusion-exclusion-এর ছোট্ট encode।
1. Extended Euclid — সিঁড়ি বেয়ে নেমে, আবার উঠে আসা¶
সাধারণ Euclid gcd(a, b) দেয়; extended Euclid আরো দুটো সংখ্যা দেয় — x আর y, যাতে a·x + b·y = gcd(a, b)। কেন দরকার? কারণ এই x-ই modular inverse-এর চাবি (section 6), আর Diophantine-এর প্রথম solution (section 2)।
পদ্ধতি: Euclid-এর সিঁড়ি বেয়ে নামো, তারপর back-substitution-এ উঠে এসো। gcd(30, 18) দিয়ে দেখি:
নামা (Euclid): ওঠা (back-substitution):
30 = 1·18 + 12 6 = 18 − 1·12
18 = 1·12 + 6 = 18 − 1·(30 − 1·18) [12-কে বদলালাম]
12 = 2·6 + 0 → gcd = 6 = 2·18 − 1·30
⇒ 30·(−1) + 18·(2) = 6 ✓ x = −1, y = 2
প্রতিটা remainder-কে তার আগের লাইনের ভাষায় লিখে উপরে উঠতে থাকো — এটাই back-substitution-এর সিঁড়ি। Code-এ recursion এই কাজটাই নিজে করে:
def ext_gcd(a, b):
if b == 0:
return a, 1, 0 # a·1 + 0·0 = a
g, x1, y1 = ext_gcd(b, a % b) # b·x1 + (a % b)·y1 = g
return g, y1, x1 - (a // b) * y1
print(ext_gcd(30, 18)) # (6, -1, 2) → 30·(−1) + 18·2 = 6 ✓
Update rule-এর রহস্য: a % b = a − (a // b)·b বসিয়ে দিলে b·x1 + (a − (a//b)·b)·y1 = a·y1 + b·(x1 − (a//b)·y1) — মানে নতুন x = y1, নতুন y = x1 − (a//b)·y1। এক লাইনের algebra, কিন্তু একবার খাতায় না করলে বিশ্বাস হয় না — করে নাও।
2. Linear Diophantine — ax + by = c কখন সম্ভব?¶
a·x + b·y আকারের সব সংখ্যাই gcd(a, b)-এর গুণিতক (কারণ gcd দুটোকেই ভাগ করে)। তাই:
Solvable হলে: ext_gcd থেকে a·x0 + b·y0 = g পাও, দুপাশে c/g গুণ করো — এক solution তৈরি। আর বাকি অসীম solution একটা পরিবার:
def diophantine(a, b, c):
g, x0, y0 = ext_gcd(a, b)
if c % g != 0:
return None # solution নেই
k = c // g
return x0 * k, y0 * k # একটা solution
# সব solution: x = x0 + t·(b/g), y = y0 − t·(a/g), t = যেকোনো integer
print(diophantine(30, 18, 12)) # (-2, 4) → 30·(−2) + 18·4 = 12 ✓
পরিবারের ছবিটা ভাবো: একটা solution পেলে x-কে b/g বাড়িয়ে y-কে a/g কমালে যোগফল অপরিবর্তিত — দাঁড়িপাল্লার দুই পাশে সমান ওজন বদল।
3. CRT — দুটো ঘড়ি মেলানোর ধাঁধা¶
ধরো একটা ঘটনার সময় জানো দুটো ভাঙা ঘড়িতে: 3-ঘণ্টার ঘড়িতে কাঁটা 2-এ (x ≡ 2 mod 3), 5-ঘণ্টার ঘড়িতে কাঁটা 3-এ (x ≡ 3 mod 5)। আসল x কত? দুই ঘড়ি coprime হলে 0 থেকে 14-এর মধ্যে উত্তর ঠিক একটাই:
x ≡ 2 (mod 3): 2, 5, 8, 11, 14, ...
x ≡ 3 (mod 5): 3, 8, 13, ...
▲
মিলল x = 8 — আর পরেরটা 8 + 15, তার পরে 8 + 30...
উত্তর: x ≡ 8 (mod 15)
Merge-এর সাধারণ নিয়ম: x = a1 + m1·t ধরে দ্বিতীয় congruence-এ বসাও — a1 + m1·t ≡ a2 (mod m2) থেকে t বের করতে লাগে m1-এর inverse mod m2 — আবার সেই extended Euclid!
def crt(a1, m1, a2, m2): # ধরছি gcd(m1, m2) = 1
g, p, q = ext_gcd(m1, m2) # m1·p ≡ 1 (mod m2)
t = ((a2 - a1) * p) % m2
return (a1 + m1 * t) % (m1 * m2), m1 * m2
print(crt(2, 3, 3, 5)) # (8, 15) ✓
দুটোর বেশি congruence? জোড়ায় জোড়ায় merge করতে থাকো — প্রতি merge-এ modulus গুণ হয়ে বড় হয়।
4. Totient advanced আর Mobius — multiplicative দুই ভাই¶
φ(n) = n-এর চেয়ে ছোট বা সমান কতগুলো সংখ্যা n-এর সাথে coprime। Level 2-এ গুনেছিলে; এবার formula — n-এর prime factorization p1^e1 · p2^e2 · ... হলে:
def phi(n):
result = n
p = 2
while p * p <= n:
if n % p == 0:
while n % p == 0:
n //= p
result -= result // p # result *= (1 − 1/p) — integer-এ
p += 1
if n > 1:
result -= result // n
return result
print(phi(36)) # 36·(1/2)·(2/3) = 12
φ multiplicative: gcd(a, b) = 1 হলে φ(a·b) = φ(a)·φ(b)। যেমন φ(36) = φ(4)·φ(9) = 2·6 = 12 ✓। প্রতিটা prime নিজের ভাগটা স্বাধীনভাবে কাটে — এটাই multiplicative হওয়ার মানে।
Mobius function μ(n) — inclusion-exclusion-এর sign-গুলোকে একটা function-এ ভরে দেওয়া:
μ(n) = 1 যদি n = 1
(−1)^k যদি n হয় k-টা ভিন্ন prime-এর গুণফল (square-free)
0 যদি কোনো prime² n-কে ভাগ করে
n : 1 2 3 4 5 6 7 8 9 10 11 12
μ : 1 −1 −1 0 −1 1 −1 0 0 1 −1 0
কোথায় জাদু? Level 3-এর inclusion-exclusion (problem 044) মনে করো — "1 থেকে N-এ কতগুলো সংখ্যা 2 বা 3 দিয়ে ভাগ যায়" গুনতে যোগ-বিয়োগের নাচ। μ সেই নাচের choreography: "square-free d-এর জন্য μ(d) · ⌊N/d⌋ যোগ করো" — sign গুলো μ নিজেই দিয়ে দেয়। "ঠিক gcd = 1" বা "square-free count" ঘরানার গোনায় এটা one-liner বানিয়ে দেয়; বিস্তারিত problem 129-এ।
5. Matrix exponentiation — recurrence-এর time machine¶
Binary exponentiation (level 2, problem 029) মনে আছে — a^n O(log n)-এ, n-কে অর্ধেক করতে করতে। এবার একই খেলা matrix দিয়ে। কেন? কারণ linear recurrence-কে matrix-এ লেখা যায়। Fibonacci:
এক ধাপ verify: [[1,1],[1,0]] · [F(2), F(1)] = [F(2)+F(1), F(2)] = [F(3), F(2)] ✓ — matrix-টা শুধু "যোগ করে এক ঘর সরাও" নিয়মটা ধরে রেখেছে। এখন matrix-এর power binary exponentiation-এ:
def mat_mult(A, B, mod):
return [[(A[0][0]*B[0][0] + A[0][1]*B[1][0]) % mod, (A[0][0]*B[0][1] + A[0][1]*B[1][1]) % mod],
[(A[1][0]*B[0][0] + A[1][1]*B[1][0]) % mod, (A[1][0]*B[0][1] + A[1][1]*B[1][1]) % mod]]
def mat_pow(M, n, mod):
R = [[1, 0], [0, 1]] # identity = matrix জগতের "1"
while n:
if n & 1:
R = mat_mult(R, M, mod)
M = mat_mult(M, M, mod)
n >>= 1
return R
def fib(n, mod=10**9 + 7):
if n == 0: return 0
return mat_pow([[1, 1], [1, 0]], n, mod)[0][1]
print(fib(10)) # 55 ✓
F(10¹⁸) mod p — মাত্র ~60টা matrix গুণে! যেকোনো linear recurrence (a(n) = c1·a(n−1) + c2·a(n−2) + ...) একইভাবে matrix-এ ভরা যায় — এটা একটা সাধারণ ছাঁচ, Fibonacci তার সবচেয়ে বিখ্যাত খদ্দের মাত্র।
Fast doubling — Fibonacci-র জন্য matrix-এরও shortcut। Matrix identity খুলে দুটো formula বেরোয়:
k থেকে সোজা 2k-তে লাফ — recursion-এ n-এর bit ধরে ধরে নামা, matrix-এর 8টা গুণের বদলে 3টা। ছোট কিন্তু elegant; problem 133-এ পুরোটা।
6. Modular inverse — general রাস্তা¶
Level 2-এ Fermat দিয়ে inverse শিখেছিলে: a^(p−2) mod p — কিন্তু সেটা শুধু prime p-তে। General নিয়ম: a-এর inverse mod m আছে ⟺ gcd(a, m) = 1, আর সেটা দেয় extended Euclid:
def mod_inverse(a, m):
g, x, _ = ext_gcd(a, m)
if g != 1:
return None # inverse-এর অস্তিত্বই নেই
return x % m
print(mod_inverse(3, 10)) # 7, কারণ 3·7 = 21 ≡ 1 (mod 10) — 10 prime না, Fermat অচল!
কেন কাজ করে: a·x + m·y = 1 মানে a·x ≡ 1 (mod m) — m-এর গুণিতক mod-এ অদৃশ্য। দুই পথের তুলনা: Fermat ছোট আর মুখস্থ-সহজ কিন্তু prime-বন্দি; ext-gcd সর্বত্রগামী।
7. Miller-Rabin — সাক্ষী জেরা করে prime ধরা¶
Trial division O(√n) — n = 10¹⁸ হলে 10⁹ ধাপ, অচল। Miller-Rabin-এর idea: Fermat's little theorem বলে p prime হলে a^(p−1) ≡ 1 (mod p) — কিন্তু উল্টোটা সবসময় সত্য না (Carmichael number-রা ফাঁকি দেয়, যেমন 561)। তাই আরো ধারালো জেরা: n − 1 = 2^s · d (d বিজোড়) লিখে দেখো a^d, a^(2d), a^(4d), ... ধারাটা কীভাবে 1-এ পৌঁছায়। Prime হলে 1-এ পৌঁছানোর ঠিক আগের পদটা n − 1 হতেই হবে (কারণ prime modulus-এ 1-এর square root শুধু ±1)। নিয়ম ভাঙলে a হলো witness — n composite, প্রমাণিত।
def miller_rabin(n):
if n < 2: return False
for p in [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37]:
if n % p == 0:
return n == p
d, s = n - 1, 0
while d % 2 == 0:
d //= 2; s += 1
for a in [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37]: # এই base গুলোতে 64-bit deterministic
x = pow(a, d, n)
if x in (1, n - 1):
continue
for _ in range(s - 1):
x = x * x % n
if x == n - 1:
break
else:
return False # a হলো witness — composite ধরা পড়ল
return True
print(miller_rabin(561)) # False — Carmichael-ও পার পায় না!
Random base-এ প্রতি জেরায় ভুলের সম্ভাবনা ≤ 1/4; কিন্তু মজার fact — উপরের নির্দিষ্ট 12টা base-এ 64-bit পর্যন্ত উত্তর deterministic (প্রমাণিত)। তাই contest-এ এটাই যথেষ্ট।
8. Pollard rho — birthday paradox দিয়ে factor শিকার¶
বড় composite n-এর একটা factor চাই। Idea: f(x) = (x² + c) mod n দিয়ে একটা pseudo-random ধারা বানাও। n-এর কোনো (অজানা) factor p-এর চোখে এই ধারা mod p-তে ঘুরছে — আর mod p-র জগৎ ছোট (~p টা মান) বলে birthday paradox অনুযায়ী মোটে ~√p ধাপেই দুটো মান mod p-তে মিলে যায় (collision)। তখন gcd(|x − y|, n) ফাঁস করে দেয় p-কে! Cycle ধরতে Floyd-এর কচ্ছপ-খরগোশ:
import math, random
def pollard_rho(n):
if n % 2 == 0: return 2
while True:
x = y = random.randrange(2, n)
c = random.randrange(1, n)
d = 1
while d == 1:
x = (x * x + c) % n # কচ্ছপ — এক ধাপ
y = (y * y + c) % n
y = (y * y + c) % n # খরগোশ — দুই ধাপ
d = math.gcd(abs(x - y), n)
if d != n:
return d # একটা nontrivial factor!
n = 8051 # = 83 × 97
print(pollard_rho(n)) # 83 বা 97
~√p মানে n-এর জন্য O(n^(1/4)) — trial division-এর O(√n)-এর তুলনায় বিশাল লাফ। মনে রাখো: rho চালানোর আগে Miller-Rabin দিয়ে n composite নিশ্চিত করো, নইলে prime-এ অনন্ত ঘুরবে।
9. Sieve variants — এক ছাঁকনির তিন রূপ¶
SPF (smallest prime factor) sieve: প্রতিটা সংখ্যার ক্ষুদ্রতম prime factor টুকে রাখো — তারপর যেকোনো সংখ্যার factorization O(log n)-এ, শুধু spf ধরে ধরে ভাগ:
def build_spf(N):
spf = list(range(N + 1))
for i in range(2, int(N**0.5) + 1):
if spf[i] == i: # i prime
for j in range(i * i, N + 1, i):
if spf[j] == j:
spf[j] = i
return spf
spf = build_spf(100)
n, factors = 60, []
while n > 1:
factors.append(spf[n]); n //= spf[n]
print(factors) # [2, 2, 3, 5] — 60-এর factorization, ভাগ মোটে 4টা
Linear sieve: সাধারণ sieve-এ 12 কাটা পড়ে 2 আর 3 দুজনের হাতে — দ্বৈত কাজ। Linear sieve-এর নিয়ম: প্রতিটা composite-কে শুধু তার SPF দিয়েই কাটো — প্রতিটা সংখ্যা ঠিক একবার, তাই O(n)। Bonus: একই pass-এ φ বা μ-ও বানানো যায় (multiplicative বলে)।
Segmented sieve: "[L, R]-এ prime কোনগুলো, যেখানে R = 10¹², কিন্তু R − L ≤ 10⁶" — পুরো 10¹² সাইজের array অসম্ভব। সমাধান: √R পর্যন্ত prime গুলো সাধারণ sieve-এ বানাও, তারপর শুধু [L, R] window-র ছোট array-তে সেই prime-দের গুণিতক কাটো। Window-এ p-এর প্রথম গুণিতক: max(p*p, ((L + p − 1) // p) * p) — এই ceiling হিসাবটাই পুরো problem 138-এর হৃদয়।
এক নজরে¶
technique কীসের চালাক রূপ complexity-র লাফ
────────────────────────────────────────────────────────────────────────
extended Euclid Euclid + পথ মনে রাখা একই O(log n), বেশি তথ্য
CRT ext-gcd-এর প্রয়োগ দুই modulus → এক
Mobius inclusion-exclusion encode নাচ → এক sum
matrix exponentiation binary expo, সংখ্যা → matrix O(n) → O(log n)
Miller-Rabin Fermat test + square-root জেরা O(√n) → O(log³ n)
Pollard rho birthday paradox + cycle O(√n) → O(n^(1/4))
SPF / linear sieve sieve + কে কাটল মনে রাখা factorization O(log n)
segmented sieve sieve-কে window-তে নামানো memory: O(R) → O(R−L)
পরের ধাপ: visualization-ideas.md-তে back-substitution-এর সিঁড়ি আর CRT-র ঘড়ি এঁকে দেখো — তারপর problem 125 থেকে শুরু।