Cache-aware ও Cache-oblivious — Big-O যখন মিথ্যা বলে

import numpy as np
import time


def bench(fn, *args, reps=5):
    best = min(
        (lambda t0: (fn(*args), time.perf_counter() - t0)[1])(time.perf_counter())
        for _ in range(reps))
    return best


# ---------- part ১: AoS vs SoA structured-sum ----------

N = 10_000_000
rng = np.random.default_rng(42)

# AoS: one interleaved (N, 7) array; sum columns 4,5,6 (velocities)
aos = rng.random((N, 7))
t_aos = bench(lambda: aos[:, 4].sum() + aos[:, 5].sum() + aos[:, 6].sum())

# SoA: three separate contiguous velocity arrays
va = rng.random(N); vb = rng.random(N); vc = rng.random(N)
t_soa = bench(lambda: va.sum() + vb.sum() + vc.sum())

print(f'AoS interleaved sum : {t_aos*1000:8.1f} ms')
print(f'SoA contiguous sum  : {t_soa*1000:8.1f} ms')
print(f'gap                 : {t_aos/t_soa:8.2f}x')


# ---------- part ২: matmul ijk vs ikj (pure python lists, small n) ----------

def matmul_ijk(A, B, n):
    C = [[0.0] * n for _ in range(n)]
    for i in range(n):
        for j in range(n):
            s = 0.0
            for k in range(n):
                s += A[i][k] * B[k][j]     # B column-walk: stride-n
            C[i][j] = s
    return C


def matmul_ikj(A, B, n):
    C = [[0.0] * n for _ in range(n)]
    for i in range(n):
        Ci = C[i]; Ai = A[i]
        for k in range(n):
            a = Ai[k]; Bk = B[k]
            for j in range(n):
                Ci[j] += a * Bk[j]          # all walks row-sequential
    return C


n = 220
Arng = [[float((i * j) % 7) for j in range(n)] for i in range(n)]
Brng = [[float((i + j) % 5) for j in range(n)] for i in range(n)]
R1 = matmul_ijk(Arng, Brng, n)
t_ijk = bench(matmul_ijk, Arng, Brng, n, reps=1)
R2 = matmul_ikj(Arng, Brng, n)
t_ikj = bench(matmul_ikj, Arng, Brng, n, reps=1)
assert R1 == R2
gflops = lambda t: 2 * n**3 / t / 1e9
print(f'\nmatmul n={n} (pure python):')
print(f'  ijk : {t_ijk*1000:9.0f} ms  ({gflops(t_ijk):5.3f} GFLOPS)')
print(f'  ikj : {t_ikj*1000:9.0f} ms  ({gflops(t_ikj):5.3f} GFLOPS)')
print(f'  speedup: {t_ijk/t_ikj:.2f}x')

# numpy reference — BLAS blocked+tuned reality check
Anp = rng.random((512, 512)); Bnp = rng.random((512, 512))
t_blas = bench(lambda: Anp @ Bnp, reps=3)
print(f'  numpy 512² BLAS: {t_blas*1000:.1f} ms '
      f'({2*512**3/t_blas/1e9:.1f} GFLOPS — tuned blocking+SIMD)')


# ---------- part ৩: array vs pointer-chain traversal proxy ----------

import sys

M = 400_000
vals_arr = list(range(M))                    # contiguous index space
# pointer-chain: shuffled linked-list simulation via index-hopping
perm = list(range(M)); rng.shuffle(perm)
next_idx = perm[:]                           # random successor mapping

def traverse_array():
    s = 0
    for i in range(M):
        s += vals_arr[i]                     # sequential
    return s

def traverse_chain():
    s, i = 0, 0
    for _ in range(M):
        s += vals_arr[i]
        i = next_idx[i]                      # random hop = cache miss proxy
    return s

ta = bench(traverse_array)
tc = bench(traverse_chain)
print(f'\ntraversal {M:,} nodes:')
print(f'  sequential array : {ta*1000:8.1f} ms')
print(f'  random chain     : {tc*1000:8.1f} ms   ratio {tc/ta:5.2f}x')

সাধারণ ফলাফল (M-series laptop):

AoS interleaved sum :   129.2 ms
SoA contiguous sum  :    11.7 ms
gap                 :    11.04x

matmul n=220 (pure python):
  ijk :      1384 ms  (0.015 GFLOPS)
  ikj :      1040 ms  (0.020 GFLOPS)
  speedup: 1.33x
  numpy 512² BLAS: 2.0 ms (135.0 GFLOPS — tuned blocking+SIMD)

traversal 400,000 nodes:
  sequential array :     30.7 ms
  random chain     :    162.1 ms   ratio  5.28x

পড়ার বিষয়:

১. AoS→SoA gap ~১১× — numpy strided view (aos[:,4]) প্রতিটা row-এ ৬টা অবাঞ্ছিত field skip করে; contiguous SoA pure streaming। Native C-তে gap সাধারণত ২-৪× — layout-এর bandwidth tax universal।

২. Pure-Python matmul-এ ikj মাত্র 1.33× — interpreter overhead locality gain-এর বেশিরভাগ mask করে। আসল গল্পটা BLAS-এর সংখ্যায়: numpy 135 GFLOPS বনাম pure-Python 0.02 — ~৭০০০× gap টাই compiled loop + tiling + SIMD-এর combined verdict; C/C++-এ শুধু reorder-ই সাধারণত ৪-৮× দেয়।

৩. Random-chain 5.3× ধীর — pointer-chasing প্রতি hop-এ cache-miss wait; Python object heap scattered হওয়ায় এটা DRAM-latency proxy হিসেবে নিখুঁত।