Divide and conquer — Strassen, closest pair আর Karatsuba

১ · ভাগ করুন, সমাধান করুন, জোড়া দিন — কৌশলটি

মডিউল ৭-এ আমরা একটি নির্দিষ্ট ডেটা স্ট্রাকচার নয়, বরং অ্যালগরিদম ডিজাইনের কৌশল শিখব — নতুন সমস্যা দেখলে কোন দিক থেকে ভাবতে শুরু করবেন। প্রথম কৌশলটি আপনি আগেই ব্যবহার করেছেন: merge sort (পর্ব ২১), quick sort (পর্ব ২২) আর binary search (পর্ব ১৮) — তিনটিই divide and conquer. এই পর্বে কৌশলটিকে নাম ধরে আলাদা করে দেখব, তারপর এমন কয়েকটি সমস্যায় প্রয়োগ করব যেখানে এর শক্তি সত্যিই চমকে দেয়।

প্রতিটি divide and conquer অ্যালগরিদমে তিনটি অংশ থাকে:

• Divide: ইনপুটকে a টি উপ-সমস্যায় ভাগ, প্রতিটি সাধারণত মূলের n/b আকারের। ভাগটা সস্তা হতে পারে (merge sort-এ শুধু মাঝের ইনডেক্স), আবার দামিও হতে পারে (quick sort-এর partition).
• Conquer: প্রতিটি উপ-সমস্যা রিকার্শনে একই অ্যালগরিদম দিয়ে সমাধান। আকার যথেষ্ট ছোট হলে (base case) সরাসরি উত্তর।
• Combine: ছোট উত্তরগুলো থেকে বড় উত্তর তৈরি। এর খরচই সাধারণত অ্যালগরিদমের গতি ঠিক করে।

SOLVE(P)                          // P: a problem of size n
    if n ≤ n0: return DIRECT(P)     // base case, Θ(1) for constant n0
    P1, ..., Pa = DIVIDE(P)         // a pieces, each of size about n/b
    for i = 1 to a:
        Ri = SOLVE(Pi)              // a recursive calls: a · T(n/b)
    return COMBINE(R1, ..., Ra)     // divide + combine together cost f(n)

এই কাঠামো থেকে সরাসরি recurrence লেখা যায়: T(n) = aT(n/b) + f(n), যেখানে f(n) = divide আর combine-এর মোট খরচ। তারপর পর্ব ৯-এর master theorem বা পর্ব ৮-এর recursion tree দিয়ে উত্তর। এই পর্বের অ্যালগরিদমগুলোর একটি অগ্রিম তালিকা:

২ · Maximum subarray — মাঝখান পার হওয়া অংশই আসল কাজ

পর্ব ১০-এ দেখেছেন maximum subarray সমস্যা: একটি দোকানের দৈনিক লাভ-ক্ষতির তালিকায় টানা কোন কয়েক দিনে মোট লাভ সবচেয়ে বেশি? সব (l, r) জোড়া চেষ্টা করলে Θ(n²). Divide and conquer দিয়ে ভাবুন: অ্যারেটি মাঝখানে ভাগ করলে সেরা অংশটি তিন জায়গার একটিতে থাকবেই —

• পুরোটাই বাঁ অর্ধেকে → রিকার্শন খুঁজে দেবে;
• পুরোটাই ডান অর্ধেকে → রিকার্শন খুঁজে দেবে;
• মাঝের রেখা পার হয়ে — বাঁ অর্ধেকের শেষ ঘর a[mid−1] আর ডান অর্ধেকের প্রথম ঘর a[mid] দুটিই এর মধ্যে। এটি রিকার্শন দেখে না, combine ধাপকে খুঁজতে হবে।

তৃতীয় ক্ষেত্রের চাবি: পার হওয়া যেকোনো অংশ = বাঁ অর্ধেকের একটি suffix (শেষ দিকের টুকরো) + ডান অর্ধেকের একটি prefix (শুরুর দিকের টুকরো). দুটি স্বাধীনভাবে বাছা যায়, তাই সেরা suffix আর সেরা prefix আলাদা করে খুঁজে যোগ করলেই সেরা পার-হওয়া অংশ। মাঝ থেকে বাঁয়ে একটি লুপ, ডানে একটি লুপ — মোট Θ(n).

MAXSUB(a, lo, hi)                     // 0-based, half-open a[lo..hi), hi - lo ≥ 1
    if hi - lo == 1: return a[lo]
    mid = (lo + hi) / 2
    best = max(MAXSUB(a, lo, mid), MAXSUB(a, mid, hi))
    s = 0; bestL = -∞
    for i = mid-1 down to lo: s = s + a[i]; bestL = max(bestL, s)   // best suffix
    s = 0; bestR = -∞
    for j = mid to hi-1:      s = s + a[j]; bestR = max(bestR, s)   // best prefix
    return max(best, bestL + bestR)

৮ দিনের হিসাব (হাজার টাকায়) [−2, 1, −3, 4, −1, 2, 1, −5] নিয়ে ধাপে ধাপে দেখুন:

index :   0   1   2   3 │  4   5   6   7
a     :  -2   1  -3   4 │ -1   2   1  -5
          left half     │   right half

[-2 1] │ [-3 4]
left 1 , right 4 , cross = 1 + 1 = 2
best of left half = 4    ( [4] )

