Add L3

Author: hxuhack

README.md

@@ -5,9 +5,9 @@
 
 | 周 | 内容 | 资料 | 代码 |
 |---|---|---|---|
-| 1 | 使用Python编写神经网络 | [课件](notes/L1-使用python编写神经网络.pdf) | [factorial-relu](code/l1-factorial-relu.py), [factorial-sigmoid](code/l1-factorial-sigmoid.py) |
+| 1 | 使用Python编写神经网络 | [课件](notes/L1-使用Python编写神经网络.pdf) | [factorial-relu](code/l1-factorial-relu.py), [factorial-sigmoid](code/l1-factorial-sigmoid.py) |
 | 2 | Python程序的解析和运行 | [课件](notes/L2-Python程序的解析和运行.pdf) | [regex](code/l2-regex.py), [pytool](code/l2-pytool.py), [relu](code/l2-relu.py) |
-| 3 | Python与Native Code交互 | | | 
+| 3 | Python与Native Code交互 | [课件](notes/L3-Python与native交互.pdf) | [bench_dot(py)](code/l3-1-bench_dot.py), [bench_dot(c)](code/l3-1-bench_dot.c), [bench_numba](code/l3-2-bench_numba.py), [decorator(py)](code/l3-3-deco.py), [decorator(c)](code/l3-3-vec.c) | 
 | 4 | Python多线程 | | | 
 | 5 | 课堂反转 | |
 | 6 | 使用PyTorch编写和训练神经网络模型 | | | 

code/l3-1-bench_dot.c

@@ -0,0 +1,9 @@
+#include <stddef.h>
+
+double dot_c(const double *a, const double *b, size_t n) {
+    double s = 0.0;
+    for (size_t i = 0; i < n; i++) {
+        s += a[i] * b[i];
+    }
+    return s;
+}

code/l3-1-bench_dot.py

@@ -0,0 +1,77 @@
+import ctypes
+import random
+import time
+import numpy as np
+from numba import njit
+
+# -------------------------------
+# Configuration
+# -------------------------------
+N = 5_000_000      # vector length
+SEED = 42
+
+random.seed(SEED)
+
+# -------------------------------
+# Data preparation
+# -------------------------------
+a_list = [random.random() for _ in range(N)]
+b_list = [random.random() for _ in range(N)]
+
+a_np = np.array(a_list, dtype=np.float64)
+b_np = np.array(b_list, dtype=np.float64)
+
+
+# -------------------------------
+# Implementations
+# -------------------------------
+def dot_python(a, b):
+    s = 0.0
+    for i in range(len(a)):
+        s += a[i] * b[i]
+    return s
+
+
+def dot_numpy(a, b):
+    return np.dot(a, b)
+
+# -------------------------------
+# Benchmark helper
+# -------------------------------
+def benchmark(func, *args):
+    times = []
+    start = time.perf_counter()
+    func(*args)
+    end = time.perf_counter()
+    return end - start
+
+
+lib = ctypes.cdll.LoadLibrary("./libdot.so")
+
+lib.dot_c.restype = ctypes.c_double
+lib.dot_c.argtypes = [
+    ctypes.POINTER(ctypes.c_double),
+    ctypes.POINTER(ctypes.c_double),
+    ctypes.c_int,
+]
+
+def dot_c(a, b):
+    return lib.dot_c(
+        a.ctypes.data_as(ctypes.POINTER(ctypes.c_double)),
+        b.ctypes.data_as(ctypes.POINTER(ctypes.c_double)),
+        a.size
+    )
+
+# -------------------------------
+# Run benchmarks
+# -------------------------------
+t_py = benchmark(dot_python, a_list, b_list)
+t_np = benchmark(dot_numpy, a_np, b_np)
+t_c = benchmark(dot_c, a_np, b_np)
+
+# -------------------------------
+# Results
+# -------------------------------
+print(f"Python loop : {t_py:.4f} s")
+print(f"NumPy dot   : {t_np:.4f} s")
+print(f"C dot   : {t_np:.4f} s")

code/l3-2-bench_numba.py

@@ -0,0 +1,45 @@
+import random
+import time
+import numpy as np
+from numba import njit
+
+N = 5_000_000
+
+def dot_python(a, b):
+    s = 0.0
+    for i in range(len(a)):
+        s += a[i] * b[i]
+    return s
+
+@njit
+def dot_numba(a, b):
+    s = 0.0
+    for i in range(len(a)):
+        s += a[i] * b[i]
+    return s
+
+# 生成数据
+a_list = [random.random() for _ in range(N)]
+b_list = [random.random() for _ in range(N)]
+
+a_np = np.array(a_list)
+b_np = np.array(b_list)
+
+# Python benchmark
+t0 = time.perf_counter()
+dot_python(a_np, b_np)
+t1 = time.perf_counter()
+
+# Numba 第一次调用（包含 JIT 编译）
+t2 = time.perf_counter()
+dot_numba(a_np, b_np)
+t3 = time.perf_counter()
+
+# Numba 第二次调用（已编译）
+t4 = time.perf_counter()
+dot_numba(a_np, b_np)
+t5 = time.perf_counter()
+
+print("Python time:", t1 - t0)
+print("Numba first call (compile + run):", t3 - t2)
+print("Numba second call (run only):", t5 - t4)

