简介这份C工程面向天体物理模拟与数值算法学习者围绕Barnes-Hut算法高效求解N体问题。它用八叉树对三维空间分层划分并按远近阈值以节点质心和总质量近似远处粒子团从而把传统逐对引力计算的复杂度显著降低从粒子模型、八叉树构建到引力计算与时间推进完整呈现了近似算法的实现链路代码中还能看到欧拉、龙格-库塔、Adams-Bashforth等积分器与SDL可视化窗口的组织方式适合研究多体模拟、空间数据结构与性能优化。压缩包共60个文件约54KB以21个cpp源文件和22个h头文件为主体另有4个mk构建脚本、工程配置与说明文档覆盖粒子模型、八叉树、引力计算、积分推进和窗口渲染等模块。已有329人学习可作为课程设计、算法复现和C工程实践的可编译底稿便于对照理解近似精度与计算效率的权衡。1. 一万个天体、每一步五千万对引力Barnes-Hut 省掉的是哪部分一万个天体放进二维方框做引力演化每一步要算多少对引力N(N-1)/2 ≈ 5×10⁷。单步还能忍一套轨道跑十万步就是 5×10¹² 次开方和除法一台工作站从天亮跑到天黑。Barnes-Hut 算法不改变物理它把「远处的天体团」折叠成一个质点让每个粒子的受力计算从一万次配对降到几百次整体复杂度从 O(N²) 压到 O(N log N)。这套 Barnes-Hut Simulator 的内核就三件事四叉树三维是八叉树、开角判据、蛙跳积分器。它拿的是「可接受的近似误差换十倍以上速度」这笔账。适合写星团演化、暗物质晕塌缩、分子动力学粗粒化、粒子群可视化的人不适合需要长时间精确轨道的双星系统也不适合 N1000 的场景——那种规模直接求和反而更快还少一层玄学。下面按「树怎么建 → 力怎么算 → 时间怎么走 → 哪里会崩」的顺序拆开代码用 Python NumPy 写重点在结构而不在语言换成 C 或 Rust 时一一对应即可。2. 从粒子表到四叉树Barnes-Hut 的节点里到底存什么2.1 一个节点必须携带的五样东西四叉树节点的最小字段集是子树总质量 mass、子树质心 (mx, my)、节点半边长 half、四个子节点指针 kids、以及叶节点持有的粒子下标数组 body。五个字段各自对应一件事少一个都要在后面付代价。mass 和 mx/my 是远场近似的全部依据。当判据认定某个节点「足够远」我们用它整体的质心和一个等效质量代替它内部所有粒子这一步等价于多极展开只保留零阶单极项误差量级由下一节的开角判据控制。half 提供空间尺度判据里的 s 就是 2×half也就是节点边长。kids 决定还能不能继续细分body 决定叶子上要老老实实做几次精确配对。这里最容易踩的坑是把「节点几何中心 (cx, cy)」当成质心用。几何中心是建树时切分出来的方框中心质心是质量加权平均两者在节点内质量分布不对称时可以差出半个边长。用错的那个力会有系统性偏差而且不会报错只会让轨道慢慢漂——这种血泪经验后文还会专门讲。import numpy as np class QuadNode: __slots__ (cx, cy, half, mass, mx, my, kids, body) def __init__(self, cx, cy, half): self.cx, self.cy, self.half cx, cy, half self.mass 0.0 # 子树总质量 self.mx, self.my 0.0, 0.0 # 子树质心质量加权不是几何中心 self.kids None # None 表示叶节点 self.body None # 叶节点持有的粒子下标数组__slots__在这里不是炫技N10⁵ 时节点数大约是粒子数的两到三倍每个节点省掉一个__dict__内存占用和缓存命中率都会明显不同。这一点在 Python 里只是小优化在 C 里换成紧凑数组就是数量级差别。2.2 自顶向下按象限分治别逐个插入建树有两种写法。一种是逐个粒子插入每遇到已占用的叶子就分裂另一种是先算全体包围盒再按象限把粒子下标数组递归切分。前者写起来顺手后者快得多而且在 Python 这种解释执行的环境里差距会被放大。原因在于下标数组的切分是一次向量化的布尔索引一次处理一整批粒子逐个插入则是 N 次指针跳跃每次都要走一遍从根到叶的路径。用分治写法建树整体是 O(N log N)递归深度约 log₄N。def build_tree(pos, mass, leaf_cap1, max_depth24): lo, hi pos.min(axis0), pos.max(axis0) center 0.5 * (lo hi) half 0.5 * float(np.max(hi - lo)) 1e-9 # 取最大跨度保证是正方形 root QuadNode(center[0], center[1], half) _build(root, np.arange(len(pos)), pos, mass, leaf_cap, max_depth) return root def _build(node, idx, pos, mass, leaf_cap, depth): w mass[idx] node.mass float(w.sum()) node.mx float((w * pos[idx, 0]).sum()) / node.mass node.my float((w * pos[idx, 1]).sum()) / node.mass # 两个终止条件粒子数不超过叶容量或已到深度上限重合点保护 if len(idx) leaf_cap or depth 0: node.body idx return rel pos[idx] - np.array([node.cx, node.cy]) # 象限编号x 正负 2 * y 正负 → 0..3对应左下/右下/左上/右上 quad (rel[:, 0] 0).astype(np.int8) 2 * (rel[:, 1] 0).astype(np.int8) h node.half * 0.5 offsets ((-h, -h), (h, -h), (-h, h), (h, h)) node.kids [] for q in range(4): sub idx[quad q] ox, oy offsets[q] child QuadNode(node.cx ox, node.cy oy, h) if len(sub): _build(child, sub, pos, mass, leaf_cap, depth - 1) node.kids.append(child)代码逻辑顺着读先把这个节点负责的那批粒子的总质量和质心算出来这是后面前序遍历力的前提再判断要不要继续分。要分就把粒子按所在象限打标签四个子节点各自拿走一批递归下去。参数上有三个可以调。leaf_cap是叶容量取 1 时树最深、精度最高、节点数最多取 4~8 可以显著减少节点总数和递归开销代价是叶子内部仍然要做几次精确配对——因为叶子上的粒子没法再近似了。实际项目里我一般先用leaf_cap1跑通确认精度对得上再往上调。max_depth24是重合点的后悔药如果两个粒子坐标完全一致象限切分永远分不开它们没有深度上限就会一路递归到栈溢出。取 24 意味着最小格子边长是初始边长的 2⁻²⁴早就低于双精度在这个尺度上的分辨能力了到此为止是合理的。还有一处容易忽略的细节根节点的half用的是max(hi - lo) / 2而不是分别取 x、y 跨度。因为四叉树的子节点是正方形根节点不是正方形的话分割出来的子节点边长会不一致开角判据里的 s 就失去了统一含义。2.3 质心必须由子节点回算不能靠插入时累加上面这份实现里质心是在每次_build入口用mass[idx]直接对这批粒子求加权平均的所以严格来说不存在「回填」这一步。如果你采用逐个插入的写法就必须在每次分裂时把父节点原有的质量重新分配到子节点并在整棵树建完后做一次后序遍历自底向上重算。两种写法的结果应当一致但逐个插入的写法在分裂时特别容易漏掉「把已有那个粒子的贡献搬到新子节点」这一步症状是父节点的 mass 对、质心偏越靠近树顶偏得越离谱而粒子少的地方又看不出来。所以我更倾向分治写法质心由定义直接算出来不依赖任何中间状态的正确性。代价是每个节点都要对它那一批粒子做一次加权求和整体是 O(N log N) 而不是 O(N)。这个代价在 N10⁵ 时大约是建树总时间的 30%~40%可以接受真要抠就在递归时把子节点已经算好的 mass 和 mx/my 传上来合并回到 O(N)。2.4 三条建完就该跑的自检断言树建错的典型表现不是崩溃是「跑得动但结果不对」。所以在进入力的计算之前加三条断言把它摁住def check_tree(root, pos, mass, tol1e-10): total_m root.mass assert abs(total_m - mass.sum()) tol * mass.sum(), 总质量不守恒 cx float((mass * pos[:, 0]).sum()) / mass.sum() cy float((mass * pos[:, 1]).sum()) / mass.sum() assert abs(root.mx - cx) tol and abs(root.my - cy) tol, 根节点质心不对 # 总质量为零的空节点不该被当成有效源 assert total_m 0.0, 全质量为零第一条查的是质量有没有在分裂过程中丢失或重复计入。第二条查的是「质心 vs 几何中心」这个最容易写错的地方——如果mx恰好等于root.cx八成是把几何中心赋过去了。第三条是为后面力计算的分母兜底mass 为零的节点在近似公式里会给出 0 加速度看起来无害但如果它同时又是某个粒子的唯一父节点说明这颗粒子根本没进树。3. 开角判据与受力遍历一次递归同时买到精度和速度3.1 s/d θ 到底近似掉了什么对一个距离为 d、边长为 s 的节点Barnes-Hut 用它的质心代替内部所有质量。这个近似展开到单极项就截止丢掉的是偶极项质心定义保证它为零和四极项及更高阶。四极项的相对大小大约是 (s/d)² 量级所以判据写成 s/d θ 时力的相对误差大致按 θ² 走。这就解释了两件事。第一θ 从 0.5 放宽到 1.0误差不是翻倍而是涨到四倍左右实测里经常看到「θ 调到 0.8 就明显不对了」。第二θ 趋近 0 时这个近似退化成精确配对算法复杂度回到 O(N²)——所以 θ 不是「越小越好」的精度旋钮它是精度和速度之间唯一的那根杠杆。3.2 单粒子受力遍历的完整实现遍历从根开始。对当前节点如果它是叶子就对叶内每个粒子做精确配对如果它是内部节点且满足开角判据就用它的 mass 和质心算一次近似力然后返回否则递归进四个子节点。整个函数是尾递归形式深度受max_depth约束。def accelerations(root, pos, mass, G1.0, theta0.5, eps0.01): acc np.zeros_like(pos) for i in range(len(pos)): _walk(root, i, pos, mass, G, theta, eps, acc) return acc def _walk(node, i, pos, mass, G, theta, eps, acc): if node.mass 0.0: return dx node.mx - pos[i, 0] dy node.my - pos[i, 1] # 叶子节点逐个粒子精确计算不做近似 if node.kids is None: for j in node.body: if j i: continue ex pos[j, 0] - pos[i, 0] ey pos[j, 1] - pos[i, 1] s2 ex * ex ey * ey eps * eps # 软化因子在这里进分母 inv G * mass[j] / (s2 * np.sqrt(s2)) acc[i, 0] ex * inv acc[i, 1] ey * inv return d2 dx * dx dy * dy d np.sqrt(d2) # 内部节点s 2 * half 是节点边长s/d theta 时整体近似 if d 0.0 and (2.0 * node.half) / d theta: s2 d2 eps * eps inv G * node.mass / (s2 * np.sqrt(s2)) acc[i, 0] dx * inv acc[i, 1] dy * inv return for child in node.kids: _walk(child, i, pos, mass, G, theta, eps, acc)逐段看逻辑。开头的node.mass 0.0是空节点早退用而不是是为了顺手挡掉负质量这种不可能但一旦出现就会静默污染全局的输入。叶子分支里j i的跳过不能省叶节点持有的是下标数组自身一定在里面不跳过就是自相互作用会得到一个巨大的虚假加速度。近似分支里的d 0.0是防除零d恰好为零意味着粒子和节点质心重合这在实际演化里出现概率极低但一旦出现就是 NaN 扩散。近似公式acc G * M * Δr / (d² ε²)^{3/2}是标准的 Plummer 软化形式。注意 ε 是加在 d² 上而不是加在 d 上也不是加在最终的力上——加的位置不同近距离行为完全不同这一点在第 5 章会展开。3.3 θ 的取值一张可以直接抄的对照表θ 取值力的相对误差量级单粒子访问节点数典型用途0.20.5%高接近精确配对需要长时间守恒的轨道、小规模标定0.51%~3%中星团演化、暗物质晕默认起点0.73%~6%低大尺度统计量、形态可视化1.06%~15%最低只看粗略结构的快速预览表里的误差是量级估计不是可以直接引用的数字——它和粒子分布、软化因子、节点形状都有关系。稳妥做法是拿自己的初始条件用直接求和下一节的脚本算一遍精确加速度再扫一遍 θ画一条「θ — 相对误差」曲线选误差落在你可接受范围内、且访问节点数已经明显下降的那个拐点。我一般从 0.5 起步因为这个值在多数分布下误差已经压到个位数百分比而节点访问量相比精确配对少了两个数量级。还有一点θ 是全局常数还是逐节点变化对结果影响不小。有些实现会根据节点里粒子数动态放宽 θ粒子多、距离远就更敢近似这在密度反差极大的分布里有效但它会让误差不再有统一量纲做精度标定时很难说清是哪个 θ 在起作用。前期建议先固定全局 θ把其他环节调稳再考虑动态策略。3.4 软化因子 ε 不是精度参数是物理参数ε 的作用是把 1/r² 的奇点抹平。在真实的星团里两颗星靠得足够近时会有非引力相互作用或者直接碰撞软化就是在粗粒化掉这个过程。所以 ε 的选择应该跟你的物理问题挂钩模拟的是稀疏星团ε 取平均间距的 1/10 到 1/100 就够模拟的是致密系统ε 太小会让近距离交会主导整个演化时间步被迫压得很小。经验做法是先估一下初始分布的平均最近邻距离把 ε 设在它的 1/10 附近然后跑短轨道看能量曲线的斜率。如果能量快速漂移先怀疑时间步再怀疑 ε如果双星轨道被异常散射那基本是 ε 太小。4. 时间积分把加速度接成能跑十万步的轨道4.1 为什么显式 Euler 在 N 体里必然翻车力的计算再准积分器选错也白搭。显式 Euler 每一步用当前加速度同时更新速度和位置局部误差 O(dt²)而且是系统性地往系统里注入能量——在引力这种长程、近距离交会频繁的势场里几万步之后能量偏移足以让整个系统「热膨胀」到面目全非。这不是精度不够的问题是方法本身没有守恒结构。蛙跳积分leapfrogKDK 形式属于辛积分器它精确保持相空间体积能量误差不在时间上单调增长而是在一个有界区间里振荡。这个性质对 N 体模拟是决定性的你不需要很小的时间步就能让总能量在十万步后仍然待在同一条带子里。代价是它只有二阶精度补偿手段是「半步 kick、整步 drift、再半步 kick」。4.2 KDK 三步和树的重建时机def leapfrog_step(pos, vel, mass, dt, acc, cfg): vel 0.5 * dt * acc # 半步 kick pos dt * vel # 整步 drift root build_tree(pos, mass, cfg[leaf_cap], cfg[max_depth]) acc accelerations(root, pos, mass, Gcfg[G], thetacfg[theta], epscfg[eps]) vel 0.5 * dt * acc # 再半步 kick return pos, vel, acc顺序不能换。第一步用旧位置的加速度把速度推半步第二步用这个中间速度把位置推一整步第三步在新位置重算加速度再补半步速度。这样速度和位置在时间上错开半步正好构成一个辛映射。注意build_tree和accelerations都必须在pos更新之后调用——树是位置的函数位置变了树就作废了。这一步是整段代码里最容易被「优化」错的地方有人为了省时间在多个步长内复用同一棵树结果粒子已经移出原来的节点范围近似判据算出来的 s/d 完全是错的误差不会报错只会让轨道莫名其妙地发散。acc在步与步之间传递省掉一次重复计算。第一次调用前需要单独算一次acc从那之后每一次的返回值都是下一步要用的旧加速度。4.3 能量漂移是唯一可信的黑匣子读数N 体模拟出问题时画面往往还「看着像那么回事」只有能量曲线会说实话。总能量 E KE PEKE 是 Σ½mv²PE 是所有粒子对的 -Gmᵢmⱼ/rᵢⱼ带同样的软化。def energy(pos, vel, mass, G1.0, eps0.01): ke 0.5 * np.sum(mass * np.sum(vel * vel, axis1)) # 势能用直接求和O(N²)只在 N 较小或做基准时用 dx pos[:, None, 0] - pos[None, :, 0] dy pos[:, None, 1] - pos[None, :, 1] r2 dx * dx dy * dy eps * eps np.fill_diagonal(r2, np.inf) # 排除自相互作用 pe -0.5 * G * np.sum(mass[:, None] * mass[None, :] / np.sqrt(r2)) return ke pe逻辑上有两个点要抠。np.fill_diagonal(r2, np.inf)把自配对的距离置成无穷大对应的贡献是 0同时避免了对角线上的除零。前面的 0.5 是因为mass[:,None] * mass[None,:]这张 N×N 矩阵把每一对算了两遍。软化因子必须和力的计算里用的完全一致。如果力里用(d²ε²)^{3/2}、势能里用1/d能量曲线会呈现出「看起来在漂移」的假象实际上是两套势场在打架这种不一致比数值误差更难查因为它不随步长减小而收敛。判断标准很简单把 dt 减半再跑一遍如果能量漂移的形状和幅度都明显改变说明是时间步问题如果几乎不变说明误差来自树的近似θ或者积分器的守恒性质跟 dt 无关。这一步判断能省掉大量盲目调参。4.4 自适应时间步让近距离交会自己决定节奏固定 dt 的问题在于它必须迁就最极端的那个事件。一套星团里大部分粒子运动平缓个别粒子偶尔会有近距离交会用统一的小步长跑全程是巨大的浪费。常见做法是 Aarseth 判据的变体每个粒子按自己的加速度算一个特征时间 τᵢ η·sqrt(ε/|aᵢ|)取全局最小值作为这一步的 dt。η 一般取 0.01~0.05它控制的是「一步之内速度变化占特征速度的比例」。这么做能自动在系统剧烈变化时收缩步长平静时放大。代价是 dt 变成时变的严格的辛性质会打折扣——不过只要 dt 变化缓慢能量误差仍然能保持在有界振荡而不是单调漂移这一点在实践中比理论上的严格性更重要。如果观察到能量曲线变成阶梯状上跳多半是 dt 变化太剧烈把 η 调小或者对 dt 做平滑比如新 dt 不超过旧 dt 的 1.2 倍就能压住。5. 排查手册θ、ε 和退化树导致的五类典型故障5.1 力算错的三类症状以及该先看哪里Barnes-Hut 的 bug 有一个共同特点不崩溃。程序照跑图照出只有和基准对不上。所以排查顺序要固定下来别东一榔头西一棒子。第一顺位是能量漂移的形状。单调上升或下降查积分器和 dt有界振荡但幅值偏大查 θ 和 ε出现台阶状突变查自适应步长。第二顺位是加速度的逐粒子对比。拿 N200 的小规模用直接求和算出精确加速度和树方法逐分量相减看误差是均匀分布还是集中在某些粒子上。均匀分布是 θ 的问题集中在少数粒子通常是靠近节点边界的是判据实现的问题集中在质量最大的那几个粒子上是质心算错了。第三顺位才是看轨迹图。因为轨迹是最后才被误差污染的东西从它反推原因效率很低。5.2 五条翻车记录一、跑几十步后位置全是 NaN现象pos里出现nan能量变成-inf日志里没有任何异常。原因初始条件里有两个粒子坐标完全相同或者近到双精度在这个尺度上分辨不出。象限判据用的是 0两个点落在同一侧永远分不开递归一路到底直到栈溢出即使有深度上限兜住到了力的计算阶段 d≈0 又会除零。解决建树时加max_depth兜底本文实现里已有力计算里所有分母都过一遍软化因子 ε初始化生成粒子时对重复坐标加 1e-10 量级的随机扰动。三条里最不该省的是第一条因为它是唯一能防住栈溢出的。二、同一个 θ别人的代码误差 1%我的误差 8%现象对照直接求和相对误差差了好几倍但代码逻辑逐行比对看不出区别。原因θ 的比较方式没有统一约定。常见写法有三种s/d θ、s²/d² θ²等价、以及用节点对角线代替边长差一个 √2。更隐蔽的是 s 取边长还是半边长——差一个 2 的因子等效于 θ 差了 2 倍误差按 θ² 走就是 4 倍差距。解决在代码里把约定写死并注释清楚本文用的是「s 2 × half 即边长判据 s/d θ」然后用同一份初始条件、同一套直接求和基准扫一遍 θ 画出自己的误差曲线。别拿别人的 θ 数字直接用。三、总质量守恒但根节点的质心慢慢偏离现象root.mass一直对root.mx和全体粒子的加权平均对不上偏差随时间增大。原因用了节点的几何中心 (cx, cy) 代替质心做远场近似。几何中心在节点内质量分布对称时恰好等于质心所以小规模测试很容易蒙混过关一旦分布不对称偏差立刻显现而且它表现为力的系统性偏移不会报错。解决每个节点的 mx/my 必须由(mass * pos).sum() / mass.sum()算出来并且加一条断言根节点的 mx 与全体粒子的加权平均之差应当在 1e-10 量级。如果这两个数恰好相等而root.cx也相等反而要警惕——那说明初始分布太对称没测出问题。四、双星靠得很近时被弹飞能量突然跳一截现象两条轨道接近到最小距离时其中一个粒子的速度异常增大之后整团天体的能量曲线抬升一个台阶。原因软化因子 ε 只加在了力的大小上势能里用的还是1/r或者力和势能虽然都加了 ε但形式不同一个用(d²ε²)^{3/2}一个用d²ε²。两套势场不一致在近距离时表现出的等效能量不守恒。解决力和势能统一使用 Plummer 形式——力GmΔr/(d²ε²)^{3/2}势能-Gm/sqrt(d²ε²)。改完之后能量曲线的台阶应当消失剩下的平滑漂移才是数值误差。五、总能量单调上升系统整体「变热」现象几万步之后原本的团块逐渐散开速度弥散变大能量曲线一路向上。原因多半是 dt 太大遇上近距离交会。也有可能是 θ 太松导致树近似在每步引入的误差方向不是随机的而是有偏的近似总是低估远场的力长期累积就成了净能量注入。解决先把 θ 从 0.5 收到 0.3 重跑如果漂移速度明显下降问题在树近似如果几乎没变把 dt 减半再跑漂移幅度应随之减半二阶方法。两个都试过还漂检查自适应步长是否变化太剧烈。5.3 一个三十行的交叉验证脚本上面所有排查都指向同一件事你得有一个「随时能拿来对答案」的精确解。直接求和写起来很短值得常备。def direct_accelerations(pos, mass, G1.0, eps0.01): dx pos[:, None, 0] - pos[None, :, 0] # dx[i, j] x_i - x_j dy pos[:, None, 1] - pos[None, :, 1] r2 dx * dx dy * dy eps * eps np.fill_diagonal(r2, np.inf) # 排除自配对 inv G / (r2 * np.sqrt(r2)) # 1/(r²ε²)^{3/2} ax -np.sum(mass[None, :] * dx * inv, axis1) ay -np.sum(mass[None, :] * dy * inv, axis1) return np.stack([ax, ay], axis1) def compare(pos, mass, root, theta, eps): a_tree accelerations(root, pos, mass, thetatheta, epseps) a_ref direct_accelerations(pos, mass, epseps) rel np.linalg.norm(a_tree - a_ref) / np.linalg.norm(a_ref) worst np.max(np.linalg.norm(a_tree - a_ref, axis1) / np.linalg.norm(a_ref, axis1)) return rel, worst符号要留意dx[i, j] x_i - x_j而作用在 i 上的加速度指向 j所以整体带一个负号。这个负号和前面树遍历里的写法是等价的那边写的是pos[j] - pos[i]不带负号两种约定都对但混用会让误差看起来像 200% 而不是 2%。worst这个指标比全局相对误差更有用。全局误差被大量平缓的远场粒子摊薄了看起来总是很漂亮worst 会告诉你最差的那个粒子错了多少而轨道演化恰恰是被最差的那几个粒子带偏的。我一般要求 θ0.5 时 worst 控制在 5% 以内。6. 进阶给质心项补上四极矩以及我踩过的那个坑单极近似丢掉的是四极项。如果你想在同样的 θ 下把误差压下来最划算的改法不是调小 θ那是指数级的代价而是给每个节点额外存两个二阶矩分量力公式里补一项项存什么计算代价收益单极mass, mx, my已有基准四极Qxx, Qyy, Qxy建树时各多两次加权求和θ 相同下误差降一个量级八极三阶矩明显变重边际收益低不推荐四极修正的形状是把G·M·Δr/(d²ε²)^{3/2}换成G·M·Δr/(d²ε²)^{3/2} 修正项(Δr, Q)。修正项里包含Δr的三次方组合和 Q 的分量具体形式按二维展开写开有六项实践中多数人直接调用现成的多极展开库或者查表实现没必要手推。我更想讲的是另一样东西在 Python 里用对象指针建树N 超过 5×10⁴ 之后耗时的大头不是浮动运算而是属性访问和缓存未命中。那次翻车是把一份能跑 N10⁴ 的代码直接拿去跑 10⁵以为复杂度是 O(N log N) 所以只慢十倍结果慢了将近四十倍而且 profiler 显示 70% 的时间花在_walk的属性读取上。改法是把树压成几个并行的 NumPy 数组mass[]、mx[]、my[]、half[]、first_child[]、child_count[]本文还有配套的精品资源点击获取
简介这份C工程面向天体物理模拟与数值算法学习者围绕Barnes-Hut算法高效求解N体问题。它用八叉树对三维空间分层划分并按远近阈值以节点质心和总质量近似远处粒子团从而把传统逐对引力计算的复杂度显著降低从粒子模型、八叉树构建到引力计算与时间推进完整呈现了近似算法的实现链路代码中还能看到欧拉、龙格-库塔、Adams-Bashforth等积分器与SDL可视化窗口的组织方式适合研究多体模拟、空间数据结构与性能优化。压缩包共60个文件约54KB以21个cpp源文件和22个h头文件为主体另有4个mk构建脚本、工程配置与说明文档覆盖粒子模型、八叉树、引力计算、积分推进和窗口渲染等模块。已有329人学习可作为课程设计、算法复现和C工程实践的可编译底稿便于对照理解近似精度与计算效率的权衡。1. 一万个天体、每一步五千万对引力Barnes-Hut 省掉的是哪部分一万个天体放进二维方框做引力演化每一步要算多少对引力N(N-1)/2 ≈ 5×10⁷。单步还能忍一套轨道跑十万步就是 5×10¹² 次开方和除法一台工作站从天亮跑到天黑。Barnes-Hut 算法不改变物理它把「远处的天体团」折叠成一个质点让每个粒子的受力计算从一万次配对降到几百次整体复杂度从 O(N²) 压到 O(N log N)。这套 Barnes-Hut Simulator 的内核就三件事四叉树三维是八叉树、开角判据、蛙跳积分器。它拿的是「可接受的近似误差换十倍以上速度」这笔账。适合写星团演化、暗物质晕塌缩、分子动力学粗粒化、粒子群可视化的人不适合需要长时间精确轨道的双星系统也不适合 N1000 的场景——那种规模直接求和反而更快还少一层玄学。下面按「树怎么建 → 力怎么算 → 时间怎么走 → 哪里会崩」的顺序拆开代码用 Python NumPy 写重点在结构而不在语言换成 C 或 Rust 时一一对应即可。2. 从粒子表到四叉树Barnes-Hut 的节点里到底存什么2.1 一个节点必须携带的五样东西四叉树节点的最小字段集是子树总质量 mass、子树质心 (mx, my)、节点半边长 half、四个子节点指针 kids、以及叶节点持有的粒子下标数组 body。五个字段各自对应一件事少一个都要在后面付代价。mass 和 mx/my 是远场近似的全部依据。当判据认定某个节点「足够远」我们用它整体的质心和一个等效质量代替它内部所有粒子这一步等价于多极展开只保留零阶单极项误差量级由下一节的开角判据控制。half 提供空间尺度判据里的 s 就是 2×half也就是节点边长。kids 决定还能不能继续细分body 决定叶子上要老老实实做几次精确配对。这里最容易踩的坑是把「节点几何中心 (cx, cy)」当成质心用。几何中心是建树时切分出来的方框中心质心是质量加权平均两者在节点内质量分布不对称时可以差出半个边长。用错的那个力会有系统性偏差而且不会报错只会让轨道慢慢漂——这种血泪经验后文还会专门讲。import numpy as np class QuadNode: __slots__ (cx, cy, half, mass, mx, my, kids, body) def __init__(self, cx, cy, half): self.cx, self.cy, self.half cx, cy, half self.mass 0.0 # 子树总质量 self.mx, self.my 0.0, 0.0 # 子树质心质量加权不是几何中心 self.kids None # None 表示叶节点 self.body None # 叶节点持有的粒子下标数组__slots__在这里不是炫技N10⁵ 时节点数大约是粒子数的两到三倍每个节点省掉一个__dict__内存占用和缓存命中率都会明显不同。这一点在 Python 里只是小优化在 C 里换成紧凑数组就是数量级差别。2.2 自顶向下按象限分治别逐个插入建树有两种写法。一种是逐个粒子插入每遇到已占用的叶子就分裂另一种是先算全体包围盒再按象限把粒子下标数组递归切分。前者写起来顺手后者快得多而且在 Python 这种解释执行的环境里差距会被放大。原因在于下标数组的切分是一次向量化的布尔索引一次处理一整批粒子逐个插入则是 N 次指针跳跃每次都要走一遍从根到叶的路径。用分治写法建树整体是 O(N log N)递归深度约 log₄N。def build_tree(pos, mass, leaf_cap1, max_depth24): lo, hi pos.min(axis0), pos.max(axis0) center 0.5 * (lo hi) half 0.5 * float(np.max(hi - lo)) 1e-9 # 取最大跨度保证是正方形 root QuadNode(center[0], center[1], half) _build(root, np.arange(len(pos)), pos, mass, leaf_cap, max_depth) return root def _build(node, idx, pos, mass, leaf_cap, depth): w mass[idx] node.mass float(w.sum()) node.mx float((w * pos[idx, 0]).sum()) / node.mass node.my float((w * pos[idx, 1]).sum()) / node.mass # 两个终止条件粒子数不超过叶容量或已到深度上限重合点保护 if len(idx) leaf_cap or depth 0: node.body idx return rel pos[idx] - np.array([node.cx, node.cy]) # 象限编号x 正负 2 * y 正负 → 0..3对应左下/右下/左上/右上 quad (rel[:, 0] 0).astype(np.int8) 2 * (rel[:, 1] 0).astype(np.int8) h node.half * 0.5 offsets ((-h, -h), (h, -h), (-h, h), (h, h)) node.kids [] for q in range(4): sub idx[quad q] ox, oy offsets[q] child QuadNode(node.cx ox, node.cy oy, h) if len(sub): _build(child, sub, pos, mass, leaf_cap, depth - 1) node.kids.append(child)代码逻辑顺着读先把这个节点负责的那批粒子的总质量和质心算出来这是后面前序遍历力的前提再判断要不要继续分。要分就把粒子按所在象限打标签四个子节点各自拿走一批递归下去。参数上有三个可以调。leaf_cap是叶容量取 1 时树最深、精度最高、节点数最多取 4~8 可以显著减少节点总数和递归开销代价是叶子内部仍然要做几次精确配对——因为叶子上的粒子没法再近似了。实际项目里我一般先用leaf_cap1跑通确认精度对得上再往上调。max_depth24是重合点的后悔药如果两个粒子坐标完全一致象限切分永远分不开它们没有深度上限就会一路递归到栈溢出。取 24 意味着最小格子边长是初始边长的 2⁻²⁴早就低于双精度在这个尺度上的分辨能力了到此为止是合理的。还有一处容易忽略的细节根节点的half用的是max(hi - lo) / 2而不是分别取 x、y 跨度。因为四叉树的子节点是正方形根节点不是正方形的话分割出来的子节点边长会不一致开角判据里的 s 就失去了统一含义。2.3 质心必须由子节点回算不能靠插入时累加上面这份实现里质心是在每次_build入口用mass[idx]直接对这批粒子求加权平均的所以严格来说不存在「回填」这一步。如果你采用逐个插入的写法就必须在每次分裂时把父节点原有的质量重新分配到子节点并在整棵树建完后做一次后序遍历自底向上重算。两种写法的结果应当一致但逐个插入的写法在分裂时特别容易漏掉「把已有那个粒子的贡献搬到新子节点」这一步症状是父节点的 mass 对、质心偏越靠近树顶偏得越离谱而粒子少的地方又看不出来。所以我更倾向分治写法质心由定义直接算出来不依赖任何中间状态的正确性。代价是每个节点都要对它那一批粒子做一次加权求和整体是 O(N log N) 而不是 O(N)。这个代价在 N10⁵ 时大约是建树总时间的 30%~40%可以接受真要抠就在递归时把子节点已经算好的 mass 和 mx/my 传上来合并回到 O(N)。2.4 三条建完就该跑的自检断言树建错的典型表现不是崩溃是「跑得动但结果不对」。所以在进入力的计算之前加三条断言把它摁住def check_tree(root, pos, mass, tol1e-10): total_m root.mass assert abs(total_m - mass.sum()) tol * mass.sum(), 总质量不守恒 cx float((mass * pos[:, 0]).sum()) / mass.sum() cy float((mass * pos[:, 1]).sum()) / mass.sum() assert abs(root.mx - cx) tol and abs(root.my - cy) tol, 根节点质心不对 # 总质量为零的空节点不该被当成有效源 assert total_m 0.0, 全质量为零第一条查的是质量有没有在分裂过程中丢失或重复计入。第二条查的是「质心 vs 几何中心」这个最容易写错的地方——如果mx恰好等于root.cx八成是把几何中心赋过去了。第三条是为后面力计算的分母兜底mass 为零的节点在近似公式里会给出 0 加速度看起来无害但如果它同时又是某个粒子的唯一父节点说明这颗粒子根本没进树。3. 开角判据与受力遍历一次递归同时买到精度和速度3.1 s/d θ 到底近似掉了什么对一个距离为 d、边长为 s 的节点Barnes-Hut 用它的质心代替内部所有质量。这个近似展开到单极项就截止丢掉的是偶极项质心定义保证它为零和四极项及更高阶。四极项的相对大小大约是 (s/d)² 量级所以判据写成 s/d θ 时力的相对误差大致按 θ² 走。这就解释了两件事。第一θ 从 0.5 放宽到 1.0误差不是翻倍而是涨到四倍左右实测里经常看到「θ 调到 0.8 就明显不对了」。第二θ 趋近 0 时这个近似退化成精确配对算法复杂度回到 O(N²)——所以 θ 不是「越小越好」的精度旋钮它是精度和速度之间唯一的那根杠杆。3.2 单粒子受力遍历的完整实现遍历从根开始。对当前节点如果它是叶子就对叶内每个粒子做精确配对如果它是内部节点且满足开角判据就用它的 mass 和质心算一次近似力然后返回否则递归进四个子节点。整个函数是尾递归形式深度受max_depth约束。def accelerations(root, pos, mass, G1.0, theta0.5, eps0.01): acc np.zeros_like(pos) for i in range(len(pos)): _walk(root, i, pos, mass, G, theta, eps, acc) return acc def _walk(node, i, pos, mass, G, theta, eps, acc): if node.mass 0.0: return dx node.mx - pos[i, 0] dy node.my - pos[i, 1] # 叶子节点逐个粒子精确计算不做近似 if node.kids is None: for j in node.body: if j i: continue ex pos[j, 0] - pos[i, 0] ey pos[j, 1] - pos[i, 1] s2 ex * ex ey * ey eps * eps # 软化因子在这里进分母 inv G * mass[j] / (s2 * np.sqrt(s2)) acc[i, 0] ex * inv acc[i, 1] ey * inv return d2 dx * dx dy * dy d np.sqrt(d2) # 内部节点s 2 * half 是节点边长s/d theta 时整体近似 if d 0.0 and (2.0 * node.half) / d theta: s2 d2 eps * eps inv G * node.mass / (s2 * np.sqrt(s2)) acc[i, 0] dx * inv acc[i, 1] dy * inv return for child in node.kids: _walk(child, i, pos, mass, G, theta, eps, acc)逐段看逻辑。开头的node.mass 0.0是空节点早退用而不是是为了顺手挡掉负质量这种不可能但一旦出现就会静默污染全局的输入。叶子分支里j i的跳过不能省叶节点持有的是下标数组自身一定在里面不跳过就是自相互作用会得到一个巨大的虚假加速度。近似分支里的d 0.0是防除零d恰好为零意味着粒子和节点质心重合这在实际演化里出现概率极低但一旦出现就是 NaN 扩散。近似公式acc G * M * Δr / (d² ε²)^{3/2}是标准的 Plummer 软化形式。注意 ε 是加在 d² 上而不是加在 d 上也不是加在最终的力上——加的位置不同近距离行为完全不同这一点在第 5 章会展开。3.3 θ 的取值一张可以直接抄的对照表θ 取值力的相对误差量级单粒子访问节点数典型用途0.20.5%高接近精确配对需要长时间守恒的轨道、小规模标定0.51%~3%中星团演化、暗物质晕默认起点0.73%~6%低大尺度统计量、形态可视化1.06%~15%最低只看粗略结构的快速预览表里的误差是量级估计不是可以直接引用的数字——它和粒子分布、软化因子、节点形状都有关系。稳妥做法是拿自己的初始条件用直接求和下一节的脚本算一遍精确加速度再扫一遍 θ画一条「θ — 相对误差」曲线选误差落在你可接受范围内、且访问节点数已经明显下降的那个拐点。我一般从 0.5 起步因为这个值在多数分布下误差已经压到个位数百分比而节点访问量相比精确配对少了两个数量级。还有一点θ 是全局常数还是逐节点变化对结果影响不小。有些实现会根据节点里粒子数动态放宽 θ粒子多、距离远就更敢近似这在密度反差极大的分布里有效但它会让误差不再有统一量纲做精度标定时很难说清是哪个 θ 在起作用。前期建议先固定全局 θ把其他环节调稳再考虑动态策略。3.4 软化因子 ε 不是精度参数是物理参数ε 的作用是把 1/r² 的奇点抹平。在真实的星团里两颗星靠得足够近时会有非引力相互作用或者直接碰撞软化就是在粗粒化掉这个过程。所以 ε 的选择应该跟你的物理问题挂钩模拟的是稀疏星团ε 取平均间距的 1/10 到 1/100 就够模拟的是致密系统ε 太小会让近距离交会主导整个演化时间步被迫压得很小。经验做法是先估一下初始分布的平均最近邻距离把 ε 设在它的 1/10 附近然后跑短轨道看能量曲线的斜率。如果能量快速漂移先怀疑时间步再怀疑 ε如果双星轨道被异常散射那基本是 ε 太小。4. 时间积分把加速度接成能跑十万步的轨道4.1 为什么显式 Euler 在 N 体里必然翻车力的计算再准积分器选错也白搭。显式 Euler 每一步用当前加速度同时更新速度和位置局部误差 O(dt²)而且是系统性地往系统里注入能量——在引力这种长程、近距离交会频繁的势场里几万步之后能量偏移足以让整个系统「热膨胀」到面目全非。这不是精度不够的问题是方法本身没有守恒结构。蛙跳积分leapfrogKDK 形式属于辛积分器它精确保持相空间体积能量误差不在时间上单调增长而是在一个有界区间里振荡。这个性质对 N 体模拟是决定性的你不需要很小的时间步就能让总能量在十万步后仍然待在同一条带子里。代价是它只有二阶精度补偿手段是「半步 kick、整步 drift、再半步 kick」。4.2 KDK 三步和树的重建时机def leapfrog_step(pos, vel, mass, dt, acc, cfg): vel 0.5 * dt * acc # 半步 kick pos dt * vel # 整步 drift root build_tree(pos, mass, cfg[leaf_cap], cfg[max_depth]) acc accelerations(root, pos, mass, Gcfg[G], thetacfg[theta], epscfg[eps]) vel 0.5 * dt * acc # 再半步 kick return pos, vel, acc顺序不能换。第一步用旧位置的加速度把速度推半步第二步用这个中间速度把位置推一整步第三步在新位置重算加速度再补半步速度。这样速度和位置在时间上错开半步正好构成一个辛映射。注意build_tree和accelerations都必须在pos更新之后调用——树是位置的函数位置变了树就作废了。这一步是整段代码里最容易被「优化」错的地方有人为了省时间在多个步长内复用同一棵树结果粒子已经移出原来的节点范围近似判据算出来的 s/d 完全是错的误差不会报错只会让轨道莫名其妙地发散。acc在步与步之间传递省掉一次重复计算。第一次调用前需要单独算一次acc从那之后每一次的返回值都是下一步要用的旧加速度。4.3 能量漂移是唯一可信的黑匣子读数N 体模拟出问题时画面往往还「看着像那么回事」只有能量曲线会说实话。总能量 E KE PEKE 是 Σ½mv²PE 是所有粒子对的 -Gmᵢmⱼ/rᵢⱼ带同样的软化。def energy(pos, vel, mass, G1.0, eps0.01): ke 0.5 * np.sum(mass * np.sum(vel * vel, axis1)) # 势能用直接求和O(N²)只在 N 较小或做基准时用 dx pos[:, None, 0] - pos[None, :, 0] dy pos[:, None, 1] - pos[None, :, 1] r2 dx * dx dy * dy eps * eps np.fill_diagonal(r2, np.inf) # 排除自相互作用 pe -0.5 * G * np.sum(mass[:, None] * mass[None, :] / np.sqrt(r2)) return ke pe逻辑上有两个点要抠。np.fill_diagonal(r2, np.inf)把自配对的距离置成无穷大对应的贡献是 0同时避免了对角线上的除零。前面的 0.5 是因为mass[:,None] * mass[None,:]这张 N×N 矩阵把每一对算了两遍。软化因子必须和力的计算里用的完全一致。如果力里用(d²ε²)^{3/2}、势能里用1/d能量曲线会呈现出「看起来在漂移」的假象实际上是两套势场在打架这种不一致比数值误差更难查因为它不随步长减小而收敛。判断标准很简单把 dt 减半再跑一遍如果能量漂移的形状和幅度都明显改变说明是时间步问题如果几乎不变说明误差来自树的近似θ或者积分器的守恒性质跟 dt 无关。这一步判断能省掉大量盲目调参。4.4 自适应时间步让近距离交会自己决定节奏固定 dt 的问题在于它必须迁就最极端的那个事件。一套星团里大部分粒子运动平缓个别粒子偶尔会有近距离交会用统一的小步长跑全程是巨大的浪费。常见做法是 Aarseth 判据的变体每个粒子按自己的加速度算一个特征时间 τᵢ η·sqrt(ε/|aᵢ|)取全局最小值作为这一步的 dt。η 一般取 0.01~0.05它控制的是「一步之内速度变化占特征速度的比例」。这么做能自动在系统剧烈变化时收缩步长平静时放大。代价是 dt 变成时变的严格的辛性质会打折扣——不过只要 dt 变化缓慢能量误差仍然能保持在有界振荡而不是单调漂移这一点在实践中比理论上的严格性更重要。如果观察到能量曲线变成阶梯状上跳多半是 dt 变化太剧烈把 η 调小或者对 dt 做平滑比如新 dt 不超过旧 dt 的 1.2 倍就能压住。5. 排查手册θ、ε 和退化树导致的五类典型故障5.1 力算错的三类症状以及该先看哪里Barnes-Hut 的 bug 有一个共同特点不崩溃。程序照跑图照出只有和基准对不上。所以排查顺序要固定下来别东一榔头西一棒子。第一顺位是能量漂移的形状。单调上升或下降查积分器和 dt有界振荡但幅值偏大查 θ 和 ε出现台阶状突变查自适应步长。第二顺位是加速度的逐粒子对比。拿 N200 的小规模用直接求和算出精确加速度和树方法逐分量相减看误差是均匀分布还是集中在某些粒子上。均匀分布是 θ 的问题集中在少数粒子通常是靠近节点边界的是判据实现的问题集中在质量最大的那几个粒子上是质心算错了。第三顺位才是看轨迹图。因为轨迹是最后才被误差污染的东西从它反推原因效率很低。5.2 五条翻车记录一、跑几十步后位置全是 NaN现象pos里出现nan能量变成-inf日志里没有任何异常。原因初始条件里有两个粒子坐标完全相同或者近到双精度在这个尺度上分辨不出。象限判据用的是 0两个点落在同一侧永远分不开递归一路到底直到栈溢出即使有深度上限兜住到了力的计算阶段 d≈0 又会除零。解决建树时加max_depth兜底本文实现里已有力计算里所有分母都过一遍软化因子 ε初始化生成粒子时对重复坐标加 1e-10 量级的随机扰动。三条里最不该省的是第一条因为它是唯一能防住栈溢出的。二、同一个 θ别人的代码误差 1%我的误差 8%现象对照直接求和相对误差差了好几倍但代码逻辑逐行比对看不出区别。原因θ 的比较方式没有统一约定。常见写法有三种s/d θ、s²/d² θ²等价、以及用节点对角线代替边长差一个 √2。更隐蔽的是 s 取边长还是半边长——差一个 2 的因子等效于 θ 差了 2 倍误差按 θ² 走就是 4 倍差距。解决在代码里把约定写死并注释清楚本文用的是「s 2 × half 即边长判据 s/d θ」然后用同一份初始条件、同一套直接求和基准扫一遍 θ 画出自己的误差曲线。别拿别人的 θ 数字直接用。三、总质量守恒但根节点的质心慢慢偏离现象root.mass一直对root.mx和全体粒子的加权平均对不上偏差随时间增大。原因用了节点的几何中心 (cx, cy) 代替质心做远场近似。几何中心在节点内质量分布对称时恰好等于质心所以小规模测试很容易蒙混过关一旦分布不对称偏差立刻显现而且它表现为力的系统性偏移不会报错。解决每个节点的 mx/my 必须由(mass * pos).sum() / mass.sum()算出来并且加一条断言根节点的 mx 与全体粒子的加权平均之差应当在 1e-10 量级。如果这两个数恰好相等而root.cx也相等反而要警惕——那说明初始分布太对称没测出问题。四、双星靠得很近时被弹飞能量突然跳一截现象两条轨道接近到最小距离时其中一个粒子的速度异常增大之后整团天体的能量曲线抬升一个台阶。原因软化因子 ε 只加在了力的大小上势能里用的还是1/r或者力和势能虽然都加了 ε但形式不同一个用(d²ε²)^{3/2}一个用d²ε²。两套势场不一致在近距离时表现出的等效能量不守恒。解决力和势能统一使用 Plummer 形式——力GmΔr/(d²ε²)^{3/2}势能-Gm/sqrt(d²ε²)。改完之后能量曲线的台阶应当消失剩下的平滑漂移才是数值误差。五、总能量单调上升系统整体「变热」现象几万步之后原本的团块逐渐散开速度弥散变大能量曲线一路向上。原因多半是 dt 太大遇上近距离交会。也有可能是 θ 太松导致树近似在每步引入的误差方向不是随机的而是有偏的近似总是低估远场的力长期累积就成了净能量注入。解决先把 θ 从 0.5 收到 0.3 重跑如果漂移速度明显下降问题在树近似如果几乎没变把 dt 减半再跑漂移幅度应随之减半二阶方法。两个都试过还漂检查自适应步长是否变化太剧烈。5.3 一个三十行的交叉验证脚本上面所有排查都指向同一件事你得有一个「随时能拿来对答案」的精确解。直接求和写起来很短值得常备。def direct_accelerations(pos, mass, G1.0, eps0.01): dx pos[:, None, 0] - pos[None, :, 0] # dx[i, j] x_i - x_j dy pos[:, None, 1] - pos[None, :, 1] r2 dx * dx dy * dy eps * eps np.fill_diagonal(r2, np.inf) # 排除自配对 inv G / (r2 * np.sqrt(r2)) # 1/(r²ε²)^{3/2} ax -np.sum(mass[None, :] * dx * inv, axis1) ay -np.sum(mass[None, :] * dy * inv, axis1) return np.stack([ax, ay], axis1) def compare(pos, mass, root, theta, eps): a_tree accelerations(root, pos, mass, thetatheta, epseps) a_ref direct_accelerations(pos, mass, epseps) rel np.linalg.norm(a_tree - a_ref) / np.linalg.norm(a_ref) worst np.max(np.linalg.norm(a_tree - a_ref, axis1) / np.linalg.norm(a_ref, axis1)) return rel, worst符号要留意dx[i, j] x_i - x_j而作用在 i 上的加速度指向 j所以整体带一个负号。这个负号和前面树遍历里的写法是等价的那边写的是pos[j] - pos[i]不带负号两种约定都对但混用会让误差看起来像 200% 而不是 2%。worst这个指标比全局相对误差更有用。全局误差被大量平缓的远场粒子摊薄了看起来总是很漂亮worst 会告诉你最差的那个粒子错了多少而轨道演化恰恰是被最差的那几个粒子带偏的。我一般要求 θ0.5 时 worst 控制在 5% 以内。6. 进阶给质心项补上四极矩以及我踩过的那个坑单极近似丢掉的是四极项。如果你想在同样的 θ 下把误差压下来最划算的改法不是调小 θ那是指数级的代价而是给每个节点额外存两个二阶矩分量力公式里补一项项存什么计算代价收益单极mass, mx, my已有基准四极Qxx, Qyy, Qxy建树时各多两次加权求和θ 相同下误差降一个量级八极三阶矩明显变重边际收益低不推荐四极修正的形状是把G·M·Δr/(d²ε²)^{3/2}换成G·M·Δr/(d²ε²)^{3/2} 修正项(Δr, Q)。修正项里包含Δr的三次方组合和 Q 的分量具体形式按二维展开写开有六项实践中多数人直接调用现成的多极展开库或者查表实现没必要手推。我更想讲的是另一样东西在 Python 里用对象指针建树N 超过 5×10⁴ 之后耗时的大头不是浮动运算而是属性访问和缓存未命中。那次翻车是把一份能跑 N10⁴ 的代码直接拿去跑 10⁵以为复杂度是 O(N log N) 所以只慢十倍结果慢了将近四十倍而且 profiler 显示 70% 的时间花在_walk的属性读取上。改法是把树压成几个并行的 NumPy 数组mass[]、mx[]、my[]、half[]、first_child[]、child_count[]本文还有配套的精品资源点击获取