[-1 2] │ [1 -5]
left 2 , right 1 , cross = 2 + 1 = 3
best of right half = 3   ( [2 1] )

i      :  3    2    1    0
a[i]   :  4   -3    1   -2
s      :  4    1    2    0
bestL  = 4   ( a[3..3] )

j      :  4    5    6    7
a[j]   : -1    2    1   -5
s      : -1    1    2   -3
bestR  = 2   ( a[4..6] )

left  = 4
right = 3
cross = bestL + bestR = 4 + 2 = 6   ->  a[3..6] = [4 -1 2 1]
answer = max(4, 3, 6) = 6

কেন সঠিক, আর কত দ্রুত

সঠিকতা (induction, আকারের ওপর): আকার ১-এ উত্তর একমাত্র ঘরটিই। ধরুন ছোট সব আকারে MAXSUB সঠিক। সেরা অংশ [l, r] হয় r < mid (পুরো বাঁয়ে — হাইপোথিসিসে বাঁয়ের রিকার্শন অন্তত এর সমান পায়), নয়তো l ≥ mid (একইভাবে ডানে), নয়তো l < mid ≤ r. শেষ ক্ষেত্রে এর যোগফল = sum(a[l..mid−1]) + sum(a[mid..r]) ≤ bestL + bestR, আর bestL + bestR নিজেও একটি বৈধ পার-হওয়া অংশের যোগফল। তাই তিন প্রার্থীর সর্বোচ্চই উত্তর।

সময়: দুটি অর্ধেক-আকারের ডাক আর দুটি লুপ মিলিয়ে Θ(n) — T(n) = 2T(n/2) + Θ(n), master theorem-এর case 2: Θ(n log n). জায়গা: কোনো বাড়তি অ্যারে নেই, শুধু রিকার্শনের stack — Θ(log n).

৩ · Combine ধাপের কারিগরি — inversion গোনা আর দ্রুত ঘাত

আরও দুটি পরিচিত উদাহরণে combine ধাপকে আলাদা চোখে দেখি।

Inversion গোনা — merge-এর ভেতরে লুকানো হিসাব

অ্যারেতে i < j অথচ a[i] > a[j] — এমন জোড়াকে বলে inversion. পর্ব ২১-এর শেষ দিকে এটি merge sort দিয়ে Θ(n log n)-এ গুনেছেন। Divide and conquer-এর ভাষায় পুরো ব্যাপারটা এক লাইনে: মোট inversion = বাঁয়ের ভেতরের + ডানের ভেতরের + পার-হওয়া (একটি বাঁয়ে, একটি ডানে). প্রথম দুটি রিকার্শন গোনে। পার-হওয়াগুলো সরাসরি গুনলে Θ(n²) লাগত — কিন্তু দুই অর্ধেক যদি আগেই সাজানো থাকে, তাহলে merge করার সময় ডান থেকে a[j] নেওয়ার মুহূর্তে বাঁয়ে অপেক্ষায় থাকা সবাই (mid − i জন) তার চেয়ে বড় — এক যোগে গোনা শেষ। অর্থাৎ রিকার্শন শুধু উত্তর নয়, combine-কে সহজ করার মতো বাড়তি তথ্যও (সাজানো ক্রম) ফেরত দেয়। এই কৌশলটি — “সমস্যাকে একটু শক্তিশালী করে নিন, যাতে জোড়া দেওয়া সহজ হয়” — closest pair-এও আবার দেখবেন, যেখানে রিকার্শন বিন্দুগুলো y অনুযায়ী সাজিয়ে ফেরত দেয়।

দ্রুত ঘাত — একটি উপ-সমস্যাই যথেষ্ট

3¹³ বের করতে সরলভাবে ১২টি গুণ লাগে। কিন্তু x^e = (x^e/2)² (e জোড় হলে) আর x · (x^(e−1)/2)² (e বিজোড় হলে). এখানে a = 1: একটিমাত্র উপ-সমস্যা, আর তার উত্তরকে নিজের সাথে গুণ করাই combine. পর্ব ৭-এ এটি binary exponentiation নামে দেখেছেন। 3¹³-এর হিসাব:

Recurrence T(e) = T(e/2) + Θ(1), তাই Θ(log e) গুণ — প্রতিটি ধাপে বড়জোর ২টি। e = 10¹⁸ হলেও মাত্র ৬০টির মতো ধাপ। প্রতিটি গুণের পর mod নিলে (modular arithmetic) সংখ্যা ছোট থাকে — কনটেস্টে “উত্তর 10⁹ + 7 দিয়ে mod করুন” ধরনের সমস্যার এটি রোজকার হাতিয়ার, আর RSA-র মতো cryptography-র ভিত্তিও। একই কৌশল matrix-এ খাটে: 2×2 matrix-এর ঘাত দিয়ে Fibonacci-র n-তম সংখ্যা Θ(log n) গুণে।

৪ · Karatsuba — চারটির বদলে তিনটি গুণ

স্কুলে শেখা গুণে দুটি n-অঙ্কের সংখ্যা গুণ করতে প্রতিটি অঙ্ককে প্রতিটি অঙ্কের সাথে গুণ করতে হয় — n²টি এক-অঙ্কের গুণ। একে বলে schoolbook multiplication. ১০০ অঙ্কে ১০,০০০টি গুণ, আর cryptography-তে ব্যবহৃত হাজার অঙ্কের সংখ্যায় ১০ লাখ। এর চেয়ে কম কি সম্ভব?

