Chia để trị

/* Dat  c = log_b(a).  So sanh f(n) voi n^c:

   Truong hop 1:  f(n) = O(n^(c-e)) voi e > 0
                  Cac bai toan con chi phoi.
                  T(n) = Theta(n^c)

   Truong hop 2:  f(n) = Theta(n^c)
                  Hai ben ngang nhau.
                  T(n) = Theta(n^c log n)

   Truong hop 3:  f(n) = Omega(n^(c+e)) voi e > 0, va thoa dieu kien deu
                  Buoc chia va gop chi phoi.
                  T(n) = Theta(f(n))                                 */

#include <stdint.h>

/* Tính (a mũ n) chia dư mod, trong O(log n).

   Ý tưởng:  a^n = (a^(n/2))^2        nếu n chẵn
             a^n = a * (a^(n/2))^2    nếu n lẻ

   Bản lặp đọc n theo bit từ thấp lên cao. */
uint64_t luy_thua(uint64_t a, uint64_t n, uint64_t mod)
{
    uint64_t kq = 1 % mod;      /* xử lý đúng cả khi mod bằng 1 */

    a %= mod;

    while (n > 0) {
        if (n & 1) kq = kq * a % mod;      /* bit hiện tại là 1 */

        a = a * a % mod;                   /* bình phương cho bit kế */
        n >>= 1;
    }

    return kq;
}

/* Vet voi a = 3, n = 13, mod = 1000000007

   13 trong nhi phan la 1101, tuc 13 = 8 + 4 + 1
   nen 3^13 = 3^8 * 3^4 * 3^1

   buoc  n     bit  a hien tai    kq
     0   13     1   3             3
     1    6     0   9             3
     2    3     1   81            243
     3    1     1   6561          1594323
     4    0     -   -             ket qua 1594323

   Kiem: 3^13 = 1594323. Dung. Chi 4 vong lap thay vi 13. */

#include <stdint.h>

/* Quan sat:

     [ F(n+1) ]   [ 1 1 ]^n   [ F(1) ]
     [ F(n)   ] = [ 1 0 ]   * [ F(0) ]

   Nen tinh F(n) bang cach nang ma tran 2x2 len luy thua n,
   dung luy thua nhanh -> O(log n) phep nhan ma tran. */

typedef struct { uint64_t m[2][2]; } M2;

static M2 nhan(M2 x, M2 y, uint64_t mod)
{
    M2 z;

    for (int i = 0; i < 2; ++i)
        for (int j = 0; j < 2; ++j) {
            __uint128_t s = 0;

            for (int k = 0; k < 2; ++k)
                s += (__uint128_t)x.m[i][k] * y.m[k][j];

            z.m[i][j] = (uint64_t)(s % mod);
        }

    return z;
}

uint64_t fib_ma_tran(uint64_t n, uint64_t mod)
{
    M2 kq   = { { { 1, 0 }, { 0, 1 } } };      /* ma trận đơn vị */
    M2 nen  = { { { 1, 1 }, { 1, 0 } } };

    while (n > 0) {
        if (n & 1) kq = nhan(kq, nen, mod);

        nen = nhan(nen, nen, mod);
        n >>= 1;
    }

    return kq.m[0][1];      /* F(n) */
}

/* Nhan hai ma tran n x n theo dinh nghia: O(n^3).

   Strassen chia moi ma tran thanh bon khoi n/2 x n/2,
   roi dung BAY phep nhan thay vi TAM, bang cach cong tru kheo leo.

     T(n) = 7 T(n/2) + O(n^2)  ->  O(n^log2(7)) = O(n^2.807)

   Voi n = 1024:
     n^3     = 1073741824
     n^2.807 =  454875234    <- it hon 2.4 lan

   Nhung hang so an lon hon nhieu, va no khong on dinh ve so hoc,
   nen chi dang dung tu khoang n = 1000 tro len.                */

#include <math.h>
#include <stdlib.h>

typedef struct { double x, y; } Diem;

static double kc(Diem a, Diem b)
{
    double dx = a.x - b.x, dy = a.y - b.y;

    return sqrt(dx * dx + dy * dy);
}

/* px đã sắp theo x, py đã sắp theo y, cùng chứa n điểm. */
static double gan_nhat(Diem *px, Diem *py, size_t n, Diem *tam)
{
    if (n <= 3) {                          /* trường hợp cơ sở, vét cạn */
        double d = 1e308;

        for (size_t i = 0; i < n; ++i)
            for (size_t j = i + 1; j < n; ++j) {
                double t = kc(px[i], px[j]);

                if (t < d) d = t;
            }

        return d;
    }

    size_t giua  = n / 2;
    double x_mid = px[giua].x;

    /* Tách py thành hai nửa theo x_mid, giữ nguyên thứ tự theo y */
    Diem  *py_t = tam;
    Diem  *py_p = tam + giua;
    size_t nt = 0, np = 0;

    for (size_t i = 0; i < n; ++i) {
        if (nt < giua && py[i].x <= x_mid) py_t[nt++] = py[i];
        else                               py_p[np++] = py[i];
    }

    double dt = gan_nhat(px, py_t, giua, tam + n);
    double dp = gan_nhat(px + giua, py_p, n - giua, tam + n);
    double d  = (dt < dp) ? dt : dp;

    /* Dải giữa: chỉ những điểm cách đường chia dưới d mới đáng xét */
    Diem  *dai = tam + n;
    size_t nd  = 0;

    for (size_t i = 0; i < n; ++i)
        if (fabs(py[i].x - x_mid) < d)
            dai[nd++] = py[i];

    /* Điểm mấu chốt: với mỗi điểm, chỉ cần so với NHIỀU NHẤT 7 điểm
       kế tiếp trong dải, vì trong hình chữ nhật d nhân 2d không thể
       có quá 8 điểm cách nhau ít nhất d. */
    for (size_t i = 0; i < nd; ++i)
        for (size_t j = i + 1; j < nd && py[j].y - py[i].y < d; ++j) {
            double t = kc(dai[i], dai[j]);

            if (t < d) d = t;
        }

    return d;
}