code/l3-3-deco.py

@@ -0,0 +1,117 @@
+import ast
+import inspect
+import ctypes
+import numpy as np
+
+# ============================================================
+# Native primitive kernels (precompiled shared library)
+# ============================================================
+
+lib = ctypes.CDLL("./libvec.so")
+
+lib.vec_elem_mul.argtypes = lib.vec_elem_add.argtypes = [
+    ctypes.POINTER(ctypes.c_double),
+    ctypes.POINTER(ctypes.c_double),
+    ctypes.POINTER(ctypes.c_double),
+    ctypes.c_int,
+]
+lib.vec_elem_add.restype = lib.vec_elem_mul.restype = None
+
+
+# ============================================================
+# @kernel decorator: Python AST -> lowering -> native execution
+# ============================================================
+
+def kernel(func):
+    """
+    Treat the function body as an embedded DSL.
+    The body is parsed, not executed.
+    """
+    # -------- Parse Python AST --------
+    src = inspect.getsource(func)
+    tree = ast.parse(src)
+    func_def = tree.body[0]
+
+    # Expect:
+    #   return <expr>
+    return_stmt = func_def.body[0]
+
+    expr = return_stmt.value
+    arg_names = [arg.arg for arg in func_def.args.args]
+
+    # -------- Lowering: AST -> primitive ops --------
+    def lower(node, env, n):
+        if isinstance(node, ast.Name):
+            return env[node.id]
+
+        if isinstance(node, ast.BinOp):
+            left = lower(node.left, env, n)
+            right = lower(node.right, env, n)
+            out = np.empty_like(left)
+
+            if isinstance(node.op, ast.Mult):
+                lib.vec_elem_mul(
+                    left.ctypes.data_as(ctypes.POINTER(ctypes.c_double)),
+                    right.ctypes.data_as(ctypes.POINTER(ctypes.c_double)),
+                    out.ctypes.data_as(ctypes.POINTER(ctypes.c_double)),
+                    n,
+                )
+            elif isinstance(node.op, ast.Add):
+                lib.vec_elem_add(
+                    left.ctypes.data_as(ctypes.POINTER(ctypes.c_double)),
+                    right.ctypes.data_as(ctypes.POINTER(ctypes.c_double)),
+                    out.ctypes.data_as(ctypes.POINTER(ctypes.c_double)),
+                    n,
+                )
+            else:
+                raise NotImplementedError("Unsupported operator")
+
+            return out
+
+        raise NotImplementedError("Unsupported AST node")
+
+    # -------- Runtime wrapper --------
+    def wrapper(*args):
+        arrays = [np.asarray(a, dtype=np.float64) for a in args]
+        n = arrays[0].size
+
+        env = dict(zip(arg_names, arrays))
+        result = lower(expr, env, n)
+        return result
+
+    return wrapper
+
+
+# ============================================================
+# User-defined kernels (pure Python, no strings)
+# ============================================================
+
+@kernel
+def vec_elem_mul(a, b):
+
+    # element-wise multiplication
+    return a * b
+
+@kernel
+def vec_elem_add(a, b):
+    # element-wise addition
+    return a + b
+
+@kernel
+def vec_elem_fma(a, b, c):
+    # element-wise fused multiply-add: a * b + c
+    return (a * b) + c
+
+# ============================================================
+# Test
+# ============================================================
+
+if __name__ == "__main__":
+    a = [1.0, 2.0, 3.0]
+    b = [4.0, 5.0, 6.0]
+    c = [10.0, 10.0, 10.0]
+
+    print("mul:", vec_elem_mul(a, b))
+    print("add:", vec_elem_add(a, b))
+    print("fma:", vec_elem_fma(a, b, c))
+

code/l3-3-vec.c

@@ -0,0 +1,13 @@
+#include <stdio.h>
+
+void vec_elem_mul(double* a, double* b, double* c, int n) {
+    for (int i = 0; i < n; i++) {
+        c[i] = a[i] * b[i];
+    }
+}
+
+void vec_elem_add(double* a, double* b, double* c, int n) {
+    for (int i=0; i<n; i++) { 
+        c[i] = a[i] + b[i];
+    }
+}