কৌশলটি দেখুন। দুটি n-অঙ্কের সংখ্যাকে অর্ধেক করি, B = 10^n/2:

তাহলে xy = x₁y₁·B² + (x₁y₀ + x₀y₁)·B + x₀y₀. সরাসরি করলে চারটি অর্ধেক-আকারের গুণ: T(n) = 4T(n/2) + Θ(n), যা master theorem-এ Θ(n^{log₂ 4}) = Θ(n²) — কোনো লাভ নেই। Karatsuba-র চালাকি: মাঝের পদটির জন্য আলাদা দুটি গুণ লাগে না।

কারণ (x₁ + x₀)(y₁ + y₀) = x₁y₁ + x₁y₀ + x₀y₁ + x₀y₀ — এখান থেকে z₂ আর z₀ বাদ দিলে ঠিক মাঝের পদটি থাকে। মোট গুণ তিনটি; বিনিময়ে কয়েকটি বাড়তি যোগ-বিয়োগ, যার খরচ Θ(n). শেষে xy = z₂·B² + z₁·B + z₀ — B-এর গুণ মানে শুধু ডানে শূন্য বসানো (shift), আসল গুণ নয়।

বিশ্লেষণ

T(n) = 3T(n/2) + Θ(n). Watershed n^{log₂ 3} ≈ n^1.585, আর f(n) = n তার চেয়ে polynomially ছোট — case 1: T(n) = Θ(n^{log₂ 3}) ≈ Θ(n^1.585). Recursion tree-তে দেখলে: তলা j-তে 3^jটি গুণ, প্রতিটি n/2^j আকারের, তলার মোট খরচ (3/2)^j·n — প্রতি তলায় দেড় গুণ বাড়ে, তাই পাতাগুলো জেতে, আর পাতা ঠিক 3^{log₂ n} = n^{log₂ 3}টি। ১০২৪ অঙ্কে স্কুলের পদ্ধতি প্রায় ১০ লাখ এক-অঙ্কের গুণ করে, Karatsuba করে 3¹⁰ = ৫৯,০৪৯টি। জায়গা: প্রতি তলায় অস্থায়ী অঙ্কের তালিকা, মূল থেকে পাতা পর্যন্ত পথে n + n/2 + n/4 + … = Θ(n).

নিচের প্রোগ্রাম দশমিক অঙ্কের স্ট্রিং হিসেবে দেওয়া দুটি বড় সংখ্যা Karatsuba দিয়ে গুণ করে, এক-অঙ্কের গুণ গোনে, আর ফল যাচাই করে — C++-এ স্কুলের O(n²) পদ্ধতির সাথে, Python-এ ভাষার নিজস্ব বড়-সংখ্যার int-এর সাথে। অঙ্কগুলো উল্টো ক্রমে (little-endian) রাখা, আর দৈর্ঘ্য ২-এর ঘাতে পূরণ করা হয়েছে, যাতে প্রতিবার ঠিক অর্ধেক হয়। Carry শুধু একেবারে শেষে।

#include <cstdio>
#include <string>
#include <vector>
using namespace std;
using V = vector<long long>;
long long calls = 0;                            // one-digit multiplications done
V kara(V x, V y) {                              // little-endian digits, size = 2^k
    size_t n = x.size(), m = n / 2;
    if (n == 1) { ++calls; return {x[0] * y[0]}; }
    V xh(x.begin() + m, x.end()), yh(y.begin() + m, y.end());
    x.resize(m); y.resize(m);                   // x, y now hold the low halves
    V z0 = kara(x, y), z2 = kara(xh, yh);       // low*low, high*high
    for (size_t i = 0; i < m; ++i) x[i] += xh[i], y[i] += yh[i];
    V s = kara(x, y), r(2 * n - 1, 0);          // (low+high)(low+high)
    for (size_t i = 0; i + 1 < 2 * m; ++i) {
        r[i] += z0[i], r[i + 2 * m] += z2[i];
        r[i + m] += s[i] - z0[i] - z2[i];       // middle = s - z0 - z2
    }
    return r;
}
string digits(V c) {                            // carry coefficients into a string
    string out; long long carry = 0; c.resize(c.size() + 2, 0);
    for (long long v : c) { carry += v; out += char('0' + carry % 10); carry /= 10; }
    while (out.size() > 1 && out.back() == '0') out.pop_back();
    return string(out.rbegin(), out.rend());
}
int main() {
    for (string a : {"1234", "31415926535897932384626433832795"}) {
        string b = a.size() == 4 ? "5678" : "27182818284590452353602874713527";
        size_t n = 1; while (n < max(a.size(), b.size())) n *= 2;
        V x(n, 0), y(n, 0), plain(2 * n, 0);
        for (size_t i = 0; i < a.size(); ++i) x[i] = a[a.size() - 1 - i] - '0';
        for (size_t i = 0; i < b.size(); ++i) y[i] = b[b.size() - 1 - i] - '0';
        for (size_t i = 0; i < n; ++i)          // school method, O(n^2), as a check
            for (size_t j = 0; j < n; ++j) plain[i + j] += x[i] * y[j];
        calls = 0; string p = digits(kara(x, y));
        printf("%s x %s\n= %s\n", a.c_str(), b.c_str(), p.c_str());
        printf("digit mults: karatsuba %lld, school %zu, check %s\n", calls,
               a.size() * b.size(), p == digits(plain) ? "ok" : "FAIL");
    }
}

calls = 0                                   # one-digit multiplications done

def kara(x, y):                             # little-endian digits, len = 2^k
    global calls
    n = len(x)
    if n == 1:
        calls += 1
        return [x[0] * y[0]]
    m = n // 2
    z0 = kara(x[:m], y[:m])                 # low  * low
    z2 = kara(x[m:], y[m:])                 # high * high
    s = kara([a + b for a, b in zip(x[:m], x[m:])],
             [a + b for a, b in zip(y[:m], y[m:])])   # (low+high)(low+high)
    r = [0] * (2 * n - 1)
    for i in range(2 * m - 1):
        r[i] += z0[i]
        r[i + 2 * m] += z2[i]
        r[i + m] += s[i] - z0[i] - z2[i]    # middle = s - z0 - z2
    return r

def multiply(a, b):                         # a, b: decimal strings
    n = 1 << (max(len(a), len(b)) - 1).bit_length()     # next power of two
    x, y = ([int(d) for d in reversed(s)] + [0] * (n - len(s)) for s in (a, b))
    out, carry = [], 0
    for v in kara(x, y) + [0, 0]:           # carry coefficients into digits
        carry += v
        out.append(carry % 10)
        carry //= 10
    return ''.join(map(str, reversed(out))).lstrip('0') or '0'

for a, b in [("1234", "5678"), ("31415926535897932384626433832795",
                                "27182818284590452353602874713527")]:
    calls = 0
    p = multiply(a, b)
    print(f"{a} x {b}\n= {p}")
    print(f"digit mults: karatsuba {calls}, school {len(a) * len(b)},",
          "check ok" if p == str(int(a) * int(b)) else "check FAIL")

1234 x 5678
= 7006652
digit mults: karatsuba 9, school 16, check ok
31415926535897932384626433832795 x 27182818284590452353602874713527
= 853973422267356706546355086954668447174445893592492906242717965
digit mults: karatsuba 243, school 1024, check ok

৫ · Strassen — matrix গুণে আটটির বদলে সাতটি

দুটি n×n matrix গুণ করার সাধারণ পদ্ধতিতে তিনটি নেস্টেড লুপ: C[i][j] = Σ_k A[i][k]·B[k][j], মোট n³টি গুণ। Divide and conquer দিয়ে প্রতিটি matrix-কে চারটি n/2 × n/2 টুকরায় (block) ভাঙা যায়:

এখানে ৮টি অর্ধেক-আকারের matrix গুণ আর ৪টি matrix যোগ (প্রতিটি Θ(n²)): T(n) = 8T(n/2) + Θ(n²) = Θ(n^{log₂ 8}) = Θ(n³). আবারও কোনো লাভ নেই — Karatsuba-র মতোই, সরল ভাগে ঘাত বদলায় না। ১৯৬৯ সালে জার্মান গণিতবিদ Volker Strassen “Gaussian elimination is not optimal” প্রবন্ধে দেখান যে ৭টি গুণেই চলে — এটিই Strassen-এর অ্যালগরিদম. তাঁর সাতটি গুণফল:

M1 = (A11 + A22) · (B11 + B22)
M2 = (A21 + A22) · B11
M3 = A11 · (B12 - B22)
M4 = A22 · (B21 - B11)
M5 = (A11 + A12) · B22
M6 = (A21 - A11) · (B11 + B12)
M7 = (A12 - A22) · (B21 + B22)

C11 = M1 + M4 - M5 + M7        C12 = M3 + M5
C21 = M2 + M4                  C22 = M1 - M2 + M3 + M6

মোট ৭টি গুণ আর ১৮টি যোগ-বিয়োগ (M-এ ১০টি, C-তে ৮টি). Matrix গুণ commutative নয় (AB ≠ BA সাধারণত), আর সূত্রগুলোতে কোথাও গুণের ক্রম উল্টানো হয়নি — তাই A_ij নিজেরাই matrix হলেও সূত্র খাটে, আর রিকার্শন চলে। ছোট সংখ্যায় হাতে মিলিয়ে দেখি:

A = | 1  2 |    B = | 5  6 |
    | 3  4 |        | 7  8 |

C = A·B = | 1·5+2·7  1·6+2·8 |  =  | 19  22 |
          | 3·5+4·7  3·6+4·8 |     | 43  50 |

M1 = (1 + 4)(5 + 8) = 5 · 13 =  65
M2 = (3 + 4) · 5    = 7 · 5  =  35
M3 = 1 · (6 - 8)    = 1 · -2 =  -2
M4 = 4 · (7 - 5)    = 4 · 2  =   8

M5 = (1 + 2) · 8    = 3 · 8  =  24
M6 = (3 - 1)(5 + 6) = 2 · 11 =  22
M7 = (2 - 4)(7 + 8) = -2 · 15 = -30

C11 = M1 + M4 - M5 + M7 = 65 + 8 - 24 - 30 = 19
C12 = M3 + M5           = -2 + 24          = 22
C21 = M2 + M4           = 35 + 8           = 43
C22 = M1 - M2 + M3 + M6 = 65 - 35 - 2 + 22 = 50

n = 2      :  7 multiplications      (plain: 8)
n = 4      :  7^2  =     49            (plain: 64)
n = 1024   :  7^10 ≈ 2.8 · 10^8        (plain: 1024^3 ≈ 10.7 · 10^8)
general    :  7^(log2 n) = n^(log2 7) ≈ n^2.807

বাস্তবে Strassen — crossover আর সতর্কতা

তত্ত্বে Strassen জেতে, কিন্তু Θ লুকিয়ে রাখে ধ্রুবক। ১৮টি matrix যোগ-বিয়োগ, অস্থায়ী matrix-এর জন্য বাড়তি মেমরি আর cache-এর কম-বান্ধব অ্যাক্সেস — ছোট matrix-এ এগুলোই বেশি খরচ করে। তাই বাস্তব বাস্তবায়ন hybrid: আকার একটি crossover point-এর নিচে নামলেই সাধারণ পদ্ধতি। এই সীমা হার্ডওয়্যার আর বাস্তবায়নের ওপর এতটাই নির্ভর করে যে ভিন্ন গবেষণায় কয়েক ডজন থেকে কয়েক হাজার পর্যন্ত মান পাওয়া গেছে — মেপে ঠিক করতে হয়। দ্বিতীয় সমস্যা দশমিক সংখ্যায় (floating point): বেশি বিয়োগের কারণে rounding ভুল সাধারণ পদ্ধতির চেয়ে কিছুটা বেশি জমতে পারে, তাই বৈজ্ঞানিক গণনার অনেক লাইব্রেরি (BLAS ধাঁচের) মূলত খুব যত্নে অপ্টিমাইজ করা Θ(n³) পদ্ধতিই চালায়। পূর্ণসংখ্যা বা modular matrix-এ এই সমস্যা নেই — সেখানে Strassen বেশি স্বচ্ছন্দ।

Strassen-এর পর শুরু হয় matrix গুণের ঘাত ω কমানোর দীর্ঘ প্রতিযোগিতা — ω হলো সবচেয়ে ছোট সংখ্যা যাতে n×n matrix গুণ O(n^ω+ε) সময়ে করা যায়। জানা আছে 2 ≤ ω (প্রতিটি ঘর অন্তত একবার লিখতে হয়) আর ω ≤ log₂ 7 ≈ 2.807 (Strassen). ১৯৯০ সালে Coppersmith ও Winograd ঘাতটিকে প্রায় ২.৩৭৬-এ নামান; ২০২০-এর দশকের গবেষণাপত্রগুলোতে ওপরের সীমা প্রায় ২.৩৭১-এ এসেছে, আর সামান্য উন্নতি এখনো চলছে। ω আসলে ২ কি না — এটি অ্যালগরিদম তত্ত্বের বড় খোলা প্রশ্নগুলোর একটি। তবে এই নতুন অ্যালগরিদমগুলোর ধ্রুবক এত বিশাল যে কোনো বাস্তব আকারে কাজে লাগে না — এদের বলে galactic algorithm.

৬ · Closest pair of points — strip-এর যুক্তি

সমতলে n টি বিন্দু দেওয়া আছে — ধরুন শাহজালাল বিমানবন্দরের রাডারে n টি উড়োজাহাজের অবস্থান, বা একটি শহরের মোবাইল টাওয়ার। কোন দুটি সবচেয়ে কাছাকাছি? এটি closest pair of points সমস্যা। Brute force সব n(n−1)/2 জোড়া দেখে — Θ(n²). এক মাত্রায় (সব বিন্দু একটি সরলরেখায়) উত্তর সহজ: সাজিয়ে পাশাপাশি জোড়াগুলো দেখুন, Θ(n log n). দুই মাত্রায় “পাশাপাশি” বলে কিছু নেই — x-এ কাছের বিন্দু y-এ অনেক দূরে হতে পারে। Michael Shamos ও Dan Hoey ১৯৭৫ সালে divide and conquer দিয়ে Θ(n log n) সমাধান দেন।

ধারণা

1. বিন্দুগুলো x অনুযায়ী সাজান। মাঝের বিন্দুর x-এ একটি খাড়া রেখা টেনে বাঁ আর ডান অর্ধেকে ভাগ করুন।
2. রিকার্শনে বাঁয়ের সেরা দূরত্ব δ_L আর ডানের δ_R বের করুন। δ = min(δ_L, δ_R).
3. বাকি থাকে পার-হওয়া জোড়া — একটি বাঁয়ে, একটি ডানে। এমন জোড়া δ-এর চেয়ে কাছে হতে হলে দুটি বিন্দুই রেখা থেকে δ-এর কম দূরে থাকবে। তাই শুধু রেখার দুই পাশে δ চওড়া একটি strip-এর বিন্দুগুলো দেখলেই চলে।
4. Strip-এর বিন্দুগুলো y অনুযায়ী সাজিয়ে প্রতিটির জন্য শুধু ঠিক আগের যে কয়েকটির y-ব্যবধান δ-এর কম, তাদের সাথে দূরত্ব মাপুন।

কিন্তু সবগুলো বিন্দুই তো strip-এ পড়তে পারে! তাহলে strip-এর ভেতরেও কি Θ(n²)? না — আর এর কারণটিই এই অ্যালগরিদমের প্রাণ।

প্রমাণ (pigeonhole নীতি): q-এর নিচে y-ব্যবধান < δ এমন সব বিন্দু একটি δ উঁচু আর 2δ চওড়া আয়তক্ষেত্রে পড়ে (রেখার বাঁয়ে δ, ডানে δ). এটিকে ৮টি δ/2 × δ/2 বর্গে ভাগ করুন — রেখার প্রতিটি পাশে ৪টি। একটি বর্গের কর্ণ δ/√2 < δ, তাই একটি বর্গে দুটি বিন্দু থাকলে তাদের দূরত্ব δ-এর কম হতো। কিন্তু প্রতিটি বর্গ পুরোপুরি এক পাশে, আর এক পাশের ভেতরের যেকোনো দুটি বিন্দুর দূরত্ব অন্তত δ (δ_L, δ_R-এর সংজ্ঞায়). তাই প্রতিটি বর্গে বড়জোর একটি বিন্দু, আয়তক্ষেত্রে বড়জোর ৮টি — q নিজে তার একটি, বাকি বড়জোর ৭টি। (একই x বা y-এর বিন্দুকে কোন বর্গে ধরা হবে, সেটি সীমারেখায় একটু যত্ন নিয়ে ঠিক করলে যুক্তি একই থাকে। আরও সূক্ষ্ম হিসাবে সংখ্যাটি আরও কমানো যায়, কিন্তু ধ্রুবক হওয়াই আসল কথা।)

৮টি বিন্দু নিয়ে দেখি: A(1,2), B(2,7), C(3,4), D(6,6), E(7,5), F(9,1), G(10,8), H(12,4). নিচের ছবিতে x অক্ষে শেষ তিনটি ঘর ১০, ১১, ১২.

8  . . . . . . .│. . . G . .
7  . . B . . . .│. . . . . .
6  . . . . . . D│. . . . . .
5  . . . . . . .│E . . . . .
4  . . . C . . .│. . . . . H
3  . . . . . . .│. . . . . .
2  . A . . . . .│. . . . . .
1  . . . . . . .│. . F . . .
0  . . . . . . .│. . . . . .
                       1 1 1
   0 1 2 3 4 5 6 7 8 9 0 1 2

A-C : 2^2 + 2^2 =  8   <- best, δL = √8 ≈ 2.83
B-C : 1^2 + 3^2 = 10
C-D : 3^2 + 2^2 = 13
B-D : 4^2 + 1^2 = 17
A-B : 26 ,  A-D : 41

E-G : 3^2 + 3^2 = 18   <- best, δR = √18 ≈ 4.24
F-H : 3^2 + 3^2 = 18
E-F : 2^2 + 4^2 = 20
G-H : 2^2 + 4^2 = 20
δ = min(δL, δR) = √8 ≈ 2.83

8  . . . . .┆. .│. . .┆o . .
7  . . o . .┆. .│. . .┆. . .
6  . . . . .┆. D│. . .┆. . .
5  . . . . .┆. .│E . .┆. . .
4  . . . o .┆. .│. . .┆. . o
3  . . . . .┆. .│. . .┆. . .
2  . o . . .┆. .│. . .┆. . .
1  . . . . .┆. .│. . F┆. . .
0  . . . . .┆. .│. . .┆. . .
                       1 1 1
   0 1 2 3 4 5 6 7 8 9 0 1 2

strip, in increasing y :  F(9,1)  E(7,5)  D(6,6)

F(9,1) : nothing below                  -> 0 checks
E(7,5) : F   dy = 4,  4^2 = 16 >= 8  stop -> 0 checks

D(6,6) : E   dy = 1,  d^2 = 1 + 1 = 2 < 8  -> new best
         F   dy = 5,  25 >= 2  stop
answer : D-E,  d = √2 ≈ 1.41   (crosses the line)

CLOSEST(P)                        // P sorted by x; returns δ and P re-sorted by y
    if |P| ≤ 3: check all pairs; sort P by y; return
    mid = |P| / 2;  line = P[mid].x
    δ = min(CLOSEST(P[0..mid)), CLOSEST(P[mid..|P|)))
    merge the two halves by y                    // Θ(n), as in merge sort
    S = points of P with |x - line| < δ, in y order
    for each q in S:
        for each earlier s in S while q.y - s.y < δ:   // at most 7 of them
            δ = min(δ, dist(q, s))
    return δ

বিশ্লেষণ — কেন n log n, n log² n নয়

প্রথমে একবার x অনুযায়ী সাজানো: Θ(n log n). তারপর প্রতিটি ডাকে দুটি অর্ধেক-ডাক, Θ(n)-এ y-merge, Θ(n)-এ strip বানানো, আর strip-এ প্রতিটি বিন্দুর জন্য বড়জোর ৭টি তুলনা — মোট Θ(n). তাই T(n) = 2T(n/2) + Θ(n) = Θ(n log n), মোট সময়ও Θ(n log n). জায়গা: merge আর strip-এর জন্য Θ(n). লক্ষ করুন — প্রতিটি ডাকে strip-কে নতুন করে y অনুযায়ী সাজালে combine Θ(n log n) হয়ে যায়, আর T(n) = 2T(n/2) + Θ(n log n) = Θ(n log² n) (পর্ব ৯-এর বর্ধিত case 2). Merge sort-এর মতো রিকার্শনকে y-সাজানো তালিকা ফেরত দিতে বলা — অংশ ৩-এর সেই “সমস্যাকে একটু শক্তিশালী করা” — বাড়তি log-টা সরিয়ে দেয়।

প্রোগ্রামে D&C উত্তরকে brute force-এর সাথে মেলানো হয়েছে দুই সেট বিন্দুতে: ডায়াগ্রামের ৮টি বিন্দু, আর একটি সূত্র দিয়ে বানানো ১০০০টি ভিন্ন বিন্দু (x = i² mod 10007, y = i³ mod 10009) — যাতে দুই ভাষায় হুবহু এক ইনপুট থাকে। সাথে গোনা হয়েছে strip-এ একটি বিন্দুর জন্য সর্বোচ্চ কয়টি দূরত্ব মাপতে হয়েছে।

#include <algorithm>
#include <cstdio>
#include <initializer_list>
#include <utility>
#include <vector>
using namespace std;
using ll = long long;
struct P { ll x, y; };
ll d2(P a, P b) { return (a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y); }
bool byY(P a, P b) { return a.y < b.y; }
int most = 0;                                     // most checks for one strip point
ll rec(vector<P>& p, int lo, int hi) {            // [lo, hi): by x in, by y out
    ll best = d2(p[lo], p[lo + 1]);
    if (hi - lo <= 3) {                           // tiny: check every pair
        for (int i = lo; i < hi; ++i)
            for (int j = i + 1; j < hi; ++j) best = min(best, d2(p[i], p[j]));
        stable_sort(p.begin() + lo, p.begin() + hi, byY); return best;
    }
    int mid = (lo + hi) / 2; ll mx = p[mid].x;    // the dividing vertical line
    best = min(rec(p, lo, mid), rec(p, mid, hi));
    inplace_merge(p.begin() + lo, p.begin() + mid, p.begin() + hi, byY);
    vector<P> strip;
    for (int i = lo; i < hi; ++i) {               // in increasing y
        P q = p[i]; int k = 0;
        if ((q.x - mx) * (q.x - mx) >= best) continue;
        for (int j = int(strip.size()) - 1; j >= 0; --j, ++k) {
            if ((q.y - strip[j].y) * (q.y - strip[j].y) >= best) break;   // too far below
            best = min(best, d2(q, strip[j]));
        }
        most = max(most, k); strip.push_back(q);
    }
    return best;
}
int main() {
    vector<P> fixed = {{1, 2}, {2, 7}, {3, 4}, {6, 6}, {7, 5}, {9, 1}, {10, 8}, {12, 4}}, big;
    for (ll i = 1; i <= 1000; ++i) big.push_back({i * i % 10007, i * i * i % 10009});
    for (auto [name, pts] : {pair{"fixed 8", fixed}, pair{"formula 1000", big}}) {
        ll brute = d2(pts[0], pts[1]);
        for (size_t i = 0; i < pts.size(); ++i)
            for (size_t j = i + 1; j < pts.size(); ++j) brute = min(brute, d2(pts[i], pts[j]));
        sort(pts.begin(), pts.end(), [](P a, P b) { return a.x != b.x ? a.x < b.x : a.y < b.y; });
        most = 0; ll dc = rec(pts, 0, int(pts.size()));
        printf("%s: dc d^2 = %lld, brute d^2 = %lld, max strip checks = %d\n", name, dc, brute, most);
    }
}

from heapq import merge

def d2(p, q):                               # squared distance, exact integers
    return (p[0] - q[0]) ** 2 + (p[1] - q[1]) ** 2

def rec(p, lo, hi, st):                     # p[lo:hi]: by x in, by y out
    if hi - lo <= 3:
        best = min(d2(p[i], p[j]) for i in range(lo, hi) for j in range(i + 1, hi))
        p[lo:hi] = sorted(p[lo:hi], key=lambda q: q[1])
        return best
    mid = (lo + hi) // 2
    mx = p[mid][0]                          # the dividing vertical line
    best = min(rec(p, lo, mid, st), rec(p, mid, hi, st))
    p[lo:hi] = list(merge(p[lo:mid], p[mid:hi], key=lambda q: q[1]))
    strip = []
    for q in p[lo:hi]:                      # in increasing y
        if (q[0] - mx) ** 2 >= best:
            continue
        k = 0
        for s in reversed(strip):
            if (q[1] - s[1]) ** 2 >= best:
                break                       # too far below
            best, k = min(best, d2(q, s)), k + 1
        st[0] = max(st[0], k)               # most checks for one strip point
        strip.append(q)
    return best

fixed = [(1, 2), (2, 7), (3, 4), (6, 6), (7, 5), (9, 1), (10, 8), (12, 4)]
big = [(i * i % 10007, i * i * i % 10009) for i in range(1, 1001)]
for name, pts in [("fixed 8", fixed), ("formula 1000", big)]:
    brute = min(d2(p, q) for i, p in enumerate(pts) for q in pts[i + 1:])
    st = [0]
    dc = rec(sorted(pts), 0, len(pts), st)
    print(f"{name}: dc d^2 = {dc}, brute d^2 = {brute}, max strip checks = {st[0]}")

fixed 8: dc d^2 = 2, brute d^2 = 2, max strip checks = 2
formula 1000: dc d^2 = 5, brute d^2 = 5, max strip checks = 2

৭ · Linear-time selection, আর কখন divide and conquer কাজ করে না

k-তম ছোট — median of medians

পর্ব ২২-এ quickselect দেখেছেন: একটি pivot দিয়ে partition, তারপর শুধু যে দিকে k-তম আছে সেদিকে রিকার্শন — গড়ে Θ(n), কিন্তু pivot বারবার খারাপ হলে Θ(n²). ১৯৭৩ সালে Blum, Floyd, Pratt, Rivest ও Tarjan-এর median of medians pivot-কে নিশ্চিতভাবে মাঝামাঝি করে: উপাদানগুলোকে ৫টি করে দলে ভাগ, প্রতিটি দলের median, তারপর সেই ⌈n/5⌉টি median-এর median রিকার্শনে বের করে pivot. এমন pivot-এর চেয়ে ছোট অন্তত প্রায় ৩n/10টি আর বড় অন্তত প্রায় ৩n/10টি উপাদান থাকে, তাই পরের ধাপে বড়জোর প্রায় 7n/10. Recurrence:

এটি অসমান ভাগ, master theorem সরাসরি খাটে না। কিন্তু recursion tree-তে প্রতিটি তলার মোট আকার আগের তলার 1/5 + 7/10 = 9/10 গুণ — কমতে থাকা জ্যামিতিক ধারা, যোগফল n·(1 + 0.9 + 0.81 + …) = 10n. তাই Θ(n), worst case-এও। ৫-এর বদলে ৩ করে দল বানালে ভাগ হয় 1/3 + 2/3 = 1 — ধারা আর কমে না, ফল Θ(n log n). দলের আকার বাছাইয়ের পেছনে এই হিসাবটিই কাজ করে। বিস্তারিত প্রমাণ পর্ব ২২-এ আছে; বাস্তবে ধ্রুবক বড় বলে randomized quickselect-ই বেশি চলে।

যখন উপ-সমস্যাগুলো মিলে যায়

Fibonacci সংখ্যার সরাসরি রিকার্সিভ সংজ্ঞাটি দেখতে নিখুঁত divide and conquer: F(n) = F(n−1) + F(n−2). কিন্তু recursion tree আঁকুন — F(n−2) একবার সরাসরি ডাকা হয়, আবার F(n−1)-এর ভেতর থেকেও। F(n−3) তিনবার, F(n−4) পাঁচবার… মোট ডাক প্রায় φⁿ (φ ≈ 1.618) — exponential. এখানে overlapping subproblems: উপ-সমস্যাগুলো স্বাধীন নয়, একই প্রশ্ন বারবার আসে। সমাধান সহজ — একবার সমাধান করে মনে রাখুন (memoization), বা ছোট থেকে বড় টেবিল ভরুন। এটিই dynamic programming, পর্ব ৫২-এর বিষয়; তখন F(n) লাগে Θ(n) সময়।

আরও দুটি ফাঁদ। এক: ভাগ খুব অসম হলে — যেমন সবসময় ১ আর n−1 — গভীরতা n হয়ে যায়; quick sort-এর worst case ঠিক এটি (পর্ব ২২). দুই: combine যদি উপ-সমস্যার সংখ্যার চেয়ে দ্রুত বাড়ে, তাহলে ভাগ করে কিছু বাঁচে না — যেমন অংশ ৫-এর শুরুর 8T(n/2) + n², যা সাধারণ Θ(n³)-ই রয়ে গেল। কোনো divide and conquer লেখার আগে recurrence লিখুন আর master theorem দিয়ে দেখে নিন — লাভ আছে কি না, কাগজেই বোঝা যায়।

৮ · AI-এর যুগে divide and conquer

Divide and conquer কোড AI খুব সহজে লিখে দেয় — দেখতেও সুন্দর লাগে। কিন্তু এই কোডের ভুলগুলো সাধারণত ছোট ইনপুটে ধরা পড়ে না: অসীম রিকার্শন বড় ইনপুটে, ভুল complexity শুধু বড় n-এ, আর মাঝ-পার ক্ষেত্রের ভুল শুধু বিশেষ বিন্যাসে। যেসব ভুল বারবার দেখা যায়:

AI-এর দেওয়া divide and conquer কোড হাতে পেলে — যাচাইয়ের চেকলিস্ট

• কোড থেকে নিজে recurrence লিখেছি — ডাকের সংখ্যা, আকার, আর combine-এর আসল খরচ (লুকানো sort বা কপিসহ)
• আকার ০, ১, ২, ৩-এর ইনপুট হাতে trace করে দেখেছি রিকার্শন থামে আর ঠিক উত্তর দেয়
• Combine ধাপের সঠিকতা নিজের ভাষায় ব্যাখ্যা করতে পারি (মাঝ-পার ক্ষেত্র কেন সব সম্ভাবনা ধরে)
• Brute force-এর সাথে কয়েক হাজার এলোমেলো ছোট ইনপুটে stress test চালিয়েছি
• Edge case দেখেছি: সব ঋণাত্মক, সব সমান, একই x-এ সব বিন্দু, খুব বড় সংখ্যা (overflow)
• n দ্বিগুণ করে সময় মেপে দেখেছি বৃদ্ধির হার দাবির সাথে মেলে কি না