From 81f7bd84d0dd8bc25b9034959fb5fe0380bae4d9 Mon Sep 17 00:00:00 2001 From: Owen Rummage Date: Fri, 17 Jul 2026 13:34:44 -0500 Subject: [PATCH] initial commit for alpha release and testing on MacOS --- Makefile | 152 ++++++ README.md | 124 +++++ dist/fossmark-linux-amd64 | Bin 0 -> 30440 bytes dist/fossmark-linux-arm64 | Bin 0 -> 73872 bytes dist/test_kernels | Bin 0 -> 26096 bytes src/fossmark.S | 942 ++++++++++++++++++++++++++++++++++++ src/fossmark_x86_64.S | 987 ++++++++++++++++++++++++++++++++++++++ src/main.c | 709 +++++++++++++++++++++++++++ src/test_kernels.c | 462 ++++++++++++++++++ 9 files changed, 3376 insertions(+) create mode 100644 Makefile create mode 100644 README.md create mode 100755 dist/fossmark-linux-amd64 create mode 100755 dist/fossmark-linux-arm64 create mode 100755 dist/test_kernels create mode 100644 src/fossmark.S create mode 100644 src/fossmark_x86_64.S create mode 100644 src/main.c create mode 100644 src/test_kernels.c diff --git a/Makefile b/Makefile new file mode 100644 index 0000000..e74f4f3 --- /dev/null +++ b/Makefile @@ -0,0 +1,152 @@ +# fossmark - multi-core CPU benchmark +# +# The assembly kernels are architecture-specific: +# src/fossmark.S AArch64 (ARM64) +# src/fossmark_x86_64.S x86-64 (AMD64) +# The C driver (src/main.c) is portable across architectures and OSes. A +# "binary that runs everywhere" is not possible - each OS/arch pair uses a +# different executable format and instruction set - so output is named per +# platform, e.g. dist/fossmark-linux-arm64, dist/fossmark-linux-amd64. +# +# Common targets: +# make build for the host arch (dist/fossmark--) +# make linux-arm64 build the Linux/ARM64 binary +# make linux-amd64 build the Linux/AMD64 binary +# make macos-arm64 build the macOS/ARM64 binary +# make macos-amd64 build the macOS/AMD64 binary +# make all build both Linux binaries +# make bench build for the host and run it +# make test build and run the kernel correctness tests (host arch) +# make clean remove dist/ +# +# Cross-compiling: linux-amd64 on an ARM64 host (or vice versa) needs the +# matching cross toolchain. The compiler for each target defaults to the host +# `cc` when the host arch already matches, and to the conventional GNU cross +# compiler otherwise. Override with CC_ARM64=... / CC_AMD64=... if your +# toolchain is named differently, e.g.: +# make linux-amd64 CC_AMD64=x86_64-linux-gnu-gcc-14 +# make linux-arm64 CC_ARM64="clang --target=aarch64-linux-gnu" +# +# On macOS, Apple Clang can build both architectures. The macOS compiler may +# be overridden for an osxcross or other cross toolchain: +# make macos-arm64 CC_MACOS_ARM64=clang +# make macos-amd64 CC_MACOS_AMD64=clang + +CC ?= cc +CFLAGS ?= -O2 -Wall -Wextra +LDLIBS ?= -lm +# The driver spreads each workload across all cores with pthreads. +PTHREAD := -pthread + +DIST := dist +DRIVER := src/main.c +ASM_ARM64 := src/fossmark.S +ASM_AMD64 := src/fossmark_x86_64.S + +# ---- host detection: normalise `uname -m` to our arch names ---- +HOST_ARCH := $(shell uname -m) +ifneq (,$(filter aarch64 arm64,$(HOST_ARCH))) + HOST_ARCHNAME := arm64 + HOST_ASM := $(ASM_ARM64) +else ifneq (,$(filter x86_64 amd64,$(HOST_ARCH))) + HOST_ARCHNAME := amd64 + HOST_ASM := $(ASM_AMD64) +else + HOST_ARCHNAME := $(HOST_ARCH) + HOST_ASM := $(ASM_ARM64) +endif + +# ---- host OS name for the native binary ---- +UNAME_S := $(shell uname -s) +ifeq ($(UNAME_S),Linux) + OSNAME := linux +else ifeq ($(UNAME_S),Darwin) + OSNAME := macos +else ifeq ($(OS),Windows_NT) + OSNAME := windows +else + OSNAME := $(shell uname -s | tr '[:upper:]' '[:lower:]') +endif + +# ---- per-target compilers: native cc if the host matches, else a cross gcc ---- +ifeq ($(HOST_ARCHNAME),arm64) + CC_ARM64 ?= $(CC) +else + CC_ARM64 ?= aarch64-linux-gnu-gcc +endif +ifeq ($(HOST_ARCHNAME),amd64) + CC_AMD64 ?= $(CC) +else + CC_AMD64 ?= x86_64-linux-gnu-gcc +endif +CC_MACOS_ARM64 ?= $(CC) +CC_MACOS_AMD64 ?= $(CC) + +NATIVE_BIN := $(DIST)/fossmark-$(OSNAME)-$(HOST_ARCHNAME) + +# `make` with no target builds the host binary, as before. +.DEFAULT_GOAL := native +.PHONY: all native linux-arm64 linux-amd64 macos-arm64 macos-amd64 bench test clean + +# `make all` builds both Linux binaries. +all: linux-arm64 linux-amd64 + +# `make native` (and bare `make`) build for whatever host you are on. +native: $(NATIVE_BIN) + +linux-arm64: $(DIST)/fossmark-linux-arm64 +linux-amd64: $(DIST)/fossmark-linux-amd64 +macos-arm64: $(DIST)/fossmark-macos-arm64 +macos-amd64: $(DIST)/fossmark-macos-amd64 + +$(DIST)/fossmark-linux-arm64: $(DRIVER) $(ASM_ARM64) | $(DIST) + $(CC_ARM64) $(CFLAGS) $(PTHREAD) -o $@ $(DRIVER) $(ASM_ARM64) $(LDLIBS) + @echo "built $@" + +$(DIST)/fossmark-linux-amd64: $(DRIVER) $(ASM_AMD64) | $(DIST) + $(CC_AMD64) $(CFLAGS) $(PTHREAD) -o $@ $(DRIVER) $(ASM_AMD64) $(LDLIBS) + @echo "built $@" + +$(DIST)/fossmark-macos-arm64: $(DRIVER) $(ASM_ARM64) | $(DIST) + $(CC_MACOS_ARM64) -arch arm64 $(CFLAGS) $(PTHREAD) -o $@ $(DRIVER) $(ASM_ARM64) $(LDLIBS) + @echo "built $@" + +$(DIST)/fossmark-macos-amd64: $(DRIVER) $(ASM_AMD64) | $(DIST) + $(CC_MACOS_AMD64) -arch x86_64 $(CFLAGS) $(PTHREAD) -o $@ $(DRIVER) $(ASM_AMD64) $(LDLIBS) + @echo "built $@" + +# When the host is Linux/ARM64 or Linux/AMD64, the native binary IS one of the +# linux-* targets above, so no separate recipe is defined (that would be a +# duplicate). Otherwise - e.g. macOS/ARM64 - provide the native recipe here. +ifeq ($(OSNAME)-$(HOST_ARCHNAME),linux-arm64) +NATIVE_HAS_RULE := yes +endif +ifeq ($(OSNAME)-$(HOST_ARCHNAME),linux-amd64) +NATIVE_HAS_RULE := yes +endif +ifeq ($(OSNAME)-$(HOST_ARCHNAME),macos-arm64) +NATIVE_HAS_RULE := yes +endif +ifeq ($(OSNAME)-$(HOST_ARCHNAME),macos-amd64) +NATIVE_HAS_RULE := yes +endif +ifneq ($(NATIVE_HAS_RULE),yes) +$(NATIVE_BIN): $(DRIVER) $(HOST_ASM) | $(DIST) + $(CC) $(CFLAGS) $(PTHREAD) -o $@ $(DRIVER) $(HOST_ASM) $(LDLIBS) + @echo "built $@" +endif + +$(DIST): + mkdir -p $(DIST) + +# Build for the host and run the benchmark. +bench: $(NATIVE_BIN) + ./$(NATIVE_BIN) + +# Build and run the kernel correctness tests for the host arch. +test: | $(DIST) + $(CC) $(CFLAGS) $(PTHREAD) -o $(DIST)/test_kernels src/test_kernels.c $(HOST_ASM) $(LDLIBS) + ./$(DIST)/test_kernels + +clean: + rm -rf $(DIST) diff --git a/README.md b/README.md new file mode 100644 index 0000000..b0e0c20 --- /dev/null +++ b/README.md @@ -0,0 +1,124 @@ +# Fossmark + +A single-threaded CPU benchmark for ARM64 (AArch64), with the numeric kernels +hand-written in assembly and a small portable C driver to run and score them. + +## What it measures + +Nine workloads, each a tight assembly kernel: + +| # | Test | What it exercises | +|---|-------------------------|--------------------------------------------------------------| +| 1 | Integer Math | 64-bit ALU: `madd`, `umulh`/`smulh`, `udiv`/`sdiv`, bit ops | +| 2 | Floating Point Math | scalar double: `fmadd`, `fdiv`, `fsqrt` | +| 3 | Prime Numbers | sieve of Eratosthenes to 2,000,000 (strided memory + ALU) | +| 4 | Extended Instructions | NEON/ASIMD: 128-bit integer, widening, table, float vectors | +| 5 | Compression | LZ77 match-finder over a 4 MiB corpus (branchy, cache probe) | +| 6 | Encryption | ChaCha20, 20 rounds, NEON, over 1 MiB | +| 7 | Physics | 512-body direct-summation gravity, double precision | +| 8 | Sorting | in-place heapsort of 1M `uint32` (branch + cache stress) | +| 9 | Single-Threaded | dependent-load pointer chase over 16 MiB (memory latency) | + +Each test auto-calibrates its iteration count until it runs long enough to be +timed reliably, then reports the best of several runs (the run least disturbed +by the OS scheduler). Every kernel returns a checksum that the driver verifies +across runs, so a miscompiled or non-deterministic kernel is caught rather than +silently mis-scored. + +## Scoring + +Each test's raw rate is normalised against a **reference machine** into a +unitless score, and the overall is a **weighted geometric mean** of those +scores: + +``` +S_i = TARGET * (rate_i / REF_i) (per-test score) +Overall = TARGET * exp( Σ w_i·ln(rate_i/REF_i) / Σ w_i ) (weighted geo. mean) +``` + +The reference rates are the tuning machine's own rates, and `TARGET` is 20000, +so that machine scores ~20000 on every test and overall. Scaling is linear in +performance: a machine half as fast scores ~10000, one 10× slower ~2000, and a +future machine twice as fast ~40000 — so there is unbounded room both below and +above the reference. + +The weights reflect each test's influence on **everyday, common-workload user +experience** — integer/general-purpose throughput and memory-latency-bound +responsiveness matter most; specialised floating-point and physics matter +least. This mirrors the weighted, integer-dominant approach of mainstream +suites such as Geekbench 6 (which splits integer/FP roughly 65/35 and combines +real-world workloads with a weighted mean). + +| Test | Weight | +|---|---:| +| Integer Math | 20% | +| Single-Threaded | 16% | +| Compression | 14% | +| Sorting | 12% | +| Extended Instructions | 11% | +| Floating Point | 9% | +| Encryption | 8% | +| Prime Numbers | 6% | +| Physics | 4% | + +Everything above is configurable via `#define`s at the top of `src/main.c`: +`FM_TARGET_SCORE`, the nine `FM_REF_*` reference rates, and the nine +`FM_WEIGHT_*` weights. Weights are relative — the code normalises by their sum, +so you can change one without rebalancing the rest. To re-baseline for a +different reference machine, set each `FM_REF_*` to that machine's measured +rate. + +Note: the pointer-chase (Single-Threaded) test measures raw memory latency and +is the noisiest to sample, so the overall typically varies ~1–2% run to run. + +## "Runs on all operating systems" + +The **assembly is** OS-independent: `src/fossmark.S` contains no system calls, +no libc calls, and no external relocations. Every routine is a pure function of +its arguments under the AAPCS64 calling convention, so the same source +assembles and runs correctly on Linux (ELF), macOS (Mach-O), Windows (COFF) and +the BSDs. It avoids `x18` (reserved on Darwin/Windows) and the `v8`–`v15` +callee-saved vector bank. + +A single *binary* that runs everywhere is not possible — Linux, macOS and +Windows use incompatible executable formats and system-call ABIs. So the +portable C driver (`src/main.c`) supplies the per-OS parts (timing, memory, +I/O), and you build one binary per platform. The Linux build is named +`fossmark-linux-arm64`. + +## Build + +```sh +make # builds dist/fossmark-- for the host +make linux-arm64 +make linux-amd64 +make macos-arm64 +make macos-amd64 +make bench # build and run the benchmark +make test # build and run the kernel correctness tests +``` + +Both macOS targets can be built on either Apple Silicon or Intel Macs; Apple +Clang selects the requested architecture with `-arch`. They produce +`dist/fossmark-macos-arm64` and `dist/fossmark-macos-amd64`, respectively. + +Or by hand: + +```sh +cc -O2 src/main.c src/fossmark.S -o dist/fossmark-linux-arm64 -lm +``` + +On macOS the same command produces a native binary (name it +`fossmark-macos-arm64`); on Windows use `clang` from the LLVM/MSVC toolchain. + +## Testing + +`src/test_kernels.c` is a standalone harness that validates each kernel against +an independent reference or invariant — the sieve against a C reference sieve, +the NEON ChaCha20 against a scalar reference anchored to the RFC 8439 +known-answer vector, the sort against `qsort`, the N-body step against +conservation of momentum, and so on. It exits non-zero if any check fails. + +```sh +make test +``` diff --git a/dist/fossmark-linux-amd64 b/dist/fossmark-linux-amd64 new file mode 100755 index 0000000000000000000000000000000000000000..94308bd71e13e6c628b40a89e41ccd800fc9cc6a GIT binary patch literal 30440 zcmeHwdw5jUx%Zx2AZo}SENB#wtxarD5;6#wC@C|!z}`B+RDwbghLFsVkt8#nnIL$n zsUeo#ag?@Jdp!MGo__7A?Ppu;=>@2T2{(ay38;Xg81c44iiDsN6`AvU*IsKf*-hH# z>v_KKIe)NW_PgHSdf#`g^{(q)Gg@m>KUh31ZDlO&?0LI8HF>JXC0Ye%*}Uc-K&xMUl8Ra4PiDfDsg5Lb0ID zLunrs1ec(~pOi=XtmFQ)(mYbaB`C?#MU*}Ri-`Ms>SQR@xp`UBE8*p{QiaP4D$<+( zUZmk)Nw0?MExX+`Y^7^?dV-4dYPnt`FQ1h@%=H8n<@*}+MAdMSe@J*UlWybbja2V7 ztPx6*pu9R+9gIT?RDE;jEs%b(ar_(8ZulsDZr4)xH8*_!7me3$a87OVFP}MmYEwgD zlfR{Pb>ZrgnT0c_7li{w(v!@YuIX;soK72)Ke^hiQAVaTvz+OBQ zVjTH%IQ%ao;FPR_kM)`dU^qLAN5H>70`4CHKQsbdqPXw;uW4%5C zFkE@5oehWoVg!8t2>1gd;AyCU;p|U=yn>JQdJ@2J@<)MRZ50MTX9W2_0v{_)km5m8 zN%W(aINoi+C688H?P>6ayes_Sh&NPQUC|V1@z&NYZ}LhW&x+BC z@9Ln`6j&hzTO(m9ye3>9Xz@w)O@aEAo)z9m#NSM|UobTEmS}I1*~|wn)K9W5`?A;JGW{A7ttQMZ6Gc2((6E6Fsq>=dZaAa)|JT zU?~{#w?upps9#;@@%dZon*8^8rDkt)*c%~IA6!G})HesEU?A*Y?Ewav6)h4hHd9xg zNuhc(G}o<_S9qosO)Hf|&@nw391J;bNP^c7qSK0IP+g?qAB`$$7-P~f2NFMu-6d%> zXgUki_(35NW5Y-?7Hx^fH3c92yN_d7871WdqssuFkFogCQZdKI`Tb*P#=VL2 zPq6q5sfzQzXMDO;!}-?B{TEOk7lWru%$K4a1*e<%CO(e}T;X`!f={S5<@Z?dNgVI8 z;Q1UsY{9SNc((2oZB}nxUerZndyu6 z@EW(DWx+q-xZQ#~x%~nQp1#1e>$Kn(aa^_FS8=??f){i9jTT(wzs-Wr=62Ru@WaCY z7JTCj?tcq@3%3)u;58iIW5E}5yvu^G=k^T?F6>Ku-7WmTirdMu;NRr9-GbL}`vn&K zNiOfS;6LTKYQcZS@fr&rpKki6(SmQ~@@*D;FUQwe@B!g}3%-Hdf69U%=5{t%=5xVs z44C<%H7I`+VZf#oT&&G0BA$ZhFiG>aCj}p$f_J6hm!{xi-9gK1TJwratUCy&b#U_H z`5z5{PdO|qD+Q-HQC#*EoG;;7RBj5+*4-%aTTQCqHR%JaSD!!FnN`v;1?!Q zNphy(7p3583Z9*U*QDTMQt*W-xSWD7OTq0acw-7aHU$r+;1{RhSfH_seDI|dye$PE zmx8ZL!7oX{*QemQDfs!|c@LcTz%x7+uB()4`J`cf`|HP|rrj|B86p0_xeo_@5F_TPi*48UG~lw3X7+#`y0M zPumAQLB@ZFc-lJYS;qJWiKi`-o*KsALp*Jj^f(zGBA&KLdWso;7xA<;(xWhbIq|e5 z(qm`**NLaCs2+*&w-GNB-**l`{>{YOiSK57Iq|d=($mHG8;Peake)br(0_QT-}>7^ zI78Kclvlv$(|N_9ZvAI@C5%3u=VbKRyc*C3cQh}*02WCgmm=yv(Av-1ZX4`{%GWVZC1xtCmGN)`>P?;%yV-Wa-&A?%v+{*9Ff=FNgPxcLe(yDg^6Dy ztFyQQ6Ny&lEmZZ&JOvq%*J{iLaZs(G=I7CHSDy5an(t4ajxs;ZZl~AM=BGn&o)$cJmcb-%2*ekDHOPt=6 zSL`^UZQJO3{q#?Egukr)W%GOQUeP>v;WkwtQwU{uY(~XYRL--#V;x`mv)1*Ps=M)Fe)X4|k(dm6|{ft{LTjECMbI&Fc z?pQ{F7E^8Nw?3unlfT{0YHhn)|3r<=x{k7x<%T=##@9(2-m$AuVCtsv_uc%!<~29` zSc{Dtt?CDngDX|XLFDISRqr!%*Eb}0eeK<&)%F3~W%2gzx0r=B{X^5`3Nq!M+<$Xy*7Yo&Wy@Iixv7H;ex_{}`B9_yX z&%*6D$EIJX*}AkX!rwk$k?4N^G^pyQ)%NFXYG*G>dD?Az5k)(q>W9sO9tMZ*u>;Yl zhQsUTCbZJN*1w`VHr#?5-XgM~Z5C;(+eCidv8BrTf117<#@lzZ^xcjfnr#3nzCI-V z&Xn}$ApO%wf9oI4^-TWm&$rUJp>6uZmwXYsr=(=iC4;_DZ3oFWNiVWSn3Ai1t^C_I zJTc+p+uxI()DT+v$8)d!aE7);6kOY!bjzlsTWp6k+o54<#ceyJws&Tz`Z=?qokRIg z4QXicWH;#3U1!~;iPpdN%l3(N5h2bpy11c!vr5b0B$#jke4+*u5}udRQ-_V zh^w6^Vd#)s-wy2zNPR}CXWmiu0c2~E>OathnO&RQr8=I|+B==m)b$oOF zifY8_Y4Kb>`ePKl17R9ZvfL}XOK8w9aO=m6&NGRGOMi}Lgnq17Xu4`-)787g=@U1T zR{MQ<1(N*8Mig(b)i@HcQ9L-%rw;A+(en={paFwgL3V}?;&QvP0QYKa;vz~O0ePSk)#&rm!vfG*B->_n% zInbUmuR?Pj6RNNqKjYeeft)*PMpSnsLDGQC2t4SGh9>uAjOS_ZPjN(j&^ z*hId1aB?+fbA?$X)en(9{D(t@X3KHiJ4K6(WhXKnM11?O;L6P=GrNwP)MOrF@=gLn8~ zu7O8l6Tk9RxZs4**PlpiqbCt}>}Po!8GSl$86J&{-H3tsrC+AF8Mzn(jU6U59gy)S z79RgfK3SVV1Vd@8Y(E0V5sHh>_**%}l~qF=ae6<~XNehp)6zX;sTpAtosW<^Hs&hi z&=GgXGM`|^$NYei@!w};#&;PhdlbY7nP$Rd<{An|%N|9d#ukJt_51Q>saV|GQ7VZ7 ziU_bTXM37#-#wbY-tn*2zqNG|MaNVKC#&qPtawcUbR+fWIt2#z& zq>*bpbt;iq_;qG6OKsnUwz$KoweKv}+IOKFmu)8NB$;I#g`_D*b1FjPpNGF_0zf!C z&@XR)J`7wv}FIVP$=lwnX%)b=kfkk`(Juhr=F zDE3DvHlCt3^|9J|+lkuBhA}l(*(^oJZw9aGN7c^m3{6L;d^UsHq;W0A>ziZG&Z35C zKAh=Cjqw1fSa1pvEOXq^dkfs9VY}Qh8MB^-4#1w8)~$jXoOOiePWYhUWqdkO*SwM! z%~z1IdSnb4zca5GchoEHSnK!8OSeP%6%;V@X=-wblV}1F)&y8IH39elY8{`*Yh6ee z+)71@Rw-x|J+sX`qE*O?-h?r>E~|$j@UPL(zK))TSxq8yaiF_cRl(g$_kbI717v%C zOV6CKvZx|S)Qx8urN{-+VHL6zP24CTf5v9~`D!#U7^DWa4Pm232%=J$Amtx|ze^E( zoCLpO3TFO^ny6`E#=)zYp`A$97>#(dm;U-gScqokplI%B{x3jKj&IhfE8`_@{hga* zQxw{SNgcnUH_@p@k0e{MZ27eqi^2kK$7_M#Pd*e0awdW2z7gk{)l~)b2%yEh&0mXOn#I%ODw2@AGWs)w_$WZ+f=^v`f8z z684R8Ns07@` za+diSBGl4>$au7u4{t&kdp8#(y5&i2D5SK5Q^@Zp!Pty_L?)kOGwX-5)y{!*)wWyR ze=af(2L8egWF@*W6;}@GQ~uh@{sLX8_4);Bd)bqagje5-jDeGu0#A%z2gs;LW_M>k z2Eevt$6$Vuk0I@q?DedGEYIJ3RsUR#-i?R<2eJ}h!lEtL5nr65PucaKahA&%;#L>; z#IVcdOI`0`D!}5pfvr6vAyvoJQRR5ofk)^jtUUFP-Jpu5f00~|MDNE|z>G}flhq>a zi5Z7_sa4-!u$3y!Sc6GTbL`a6FZR3jZd#{fwt+b}ZWwqdvG$z>C>i{O#2_;WE>QI^ zSQA3Q@MKT()Hf9*_#^UBfTeRU#_gS1Xqbz)N@vF(;-@ly@!`C8}+g zR{BQkTdHob60#z&60)-L5@H#{1jw0G=Jl#SvEd;bDq~e{PZyTL{Ip8j=oy7& zgR5C z!ru*Q^a?gekt++-`q%Z^ybQJTolI4K-H?yd6Wi`8%sw5vFm{l8a z25|eQHu=$;5O6(9RVa7djj$RWZCvv(|YDY#E@VdKVK`g!E~O1;m4ZME^g#Zpw& zU(;fd3t(oJ-Lk1P5^-^l2 zN=gO1D|bA`W`g6r*O}v%C>Lppvf-IJ!xowBIMTBi0}m<|9pc>^kd@esO#FrvE2iDZ ztU0`ABPFV{o=YHVECoKa>RflYZrAj4tUY0Pp(drO(WSX58}ki$B`*Dy zmtyfcb>H3Dh#B)fLze8rcE##4v7o}jeGxKgXY!Ek3}1{aO*u?GyC;JgcEawRohB-8 z`FwUSc76`Gzr&|=L%rZI6(UYkdhauAJjJ6#HU>r1jVLnQfTfu6pYPEM0}g|C;63t_ z8_PAt+(aHvEkCJ(ht53XK~lJBN}k+t0eTVh5_2)0IcR4ND^{!kPh5^71Fp?eXuIhx z6o9uG>q-X@^JxH#*4;3&VX_IAeBR%rzF{VTEBsgqq_CJ zAsxp_U$A7Ocv8<_(GSp6Y3o*hoid^BH`vl7KdQ=C zx1QYk>!VPR-12XCSKIc0yjX494|1^D_BzPXYTFJsdxoEgEkZZ#b$5@|9C&EJZr`3s zngd@FVBc>&2#VIr_;Sv!6Wa_y7vnD{M5YbW|y zDLiknqI*$vB#(YeEqdX0pDKSdo>G3u;u^p+09_eQMx<;`yE%S}(;7}C9^W?`=IMHX z-9b%#sZCI|J&`6q_#Lolqj8n?EMz3|29q{O1b>Uaii8AXuVj|JGIhJ!xxtC&d9M<+d6ws zY0!Md5!asUBqfTJ%iCtZTI0`u&+-LBxav=}#^ z6StNc+*U2a?VG_$d2u&BJ-#aWWn!R@`hU{CHR_5ebm{&b%lIxrEzUJO`A-tfsE8b0~> zDWBYZmV`+0lQ$9l#XCMZ7@q)-?x~Us&Q{5bzNnHHo~e@O@2rwD*lE~ZQFjNbzTZ9hm`ktikhIwJYfqg@xH$#ohvAH~0$Y4RUcnQ8tdmo#` zXVHMzH@2q{!2W6e$7YMV7HuYEhmR7}UQSdy;?)j(UlG9GP~Zi84#M#*tUXZk2<$bR z0zyYzTWBNEtb~(H;xrn^X*(1`2kG96J=Aklqx(n+Z6tImn?;VJ$QS=kGmAy$4zD?y zMP9}=gL_NS1P-7r1dn|G?+pEa$^sJRStxZS_Q@TmFyA46=i+puj&qLs%(){s$1>#8 z@PGn0oa{#E0fpFQcH-YAewm#COyC#;Oo?KWc!vQjXb|6v@q@mZ(n?Rs53-*mB=*W3 zbIinpyXHaWoPmUcyPh$FmGU=@l<^dQpei50t%w0NMveIZThE+_CJ@82Rlsx@l}r?rUgAcO1dzU3}1#4ru!46nj10 z?#;qvVXv^${T8e@?(D|xv%(T%YS4ZoJHUR64(iDp`o>7r!}!%EA5vbq+jb$*1n<8eN6#bbrq8m z=E8-iDexWwXDM(3xx&SO2I%^XuIp`WSL&Q?=Z%XgroisUTT$*{ce zIWmWqLgi)JulmHyI%K~BofuaepEoyh+BbW>cwT5fmlmqi^geCBf$y6lU(ynL-1tTo zJN)#QdVpFDO~`;QopoZ&v?|a*RAg}u7en&d*8y`>eyR}Xe@HD zuo~^Qb*YI?we(czI2tG#Vk9j(dWWX}_10MJ1e;cRAY7!y%57?CS2#PG{q6FbvY%L$ zuIby$qZvObk7hrC8G1~)p8XRf^+dJ)a_=6jnJ~(I61kYZ&czBWifMXp>wr3WfI0w8 zq5QXH*zN#({yOi0^By?wf%6_X?}76kIPZb;9{3;N0ULbyO0 zDL9xY5L%-YS52K?aqDeW$Z=EPZlw6u4!%~Ac=mp#oDQpVVH7G$G_Jzi);5@Lf zcZ%YeNwShSb4)5MT;&Zd4}`r^VI!DEZ&Ofe4cDzeM&PBzg{$sMN-Vihf`5Y{Z?Fz$ zfvr&P4un=x65*nvB1vfow0I@++&;-ohxkd5ZK{j-0-9Szd)A;jR!Dd=+HB^meR80t{ zzbW%to0oe-I2+FIU4>rjQ>sFBkw6$lYw?DaNI;oZEmhZJhXz)vR!68wG$`6&w!+f< z`Bl?yR2K=tyP*)Kd#m)P%DoMv=)1<|&(MJ~p#+u>x?`6#2_7hW)Nsc+kp|3=s5 zZ7g2uo6Y&zwh8GstpzCK6W^OiYywsOnn?T}(~k0HBC&z^qlrW}s1x)o=rT|{w!iv7 zouI{UB@*-(W;LMqfp&pD4cZMF2VK^kNYHtxLC_4?Z3E2(?ZWRl8$jt2XUy2{xlOXI zw%aCTkIq^T9m4UIsI>ZxL;^qD!_SmPv)?r4A&9TZBnf$@JV-hAIXSsE%NO09)h5lp z{HANK&%2tU>F17j;!}dMu=8qMIrjFn+8q2e0>34JU3_LCPMm!+1^mPKsKC~PN|iaq z>ci^DsjHhQ=Q86t1Ncw#*h95~!oByKvNIDEi(y^oIndztv;_y?I z>2>)nzf?ecunotE_j55&~S8u###C%RFs8$(jq&&$D7#>7NIc0>eo3CrtEG!D`uxAV#5s@( z>ztA#&UqB)I*M}~#kq~*oJMgjqd13AoV!?2Yvxv*t0>M<6z3+2a}ve5h~gYXaqgiw z=TMw$D9$ky=N5`{3dOmER^KGgCFJL@iE{^OToPB?t>rFpt{%g#c*<$UO?QfO2gNyq z;#|QY`fI;BR9}6P>tF2@=L(8*1jV_5q8(5=ZO|9z2C{UQ>0Ce8Db5WP=LCv#0h0?u zQn*VUn*PSkT>l!UI2TZy11Qe@OHCgmzbknzAK9~aaQ!Jxaqgcu=TDsLH?$oV57l3` zhU+U%aju^@$4{KwH&oy566f|&^3tYuu5WjWbNj?OeR19{hvq-x66f+U{f(=*et}b* z%O}p^+c1ni&A;N@J*FR@&-II)A{{!bl`g?g-~llPOA3pY@KZx}-8sjkmvMgMU1k73 z^kEmxUv%LIDeMyGuF_mdmx%xPTp~dH9&?yu;|cMrX^_LkiPXpAMY~wR^GD+?U1H8k zJ?oEfu7=-E5+}w!q{sC|JI4+ej=D->*Cm`@!|6;;=Wu#Er>3+dZ3vnn8wbQ+fzPd|xKSy<%rhAeYekMC z2Nv1WW){z!HbZC*PsvV(bShzurFLM$35DVudxZKAEJw$#S^t2nDLkK>{r zi*+A8)sws!AJ&7X>jr#mI5%7R24j=`svIoIi}9$$jF#qc`7De4VlFSnKhd6jTwdYw z3V6Dj@FDvfxxd)j4jBD-i1OpK*ng1Ai}9W2e!8OgkR35^2s@8+T+Bye9Qt1z7xSFJ zf6nplnF!!wYeLZ9;6wRQE&2IR;FNB0DFV3Ix(}4r$~KHzPw{pm`ui@9Z?xbqF?@oQ zYq9?Z$L$vU1BPEBiSwp6f~V^wJ~rtRNt`?EVi!gFNtwm$a}aQv5pe>N+v%pdG5BrzXN=@{L|Xqh81AL zA|W|Qhm-#)lfN|iJ_{lJe}L1|p?IHzz_&8_%aY^1z+VA=G3q^7VIuSc8@k?NxEW5p z({o|!d88tao)d?|uL7z zdI;1G`hO*?k@ovx!J*f})`l+5Slx{K^cSY*_gxL7Z);quBX6Yt>(bZ|V#cV3?0_@lx;Uod>i@h2J1nH2JC z;CA%)8a^(Gs8QUGcyG@&JmPm;K5l70TZQ~prXYLHMQ{hhh1Eg&<_P$y5%4rPg#2u6 zw>iLvYd@2LD}%`;Llq;)F9v>fQfe@;isaL<%2Wn13f(t?ohWeX*WKJdV!n{(*M}mG zq5xjyP`BLgiPWu-=tU1j_0nB|3*ChK8;bAxG!b~D2)=nFD}4M3^E9}dzMH|*SP zM$#7*;9A(17*ePO`H9erD`3?hZ1jd$1^~h?3u1Q~JJZDy>4hZ}c{$$4(J~+99(*;2 zC(`V}J2kLj3THF~JS&<4%j=pv4cONWd+M--OfR0zEU#3|mXK`+YS?*bWu{oQLmlCgNp z#y@(0iTRq4q%^l59s&$`OGz>YZ}LDc(Wvl$8ocRcyWxOX>tcT7?7I>ZXyxNQN3dLcSme+1px@mMHIN zRIuY&9!3YFcgC=;ntYLq`RRzDOpiLM>;=%VCl+tRU3-k(J;eTFT!oGl(FwQc-wKbEIxLXe4CPMiFWWM7%{S zT3U;ixBBs_D1U=gMBiG3f)Yb9UNuTZ4QpCp!K9Ip8G|o%!t^}`6XT2raUpM09Vzfo zuqi?okIpNTVz> zZ4LmYX%Qy#lozDL@?F>LuSZ@3$N_FeEL&z4&<>;`3`)%yriab~|2jweBBSxavV~nE zylf(sMuFAQBzigiPB>c&Las3eX-vnXaN@%=_cnfjypkV zuSMvK{TM;l@d(-vrhE#0@f@iFqbv)3vA-jz*w29oyU0fBMav}GI}!R~zeG^}*q8L1 zkP~zX^l2|f@M1qk&^DeJ<)36k{yEH^cgZ(GU+f17D)xh{{uAY21zn0O=LW_8lAvNg ziexUsN2D+G+klZz__ovF`a}VUOYFOuoSBvg6TKIoRDDGu6lN4gBqv3tTy`E7^da!6 z`fXR40)oDtm!gcIb12Kcvr5$$&liG1fnU^yY4KEon z_0iZPBdC;A2)+tbdZ>jKPi+m@>rs4EgHq~)4oDSEt1wu|RASEexA$6MXGpxg_n!Ot zoIm!Ko!?r&^?R-Fdav28owIPR(O_WWF|uDX%E`lIg{A&wY<=zp8JU%`OeQml<+C(M zrQj8EGEWIwY(ya*K{M1WGl7!cMqa+rC?iLCnFt%M{JlVhccQ7hOeHeUsdQhdniZ5O zwNd|4tVg=k@*`aqB_ZfZ9>+}-nDVF|iA{PZx!y^xCuocZCa9<%)s5Z&w;M3Z$WdM< z!bT~7FHn-DN9b+h_1IJ@D||sECC8YcLT~Y3!cI!(PEiGNr^re|&}URO6IAH^3-qXd z|Ll{^1($Gt{oz^2>n$jUNvU5BY!J-zDwna_E0?Xi%(DBhH&|T$?$ayHFYLQ>@~$T_W8UA z?qz+Htm}vG>xch*KYR*D!khfu0!Ckc-|B}yw;%rTe)!)5U&5RGoCc#Wzm@&)KkSF^ z?uY+)Km1kw@R#<(9}m8SH~D!LjK1Q3s2_f)AO8R9hrgg7{AxFJE$ZBgn4WW7#3I{zd9}9%* zgD$(%c8A}6A>#lm?1G5X-w;MP7+?+_j@J%36u5%06bQO~;acW(c|)!+5l3JJ8FqLP zlRxBcw1H>$xa)k(UgHmj84P)8WKN@$Lru+FINw@kn^ZJu1{2B2Q+kpm(}AE-`H$uR zN~PcxJR@cX;!$2rQK=E*-N4R)OY@I}_o3+H82%}2paB^?29iNq%ujoiauXZMc~2;L zBP-zY?{r&KJxqk!_2h_#pQ%Wwe_|TGwm*t~Koc9yKjIPnhWIo$U3k>~WCWjwq7-HD z{zQDLy?C-U{GO!=g}E9&RC*ps!&l}UDohLAE)9{4S#?}Ue@rn_>^n- z12yuCHT*0M|0WIp5)Hpv!ylyK*K7EfYWM*SU$i|HH){BUHS!S+e~5{_*d`^(_9nQVtOV(b)TPBC_B z@VE04*~e!^C1d{Ke982_c(V;`{a~Y<$Tl*iOv$GqKSGnwMBb{&4?@01lfMl4b(;K8 zyuSEWUCT~XmxF$aWd86Ti|ErK6p~;U$-m1wLAz!1(PeOj3CO-}Nr!@JQ$RE(; zzmNQJP2P&U(RjiCT;xY+@(Yl+YVwuH*J$!fkYA_fyI;BN-L6$kX_qrg4LkmRi*d*0 ztK|*n9Ht#Fr?JyJJgf@s_|93AVduG;anZ>v(vINTxgkgL!Hw zjYZcQFYk~{tm?g$q0WdAGJDvat>1Q69W!0AGiJ)!`BHcH6K&siPr=xKhw=rjY^sbQj7I_j=}(6Hg&s zz4Xa;V@7LxEc)}fAtTb>cJ1r&$WKqjSvpI^OexVdkcri@sbAkGC-&5`J&JC`6g3)K z+etnr7jmDNq9L;7VpGqol@pU)8aXMV$PJMo_lha%CVLJx_53P1@sdL$$JQ%yd7~ir zv?*#Mdo^ro>~1*`y9M#rT#!9>i=3dC4JnYl&lIINmO{SO&ZeFV$%(t|Y|pU654TVI z^-tOxu3OhmG5_nT-^CyQ$=j7kqBkSHD*`OB z!gnG*9Bn?d<6Qjm&G9(q@*FvJYr8!hZI86S9N$0Z*Kyf*G;X-=3^;_W|m@j!}1D? z7Bp6O#?a4(0q&Hr=_y!pV{KmmF&h7$QsL&;LCxP5(nY_Uf7V)rjNMkqdUj8s}XULw<%FXDN{-Ec@OoBX<;*l5KmU zyJ|b^Tm_t*z^NRz?gvJrfKefAJ_n4P!05OseWw9`wV3Gn%tSqBhZnVNr8n0`L zQrGG`3Q7x7SacSfefk*2z&6+?{4GPa#~J!fA$--s*E8T7Wz_%8x8hX(ApBDKRq#C< zaU7fRS$sSE*TVlO#Bz)ong_J8Xa@YJ!T$-wA~jkL7O<3NwmO5yfLMUHoWn|I8H^o< zZumdEwViwxT=#5T_$_$hy*T+3zV=@cyU3Rj_0{?kb1~8W{6?0KD7B*QHkP8)-SG9J z?S!xW>EFc##`|}E6xZTxkk9nOn(9fuwOAJbYx0?lb!_E`eX2hdYxea6?c_5ue_P!C z*018y%kRZ^!T%`K-vP{@`|5%Ah3Ll}5ByvE4>Py4?}$9o-ZuW*_;HM<{j)xdR{-}W z;4%_z@Z76Ei$9FM-tn7%kB@!!<+%IyUGd>5Y|+gjxy$g&PvaAKn>>Ry`i!w=2ij>a zUdCQ)Ux=~M{OnutogaM^pMUKe@z~&#@xP(HQqXoLwB7xdSp0jeZcARa&-^CF&l^?W z#CHz=Dn91Q&*KjUBJD-&lIB}HhFw->h`zm|^5ENRQd-{j5;e8Fy{7VDTi5JEZ1t-x z*q}Fyw{0+cqljfSlYBm_*+M$q!`GPqIlY0+1uRq1HFMpg(Cxlv^*^VNwh=KNB^|_R zYrNL0#zFB~4^b>^_3JI{{M~(mM(|yQ^ZDj|`X(XkKGC zsq^3Er`8bCpOMxdq|hYQBEG_Mzkxf z1wR_Xs(y`rAm+mx(avHnoQ`%BbKyC(k(l$kF&7!vuWyf;Qlke(m2HQQ#tYJjlgXz+hZHpb{*M0hB+PQoaQo&_gN|TUVBegMrHQ{ z-&Wd~p~H^3D^uP^AtJaO8B`PjK;eY*+k<}-J+#)<#?_*?Om{4?8~d{i&HEcicG_6^=>o*{FzU1sd&3ykC0dJd-a(orW--=b z|FXM^F}9mkPl&FB+&5*ZQ5oyh6G$uk$b@JZ>4YgOs^sOOXaISf^G@Lmv>WU3Zfo_# zXhCG6THZP_8bSF|>hk^+$eV)N1n#=F;m{oHl(j5--qNETMcYVs!2}Nq{)e+nX~wV+tFps3^@$wr zCll0Ql3#o~SiJow`#s0%SXbMCYO|DEZDtv6GsQ=^9IKPN+B&~A$sQW|vJ%rN%*zX>=#{AoXJ~QEUbog5HtlTx`qXXUMFY)f~?!F@CmlHTA?mh^6 z@-UD26lnGFX}nKNi+)6T@W0UVrRmYb$RD#|Jg|4$t8-6cEv(z3w)gbtA5cd9@qN%3 z`b`ea>(qZRS6`kdC$<{d>0L%v#nL0~XE3*Kyn(e;1J4k~*>3Vq9MzrAz2lDBm}r{ry&&wy?<4(TvpyvJax6?x%{<;sblX!DYs*rd&m!aAbOr^t!* zA|H?wKSDaOd^XM?b9Q__ZZ6Ir^Kkws$N2;O=8y7=U903BXO>!Urn2rh=P28OvsRSO zyR>&z_w11+RyO<%8|-|Baa`R6yn(}UoQsq(S3N$;s^!EZD5vudt#4nTO~kn;`OMRe zbI^hF-B&2{JnSeoB4ZRA_i!7fptOc`!Vax9XNh86`kE-#k2ol;W2Y#OHILT2zadR) zYzL@z?xnPCl2ZW*u!9%g3N!aB^P44*gRY?+o_{ zWFuL=fU!45oL_s$l67x$!FkJu{g6KnUF)fNi5f-M*rTh)`aAWkJdh>Mz<0H^Y(|8! z=|^!^p*E_vnr56CC?~pB9q75=p!y!~13QYpYn8^HvIfdK>=^$tn{~JvYd+axz>e%u z{KrujL-~X4#|+qqtTwdB=);F_9`0K8k`fo4-yi0@a zWRLjN&jnA|r}|NUqkA8U=M3Wc8nG$u5S9~|#}du(NtF6M**3^(eG^8Oa{}_8Kpt&! z%8L2HigOah7Ritk520Mtn{3C7Y!B6+?kfHco#cAc+K4$1dr9nraK}US&BJ_)-{6{G zMjdHASL>Lg)P-UZ^&wrluh8h-aw~t`0@q#S0(`sYO?+q+R)SYP8s^9k9tyCw) zM{PlN_oEFM##$T38pU*!*G}Ae} z7_r50pHyR%caW^zXzNtQ&bs_g(zzG?_5jNE^D)(iF_nY8!WrC0eGMDc_sgARrwDdv zowygV?X6+aHoSXL4{9SrinSAKKKcmj*oN{L>N$w~!bX-thxCjIn&)Xv=vwuv(oQqe z*}udb2IWj}^BFp;T)DI3AL#E4_cIga#EoaW_Yx0x$vMe!5=O1xbXQRfpYSA)@ z$*wkMW&6!v&Pw+VK|dC4NA#EpcV@8v0_=aew%m(09=j9!Fgh#5fGPG8w`2ct`_y_4eGXK)h+oiue~L|`ss?yfAE`aXbbGAtL{N8e>T>b+4^_dbFc=qjXoJ4GK3|Dj$K1) zW@XQ#?tkcvPj%}a4w)&aX|Y61IA8)G&cGJ)J=xmEKt#gt!Zoz-8)%u&MlyO3Y=?4 z;0|OR?lzvneT5PC5Og0x_Z)O@L3a|*Jl)dX^tT=HQI8&tKePQr{4XCr(f;iFzlo=^ zHoUPo*7ltR8L2JStkf1qc4|vgN@`0%xa^>HcFIA=>iq|s@}BzPnANYeELc6erGO1+ zw&Ffv9BcbvVRqXG-mJC{n%KaDO}L*ZV3!>HIg?+QnDO@J#aVA}{^Uo4TAEmPvvvG| z4F$Iz*zhXOY5UK2S4qrx5NkW$l!M9|ZpzRJj>!ekz@6K znpfs!IkWQ%@LKUY@Uoj~%>{U^cpZ4N+wA5Nyw<{i*-?0>c}ZcN`Alxm%nC#1Fw%d| z4Vw!J8_d?WYV%#P-TcadYO|xf&P;I?bm44QNV*Qlqd&0PEy=M_d`x*du`Qt|MOFyQ z;xxwjr#W5BDdXj*X60wF)AymD(z;LUAUnYpok0CvCs#eXm8Et(g>f<$_kL}R?P|hY z`wjZYR%Ys02fnR-RB02=qvyf93jJa(V>9Nauo?3(pM1VLy=9ic)KP=6x1heR^b5>K z9i{<0*Mdj7bL+>IE;pn%lkR7*LHolP{LY3h**^f81@+ESn}OYlUjsIipMO9m8G}f= z+CEZ>ylX_3W!I=2>#nQv%66^B{aF{@+dnfd`aF$g&oZ=KHx&6Rj}%m1gIC)4Ql(|% zZ!2p+O`uI#t(DKFZ1cUE^Hk+$DyfcCK53#0aV(evN?ETC*pt;(>c2ye;yT%v{(OHJf4hPnxtyXrrES*)I#hiBz zZ=*4%dDvw3?muATnrV-;fAi`uagX1IHx~ah<~}R>)ayq1#j$IK4_~k*d+E)7d= zr*Gc8b&*OBL;oDs)+3NRjCuYr=l~O!hkgv}q!H_61Z8yQqPddR+Y5Umb_Wc@ybE$yL+)xtu0S&Qtg_MPs5YVh{WsaH<~Ajh@8MK-=eB>K zFE{P3EI|LZHoZ{kz<8nd-m^^J^{d+~yEd97-(GWp@347<@0@w0Z_hFYzG2sle91Q> zf23~#ev4U-ez(?;(_uQ-UFD$hakjf^#&3R6DVeS^t9w-StU~jJn5RXX;Oh>%&nZow zG@)ii*d;A5K)V=~wo%V7+Hs*B5A!ow_O%VN_xl6nuCuWF`ll-o-HovjOT(Oxy+bzc zmj>g#3~wjOzrYwM8!UV046}G0rPb!9QoH%t(i(GJhQs_o-t7FJ=FQ3ftZ{Ij=Py`5OUhu(Z^{`A{# z&A+MU*8H=#8^L2;)Mx)}%;^}3`GR7k801vhYemdA+-CLOGz@n!##{4`rQMou9C&NK zoXSwfy2yvny$}1FccZ60Z zC#Zk#M=Ub#3AFQjfA<#j9gg;?^@l=UdvN(o311#b{#pqi5cz{Eq{7iF8zdS6p)r{x ze?e|^NOFfHpWipZ>B5IpUbhb)RJk2ep`+gASRQKdN~1lV2C3jWg{s<&iO;Eo?+Fv` zas_MrAs1&%s8D!cGdd1V=TkiAU?UVe0u-osG zt{xq_8bXWG~xgGAybwm71O4GE~IZLK@qOU90| zW=eQQdz@0C%kHR`&@5qUg?HVby#ER~0GYV)T?L zN(f8mSI%2FM0T;giRh=X6_Te+c1_up6QezIk7vnEPu9e4W z4`Z9f=%dN7%L{#1kUF~^)l$lryMjKKhq3v-u&d4$lq&3DWk#4bWkQWRELj#Vn~CZ> zo#UkjbkKSbr~9t)63YBa{VG5v!rNSo^DvG1#rVE5++zmyZNQhFGo@P9NUfq>8@dzU zNGi|bAU@udDjU2tt{?>+a=Y%ru&R~j1npse2r&CxAt~&aCRHf%6%K3@;BZc3*yVGg z+2&sWfb=WpEUJ_&OXpY2o+*_~Dpe|hfRwt8lWHkS#ZDQ=N#M5d=IPT>QG7Wn)w*Fg zI9{3}Rk*D*xC0F#R$--DGB(HO2(D0OJjDzvtGDCdq~h_?q+%)PZ}2(Mb4o~4Eni%} zBII_2*whm2Z2Zm@lG7bT(@voIx-MwH%N*cnXjVtPJ>-I|X(BRT zXrkRgJx~pkJ$Cr2lLLzefB*cB)9n>Ci)X#{COzU@%7s91r=UDik*v@{PK!A$n-NAx=;OVOMsV?d$#XBuAeP8 za(;x<*RE8HddDBOY|`rQzKiP@%G)`AAE&$f(2q0~Sz=oKwxwLZST+ndUGofMYzW?y zIEaSvqj=$8y1QROz8>@)(Cwg|phrNn5LnKq-Q5!CLeO&1deCanO`uJnyFoXBmVegW z{SoMcptN1w4w?g6c(S|uCeU(F59l7y=RscsZ38_5`X2G|BVZ2y7&;sN=`pNa!crQ| zhT&QHjCi?tohQ1xucbmux_K1Hc)aVuqcUchU>TJCq;ZSEdi3K#$Po_@=~m(04!L4t z;2|C@7QYwvkS2Sp@%|OOQc%U7^km8wlXZi!{OHGrxJ?>wgyVV0Wsn?yBv8VV0XxNy zWl(Qz6p`N;-Ym!{F(WPG{TgpJcPUADov(QrLte{YcD`&P(AAS00s9vN={yobS~ zy7&0rkhU?^`nc(_l(Hv`6@RoIejmQc2z*l!_e~td?jG294K@fHdZ>+F1MfZXJjAdJ zDtQAN{z{4;!-iHqnK1|W}sg%^e@7U9-{bds~3;_6oYpjeyIP! zPxX^oTLxMm%XmE9x*>DpfU+mjHl>#RIAycRxH$9Z$A9|Zk4N4=tki_QX}=Y30Q#HI zU$e-zWzc%Vlb3B7T(%*5sx7Bv}6!%ZLj$@0;M7I>U*f$Qrnh@rxgwrH+e%a4yzi|rnu^`tt26S#s8cv8b5 zKZ-xNQQX0c;th`G{po`bKT4ZV6;mpRiZ19CspWqYDQeDA@KPI6GtvJiMh^}7bS0(# z(Gk!Q&=JrP&=JrP&=JrP&=JrP&=JrP&=JrP&=JrP&=JrP&=JrP&=JrP&=JrP&=JrP z&=JrP&=JrP&=JrP&=JrP&=JrP&=JrP&=JrP&=JrP&=JrP&=JrP&=JrP&=JrP&=JrP z&=JrP&=JrP&=JrP&=JrP&=JrP&=JrP&=JrP&=JrP&=JrP&=JrP&=JrP&=JrP&=JrP z&=JrP&=JrP&=JrP&=JrP&=JrP&=JrP&=JrP&=JrP&=JrP&=JrP&=JrP&=JrP&=JrP z&=JrP&=JrP&=JrP&=JrP&=JrP&=JrP&=JrP&=JrP&=JrP&=JrP&=JrP&=JrP&=JrP z&=JrP&=JrP&=JrP&=JrP&=JrP&=JrP&=JrP&=JrP&=JrP&=JrP&=JrP&=JrP&=JrP z&=JrP&=JrP&=JrP&=JrP&=JrP&=JrP&=JrP&=JrP&=JrP&=JrP&=JrP&=JrP&=JrP z&=JrP&=JrP&=JrP&=JrP&=JrP&=JrP_@9h`^GaC>Z#<_49_%?zMc%*@S)5-#Kvq(h z@OKXcKeC7+b9!Zsmfo*%2Q<}mB#Z&jA}lU z=g0H-2J(C^r0gocB}CjdA^hD=kRZCe{BUzVOZ&+C>71Fx%x;|-3fp~7d(g>x8F%~KVHj4-@w8%!r+dl5 zLfGw?5OjI`4mFr!DT$}pN#ZH?l6Z{ze>Q;*z~@;;87A9Hm^@4-4qMG>Uox9`+4J8a zBj%V!-uL>lBj%miV^yDh$+U4Fo-6y85q1Ra%Z?}$v@aPEYtwhcEMi^EW44S0Sauis zhppCS4>tr|HZ-xV#vKk3ta!z2flMEd2@(OLd{K!BdK>%Xn|kI~t&)kQup@oSr}oUh zz2(zbU-OoUr86<#^_I_IV&3e{AHc-?+MA!rL?7tQAIQW!)SHjhX`rlSDQOeCgo!rO z7Mj=~CfcAk|5BDaMb@&EvuY{CvCG+O ze7))|Ka`1etT+D(ww+nEEG2DXxvZ^EzF4(;tC*Nrt&{n~S+Sun`(ib}AW!`_g`Muc zz)@2)b{Na*nJZK-nve6te9Q{Ix(;Wv_mc2e`4{6)ksH_rzQhfw@)x(eh$qsI-;+sx zN5S5A)l1B8+$|f2gc0=>e6i0^{a)P9RsO~Ow)8vwW->AV2*1LPm~#Yw)_2;;Vi7)9 z3i-K6Q2c^l%nKHC{%$Ux%lS8R{%f4Sp7R}?Up-t_IO<=X4qJ?9_ce3{!8V6&CIDa3P7kt%DE{i4E`B;&^z!&l-IsasmoxgMbDbAOWr{^2I2FyR%eBKat zQh*=vbCdW3IA2QQ59R#AB)%GdHfx+Jle?a$0P^28PQ~6pD9|~A`cE>yi1SbNIX?_y zH;t8r358m?yf{~h!W+S-`o6>QK~{Ox_%B`;9E$wK{2y}r;`}D;t9CA4cOF;lT+A1^ zWhddL@-M_G?C*u0@4_=C$wE}2svV5S5xx0p`(wN(@q6S)^(`Mr;tTn})qTrHlK4XY zRAJxp8Do0$+xn@mYA2g*EbT4-I^+%5CnWP3p1yGVqauItx*r3dv5WU1XB7F1_X8@P z7w?M(V4k>)OeHwi4T2{O_@p^s4?I-VBURC&4T;%sHZYNvx2X5zm;u|ks zmp|;s{%67tpC1a9s2DrN`Jz8Nx&8CtQ@vhGnpXz_x4zm@>W4qEAO3vsv-`mR2Jj^w zo0ho1A8fpMehcN&16Tj-%jED&0Jp08y|80bN$%=6_;-vAuE!! zUeeRX`E3QVlHa8KW$cZ9?0gJ9^|$2yaJ(P+IPr}bFMR$N`bm z%f8%J>m!~IJQBuL#~uC}o68ptu3)&-u{rAPA(xU0`Ga945vX4gayvq5ad1T->{pY1 zZy@Ljsr-=J>r{n>pvUFY7P=bUVQt=K^SEmqggXGIUeQxHc?z@5ud1+JK(DHz zjIgV+*LYkAuqK37QS@vKpgNY@>RjP4x18}Mz5rV=yb5@Z|GGEfDm%^;&}XZJx(XhU;%%~-DrcD-Q%viP)7uaRz?VF zMUbiU`h7jMLRrXftH41 zv^nrMe8CEuw6@mM5USTY3#|w_{JvT(%VtAUiyldRK#ZSWJwQ0Q!{58tpwejFUTv-n z1y+l?`@+3Exf(B29L>avLMyysdkttfsM30oQm%^Yd<{i44Q`Kfg4@Z8V9gZ_up;7O z?y6@+&J{jbP-!@*mfYnEhTMK%uZ#_40MkwiJQeVS$)X#7!>&gBrIC&@ztTmET=jfx z*E^vra;!)R%VM2Bfo5DqUmfzQzDUE3-RpKB20tw0Pm$tIHI0G9iqN~gn2K2uzdu(D zgWpr*B6~O-bk{Us)?`JNB^BDdcBT5iOAm({{g#96mhuFQGN%e*qg2xK7t2zwrwBO% zUZF4USvPV9+1EZKt9_;+a~8k-2>sngSt;Jl8DgI({EGdpB7{FAO?PEN|A>(|4EMJR z!F{kwQ$fhL&=>oba`3Qq>CqRzua%Pm9`Spcnp16vG||O)>9;K5U)%$S(^RP`D#4jY zeT1B#w}DN+jfuRtFW%1e$v!<44#|l4iAN9J>yT}sFYc8)xqc%F;Sv6YzNr6Qyv(Xd zv7Y5KdE(y^+bz)7ej$NwnT7x;^N^u=8NHSt3zH;^3v z<|KV_-=51m&NIR>r^)`eC+UlO`8KY9ggUq85$z{fFCd{E|Kk3BQ>iK;2|-Si{kMTn zTXf+@-18slQWcW@iM*h1LhMqFzW6;L_D|{`fiCqu;a~i2aFXlCc>m$Uxu^e%b~pxI z5)u01_l4K^O}&IN?IYqA^ivdS^~LWHC%Jx`RvBrbC+Js6`r`MB?TS@~lUeUa=m{>} zGZIh4FV0hKTt6TLITiI5?0+CZ@fQmQr$V3RckM&6^ayG~qA&d~xx^SDsH*ko3HD`3 z@To?T5a-Zbqbf052y*@8{#S@{)IGVH>{bd4&nT`ZNRl}v9x-->4ABVJS2`#1X1rd! o0)NpiG`t1&NcKDtsu-TbNdYJO9}Jnk^l!LAl`KtCNXoMR0uyBuIRF3v literal 0 HcmV?d00001 diff --git a/dist/test_kernels b/dist/test_kernels new file mode 100755 index 0000000000000000000000000000000000000000..0abdfd4090519e5e6f9446bc9c7ff8664806c974 GIT binary patch literal 26096 zcmeHw4Rlo1z3-WPK!lhxAZQekqaE7BN}2(r1fgca1kT_DqeMjkXNF{kM3T%nA1P_Q z7!CNChE#4#E4Ryg%j;|VY+rASF0V_dOJ)d!fImcdhpp8_)Ug|L6Yi?~jw&zFcgdm0~bRBBV;+lL%|LL#GsSM(-jX0VtHFN=CdD zX`D0!cp9VW{6ZaotA~cv{xE~niA;LOIX@ZR%jE63Z>g&vVfubu0)b3Flj zc|iy$%9F~3o(Q*>3{Bus5Rhak5Ts|_9Kvpfg^Wsdu7~9G4syL@sBn1!MS4R)r}$vy zBIJ_cyc59+qW!P>S<3YqWPV&&~cK;8^kS&l+tRPD3pJ|(IDJv@DA!3$GAufLFSVehOJlN)pAf2X=?@s!Em zsjiq%T~!-eHeuP+DHEnl&I@?+CQ4-XE<9u}?zh3nXdp?sV^@nVNa#(W0+f+AQ6=fO zgU~4*WdM1qo`LMhP#B0lZxH#@gV2u+LjQaa`u;)aV?kH&B!`Sa{5%iAf&BmFAoRLH z=#LFTKQ##bADtv2J)wB5PB=<4db6bM#$L%}$;CKPlpb5wiFmr6CR zYCxDoS1IR?)VN^5Rb92D)?MKc{PJZkho`C*u z=zPj!%=QHO@IGM<>XO_DgF98Xx1qsmEv@<5^{dj0C!vLx^EaptJP{ z5fQST3`k^$@+N}f)MG@aM}!=w9wVA4Q=;TM>8;D60^(t+Q;(ATusBtc3KQs=33N4q zj*5iM6YZTQgdiwb;WuL5_|8 z7;kx8A!3AA60fT>hA@5$;WQO?wKIGZ;WSluMHv1L;i-ggVEFF|rzx^)J;Q%RI8BXR zYZ(3u!f8tEYGC+J2&ZYF%g68^5l&NJmy_YI5>8WJS1H4PKsZf#U4;zy6F!{qe1^Y3 zI8A+B3d0u@PE%f&k>Sq~PE%c%#PBBxmkIB^0mAq>gwxd4)d_sw<}yT4f~6H|)4k(! z2njDVscQ2=lcF|TP5ElG-IU`|<=S>ar&4SOy&@jDdZ(OgHyH`*0M;VsrkS+okjN$~ zPqlf5XA-LtRi3*=4X0~!B}r>SqTmF<(O#x_)5GnctIg>jU@V=H4NCzD~cJ=%aM{pNO8V(~p3pT>}&Hza_e)(|?Eg zSv#)Ne?xTbI&;yC`CR*0r~iWJ?K-_4m7~3@(|=6#_jUSNRI2u-PJfx`>vj6?vJqdQ z(|tr=1-eJ}Qpw5QG?N-W{2{d`DoaJqMUiV$Op;{21V5jJJlb@jytY+4i$(=v_^@`D zU_X<<-V64JtlS}S@f+O0ERtxvlEEze2_i6BIuTB!!ilEQzr!(&Inwm-KLdi3Pww4K zob=a0Ss9Uk7EznZOlpPt?NO|8gluf<~2|IBA>;XTi*;iGEf@!rz1 zVspg21IGRi4SeUw1KmEdS}-5ZY~j;?!KtWzG>Kyob{ z!p=q}|L{8ma&M^NUN!pFqiXa@s%mIe_uL4Mf`MJ!Kt`<7BagH8>63pAGyVa(QfSgs zYU7OGLlRB=MDTXCX~ta8W24uBq?NS5DIg_Po42>%3K7J~u-RzhufrD@Z?% z#t-uIMxxKk_E==^9<@1WD)Y!M?#)2ooX-2^2d1}R!5pK8qChi>ryD};8F=l_2$z{k z)#mz4kGy0H8EIF;JDCZO?2ABT$tlRU135$Td-3v~Kx7}W_Q8Si)C<8TPRD6HpX1r; zvC>OkvJJ?`gnX>0N+vDx5{+IvKsv__`*uLxc^Q047Us3Mc zAh#ECE-1LjL*^5GFR*_M%?7sj}=8FF;X`aK2IYX-c^I)PBWEb*ATHT+3 zi`+89-CLN^LcFaeRerC?Aip;!MQz-YrZ%>wMUNPxof)^afp(3(W6>kox3z_5T~iA_ z4efzls|l)JjON>|rcpGn04qb*XaFjea~GPN)~1D~3M;yKCGxxckes{1)BrhaligI^ zw8FHS=`7zb=a!j#5Dhj3P0Kv;^6hf&YSS7fxWYueyJ@v)ErGSBbpT;2jc?2CrVMM- zuT86nVWs#>rkAWu@0#8OK^`Yr3}-Bc4vXQO)vQ^~pIOa&qY+~?YBX%O7}}y+jOOcB zbBo2WFS;!=+LCECpVSqit%#EXG^3G>Xfy+wOou6;S7k)E4UM)8wZe_e=$4FP!&Qsn zP_cP8Y+1vXt;mBB)}q^u(bl2G&1v7W7|s`)w_D8zEz#|n@Bs=|^IxLdGc1NHR&x&o z!L}G87Q;2GxfMb~t>$BJ30A?jgl|})tw`kpqYj0C#fGmfhAx=3n)g`EJCP?CM~>jG zrFl9^7mk^4z*8fN1gW&c5yM}s=4}>p47nvYpbslhh6&^WxrA7s%uw>8vD!4sAg_EE zb6I!=s}2}8TEo?-*RWp9*05gwnAd|=yeh(|Kr2Ry%CKHz;SWp=s17TfdC3|^ty^J= zY7|y{VZAzO8psBjT7x-4?J7<2$TzHMPI*Kcnq$O>X{IoAQbmiQ5ehR#nm%KURA1Rj zquzo=mW7r@mgg;w7VS5eXavPL-6S{NX^s@^lz$Q_*wc8B#{GWdt>tN3_y#6{eX}HK zUidsM^|Yy%W3h@ReaywAF}X4Q5V!_=_|nRb;9LlIk0&3sQY1oSF(^NSzyizjUCkoh zl^vnuPyp+>947n?a2nIsk)iILed9d@ey>|gVFKT00>8r~L*;{+Wfp*=vkMYBeB%S=X0Zm)o)Hm@-C!V3DnF$0i{ zWV2-*jRyIgqGq%3O01WsC14v zIqY4^-j$RaqgH_R738sQ!W$)wQu*UICFyPfXQBK?tM)3EFU%l@%ASb%5*TfDnG9M@J`#b(d(hD4>yq$= zY#)U!k}6AA5c9}+ZT@V~(He-B$+>xLh4N3y!fLW`CoE9Mp>eZ+sZ-&$Xm#gl#Yp|d z@)TQJ^Tos!toP6r(+}h%6}kq7N}cXDJag3MH%(5Wpx~V-QMTOygGT!GfMnvo!Gs$A zOl^Lbk%&djQAL+0Rc#VFT@wZsxV#D-!V*4X37-yYh4OnBEzvVZOZ0TM{N6ruq&VDZ zk>A^;9m9$!tU>x6NH<(NE3b`M8m@mSuibKgdtt+MBO>NY_qSVKxpo}M(DL;DRxJxk zLi|;BXCx>p9%O zVT8kW2C=F6H|zzd^c3@6wRtsjL6nqsH##2Q@Su2aVtGANszIw zVM{tM7AK6l(3|97n^H&&R3$x00tji7}QFsW=!+3QcBPO{|X(F>!@{qVh z5=(ZG8JOtCI+7R=qXPzPmvV^lYxK5y0?fy<~rHlOM>Jy zEh*u-?8PDrkQqYxeNE&&lnkZt;-_q3`&UvwL`2|nr3MBT#xPDMDQhSoNWfAhqnM|GeY*x!=CdBiEAWv1r|uJn+gxgFuD(~s*kL4?sT-A zDq3B(fZ{@0-b@pu^kwXYY?8eKyc<>QZH4H0Bz=-nMCE(21_yyGF#}M$Y-Gx$=jx^? ztwsHpjFkGPWaljDegt-RyRHa9Z13wTFoH2Q5=2v~3z~d7zgVt|Ay$W}kDG>7tgrF@ zbR6CtM^Vyfa$Z34pEBb=XT%GuV&*t3!dbY(9AzZ+;B}z%sKC5vny@9i{-#?( znNeoCFQQ;s2~i2}hnhZA z7IT$e$Vc5?P;dX5Q>mv^AE#FbtqM9x2?+Iz>K*io;0OmJRW;_WE7)^XpJ7We*^3<@ zYmDC0w8}%)_I(cVeWy@j-V1W;C#%5!9LT;)#Otugt%p_-4?9iNHp(5njtvS`)RSEM ziFcauNDr21`?wE1vSUwN7yNxE$YxKm-1-$WN0sY31o^6qWTBg4w8?XQp)**g#z@2p z`Wj098X7G#_zIJe6ZinX!lZD5KK8*zJ|}DdWHt&pAtDH+oX`%a<>0h$ zQv*DrDX5bmPJ9}>V=;DQ^sD>m;HFUiS*x}gi#0ZlCTCGQfi{z%&8%sK!o?!VYJZai`P-kRC-H>=@ow)eY7nD6+rt$nS!xo1}-vt|F8 zro=^O!ur#0*n@4b0jC^l^wfQ-d8bYO(;hY4u12*ns<|ClL=A5NZ=AV()vTCWa5eO? z8a_p5Zse!-G5QxfRNF+TExb*85P3!t zs3>O~zMEDtSU1kWT8$Xj)`0@Xub@>Z|GZT*gOO7eGS+zyBKq4~}M0$aI2$c9<&X&NSjt z@Z{qu#DiV`OeY@fMH(Am2w;s4*6Uz{4k9|Bwhxg`9rOZd9*jlQJ$9p^AT6H$w*Onn zUB>Y_Yuj&Qsjn7nVwqxx%b58VX1)>nw4A>KS;TT4BPwQ?vUEGOf#w-YAW2JmTEuJfAqqF15N683 zT)_S4dIZdb_G&3($v!ntn}{(!?!tIN+K~@5E-f2~`I5ch7!rDrTng(s?(Jn=yr3KO zhd^h2yn~1_jsE`7d}#saAaQKADIZ;b20tNvw3(^chAsLoHENu2*_$h4_HYIEW^`s! zvhvK!_r8s7ds%3SKkN{cpG3F+gU4TDR&D*tJVI^$7je}6*5^t+_cm;8LU zHEpOtHMDcyo_=WE(e8-aJUc~geuL6QnR?bVOf}dV_aWmEgIdrY>eZH>KoN(JU|=vu zJ|rv|sg5ZQBb@fbd5)xL`Ow^A&6Kn`u>`6|!~$b_so4XuOif8tzUmLK)M> z3Xw-hgm;a%m;}0oWntIpcwVS4M&@brM~l_@D>a3f_DExTjRDgdqH>T3#rt4TNwgBGoW3sL%fSWIMO-9Q$U+R8G zZ<^}Ln7rZ;vSEAw=Bd-RNz9wJKaV|B{rEp#C?2nV*z@7PegOr^F8`^u#Lxz?yTq^u z;6RDtuK-_^7+UOX#vE%iN9^H)YP55N%^bCdBWkp5oXxz&9&W?;mZ9q9cH8|Bjo2<9 zwHew$wO#JE!4&n6Hm*yNs2V=2HtxWQTm<{&%V*JhqPARS=AYEdS8awCTXe_p68YC{ zcd4Gx*J(X@=-jW;MszS zW_ctb|Bywv@~95e3U!#tVJD}b<*<}PiSv6Og?S1uvp1lwFEt3LHpWuqSAPVQ&Pr5# zyb51c8;^9V1#N*5=jE5vhTvr#g2Ua2(tQfw>jJm6rVWwQSfu-0>kvchEQ7(eJwjs* zwF_EpK}0E%-y1b##oo|JTR3VfxEdH%)S8xMC^1C4&)T4gV!^+W*dJVP_sC~DJ@O%(estla;j_QdejKo=)KarSlr;km4oU<6J*RtFJLA(Y9ha=UTD+ z?Db-K{^eqM-qvFI$=!(ef_8*8fcdBy2*S@E_|`+8u&{OOfzLfSS?Sq&zF1yB8>=2{ zPJM0oCUgU$zP7jZh~8Eb#@nanS;{bYKW%x&7QSW+W0%Tmbt5go&s$KH=^Pk!Owj@6 z_DCV>+~0w3?ne35_fR(_;gcodFKx~DPPRACyI?+3QgA?SYy!j5G^8r>HWuws&1h+| zuSC_Kv*#YQgv(Y+w&uy1_JT`+=k54To7=0-pk;NSv}uUtks%=)YMqAQ$C}f!)PjQn zdFGEYvx=r=W`)wJ<~8WpbFZj|ZEEATR7_~HAB|$0t?}v|EK*|FT@t?B6-0^IVy#`( zfNJ<)Tm$p8)}5v{&oWpX_{xA`bitw)wee&`UU7}q$18WpE6yW=h3oi>dW(7xZN`uD z0(><@dpT2Lj+B^jwBg-~hfYrppx|p+d!XilQxClM(yPOwZ_HV3^e~C@XdL*Ehjwy+ zooZkOeuHXsC&{CY_|IlC<|w?o@qLl%F+r2goBRjZIm6AGLmBgKuIbxRfF?lSHGLgF zy~)sjk_$+fyHH?>?UYxZ#nC$ayFsU0H>l&>pgwb>V@`9LjBnarI$-JPMC@gSP$!*i z_8cMHN#DG^P9|`a5p;=sk~qc)EJ_DuJ>3{TY>m-ETftfR)hFo>IAS~Gm9zB3ecR{@ z@eMjp@NIirkK%K>mN1^Wm-$?^70RTE6BgG8-O5sRpkKuE`!d^UAk<`L1l->cno~4)ZFvVzQ!E^xGgLpklFNV%PfY=jziNPly#Ey|j z-|GmHd-ejRk}oK74UsQNU@Nfsnys;=&=%WbD>x^=N`HMsuDz(xyFQ1kTJlK=bK(5+ z6gh#&HHw^pSC~32@(3n9f=Q3CPg(6U^GH#;Pd7N4tpP66aV{;8!3L5E!=C{qm^K*bB8U~W)EY5b545_ ziXeqFTNtz2j*^u-g7-EJoz*nzb!+&pqHx;l<{kF%PVL(e1H1b$=2r6#E7(}%U|}`X zj%`CMsurB}pF#sgLky-g4Q;W7KYF~mYz($Vdjom4<|5iA1Tvd4|4Y#$nQzpk+QM6k zn$mt*)Rg%KX6W0C!kNEBQg4)mKk076nhB%arQk?4e1%qEO_-*4gs!N$SEvKf6w3cr z4%=d&x5BLwxK#qTO5j!r+$w=vC2*?*{+}!XgOnmE$^yl^R9U2q4-6w@mc?e*5np9} zpsGBe;P;jRx4+IEP-?t2?%H6eM&hh#%3^Ot6&S8y(C;e8k1M%KovS(oGJc|2uiyuh zRrvX3xwppWcLxHBzKIAZerzSYHA;mm=;~uu1(XUmej!=|r9iN%TnbdxR44;dbMkYE z@dLNtD|vj58dtFLWHWfc`)p=fz;;!)p>lK&Z zj}jEQwABXPOWb~?1X-OX&8qggf>pIklv4cIHOPsjl(sTARI?b~O_Pe31>rBO+iC+r zf2cfIg}R(3&E$0f3^KLl{(2vArb(rG{Y{hR!712A1n+A1gfjNKS{TYHuXF|6xzi*D z(Kc#))u=Eu0mT!lu2#zH%ki6D#!M^|T{^&wDV1KI;`Mmw_p;+Fl=4b<`O?5c3ZnEw z-Vmy47@g1~Bs51P?5!;i`Tb~sNY(GoBRPRk zjS{HzhN>%+#cpMJKIwBUMe(~m%2Ibd%E;}iS&*NXH))YFLm8hxaar;!=gC%O11h{;uMlEPM)-yQ*D&WoD%d|4qzS&^D3FfD+fOgQnLn%`JX>?#)tM z8$RI&)3xPpC1;|tn0`MT$YrG>QSstG7)GB9;b+M8C`4E75~RR-Qy<^&4b@gm2>Pph z0VQYt<4;mks;aH?R)>OBb#7s@hk2S2X z`vxlU>*uBVkP&pDrgVSeMThcoT zvf@&F7%Dpdo)~l|dN*_U=+6#(QyZQ);(NN= zJue%^q&`A_^iIDCUcVz2YXgkzjKvnFBC;zM+W;W#j>YJ+xdNDhj1_(oi_Hgg0#;K- z_r+rC2=0%?-Uaj>h{g5+8V|=}qi{;@>x{)p0XH0p#a;r;KMHxm0k;CK!H1$=z;=AV zx|h;L!1pB<3_o~MGAuJ1#$*o7SOc3xzZZ}CaV&;ElEoi=4>b}n;hBu5y)70~$d@eR ztgP%g^6>9xG)Rx$^T>mfOb>v{^cMnVAaD363H>Z%SyqN*M>L4X2>KX&X_bnyjIXD_ z#E(;pvXthuqO6?N>6Wbg*M?ZKrZx_poz-hdy_7aAYbt1#tQ^P`Wo1LBC@aGUVi!o3 zQcU@B7@4R1Qy#4-{}(8)v$HlBQWqH*6X9wC`}zac9dC2m29>v^QV%XtXxBm|sf8&gs$iAi?SsBH)yd{|BO@fzi! zF)iiktn9~FY3Z{0kUb9B&&iH8tI*z@^2F*?2P=~2klWYNQkJw7w=%m6A-fK3ax%|b zqTSSQz(e?fdSvCe3bKDdyZr^1{V|neGqs7;RF1Dvn`j)88bIMwQN$ZTG3ro(D{atU zxIGpt1l03q%W|eZlrk*KXyr`>3`sirMQVDgbOrZA==Tv&d!~?&=i0Uy3)H4)%8ldj z3BL)Hrla3dVEq9&C50PlGUy%(a$8C@;^K+y6m)W!NX-}y0!a3UYnE# zDMB=+!=QoT_#n&zto7rQEcr|Fq)&lp^xO)!O5p#85}>v%?p0A}VnPu8V>FKlobGy| zK!5AX0(A}q{3#?058)4*D6lyRP~6`#hEe>@i%FGyYnP;cc9MmE*#1J(d6b zx41tpIprz59n9kJX%0Oc)^Ye6hd<}=w;aC5VJn9RI6TSWB@T!1MaZ2Tev`u~VsXM9 zTjSItqC98MoHKcmIF=Y^awW~vR%sQot*v# za0>JA7;ukP#7UUG#vr%=579LzVSV%2V#JAl+DX(t`clN1|GeW~U|ftFg3dcQTgQR6 zoYQ-G$7bs{(3&`X&3Hr**t!hx7kJ2iFQ5F`nhUhwaC&|YA_#1K1^7oiB)=hv{ywMY z^BpN$M}Zc_W03BV);SSD5dHNd#L15FK}KfxR8frjG0{uEt<%}M3DIMC25JvITz&&J zfP`=sakA5yM5i?|(ZzTxRPF(t%Kt5=PNDWk;hT(pXZ+qPor({i+>RK(1^rnjKRQ1C z3i>k82lD?_CXY@wM@aU;f#hiqW5B8{c|711Te?UI;3g+eWtH}ZK^%$wh5 zbRpFTw_<`A$p4=RJ7!(hi1MQFFN4^5dk{M9a|g2jDd?2nG8oB)XB>C5wee?@Hz!yTF z!dgbxnygsr=VF6Rlu;16}D)V%YUT?A(b0AE>_aK!0EWJ2sL}!D@}q z=R)N1LF`EL%Kbreo)=dtx)xVCg03YJouWDP)3rcex#X#FeBbK_PBA*}bNTBfob%y6 zM>^Vcu=7p{7b`kKlO{^urBanwzpGI{EtbSdu~ds=J)C#qh|R%2K5&d3YI_C354<@_ z9~F4iIafSn{IDJ1{dQt%i4L!W#7e?j$ndT;?dV6xoX%tXA}I0qkPULU9`<-8OmCJ(+PXuBH6O`>6q%cy+JZ!fxDqQY@)ZIBx=FPEouF|ulbMWY1Nv2+aVEK;U*g%1gd=>M zx%VTXO7XEuzaSKe^}n{%Q4#Pua0*|6&mUyYS0C3kv-DY`bJ!k-ennNSBNT8~kXzyl z5-W1{<%nYZ^VERXh{o+Q&j)Z5Z2;~q4g|QX>@HJrHL89aYTO2@j9;l5Cy`lkk7+-m z!-0!iMU$c%T?f3Twcn+waedrqDmtqSGLFQ1RpT6TgEjHEV&9dm{e+~vKz&Wn zwHPqy*I}iI(eRMBq&AedI8=ozkgF=BJQTRw?~@o6&8AYytEjJq1sw+cI;Rfb&+t*9 zAL0Pd@2+-{0+0EsgJiJ^Z(JvfH+5g|ysVw(xhr{xt*n5qz~CtB#!_C~J$wzsr;K1d zCs~SY6KK@T3+i_Idb*y1-s_pit{QwL&BM18_R3=p=}Kq-QXW3g(KlN$V*j5(+Ah%6 zmQQOp-8XJZGLPw@U&sSSO8}uyzq_I!pg2#^sk*F)6Qt!A?hEh=eJqOgz@bQAki~w1 zOG$-@496q%#W{q4V!ozx3-VLwi}L^#G;&$!i*pSDoqT-*5f+kAiS~PRW+3#%Ie~!r zTu=Bf5Fy0fZ{w&s3rUV6Rt1HUn=C195%!YMyCt^xfI)PQ~E-Gu>0(V zd~*74fsv@6JzJL)@U$?@4T$;^_TE9{LPGmXz8Bvo1P9UI1bX85k$j(j;os5Um837u zaRgM7^po@70Y2qFl{3Wuk3rnKpDdq@AB8S~ zd+~FDw*R31*U+UVBJ{;RUi|+edMR@V)c+}nbVWN#1)t?#XG`LK?a1RgA&1ii5c+YL zq%Y24HLky2C~_Frg8>$yG$LR?bfF+FlG zPd`72F7%VfABeM%2$Q-9<41VH>-f1Ze)TNTP$K*n`J(e>;X700ule)!NNrMs6X^> 16) & 0xffff), lsl #16 ;\ + movk reg, #(((val) >> 32) & 0xffff), lsl #32 ;\ + movk reg, #(((val) >> 48) & 0xffff), lsl #48 + + .text + +/* =================================================================== + * uint64_t fm_int_math(uint64_t iters) + * + * Four largely independent accumulator chains to expose instruction-level + * parallelism, mixed with high-latency serialising ops (udiv/sdiv) and the + * bit-manipulation instructions. Returns a checksum so the compiler and the + * driver cannot elide the work. + * =================================================================== */ +FN_BEGIN(fm_int_math) + cbz x0, .Lim_zero + + MOV64(x1, 0x9E3779B97F4A7C15) /* a */ + MOV64(x2, 0xBF58476D1CE4E5B9) /* b */ + MOV64(x3, 0x94D049BB133111EB) /* c */ + MOV64(x4, 0x2545F4914F6CDD1D) /* d */ + MOV64(x5, 0x00000000DEADBEEF) /* odd multiplier, never zero */ + mov x6, x0 /* trip count */ + +.Lim_loop: + /* four independent multiply-accumulate chains */ + madd x1, x1, x5, x2 + madd x2, x2, x5, x3 + madd x3, x3, x5, x4 + madd x4, x4, x5, x1 + + /* cross-mix with shifts and logic ops (free shifter operands) */ + eor x1, x1, x3, lsr #29 + eor x2, x2, x4, lsl #17 + eor x3, x3, x1, ror #31 + bic x4, x4, x2, asr #7 + + /* wide multiplies: umulh/smulh are the long-latency multiplier path */ + umulh x9, x1, x3 + smulh x10, x2, x4 + add x1, x1, x9 + add x2, x2, x10 + + /* bit manipulation */ + rbit x11, x1 + clz x12, x2 + rev x13, x3 + eor x4, x4, x11 + add x4, x4, x12 + eor x1, x1, x13 + + /* division: fully serialising, ~10-20 cycle latency, not pipelined */ + orr x14, x5, #1 /* guarantee a non-zero divisor */ + udiv x15, x1, x14 + sdiv x16, x2, x14 + msub x3, x15, x14, x3 + add x4, x4, x16 + + /* bitfield ops */ + ror x2, x2, #11 + extr x1, x1, x2, #23 + + subs x6, x6, #1 + b.ne .Lim_loop + + eor x0, x1, x2 + eor x0, x0, x3 + eor x0, x0, x4 + ret + +.Lim_zero: + mov x0, xzr + ret +FN_END(fm_int_math) + + +/* =================================================================== + * uint64_t fm_fp_math(uint64_t iters) + * + * Double-precision scalar FP. Four fmadd chains cover the pipelined + * multiply-add path; fdiv and fsqrt cover the non-pipelined divide/sqrt unit, + * which is usually the real differentiator between cores. + * Returns the result bit-cast to u64. + * =================================================================== */ +FN_BEGIN(fm_fp_math) + cbz x0, .Lfp_zero + mov x6, x0 + + /* Constants come from a table rather than fmov immediates: the AArch64 + * 8-bit FP immediate can only encode a narrow set of values, and + * several of the ones we want fall outside it. */ + adr x7, .Lfp_consts + ldp d0, d1, [x7] /* a = 1.5, b = 2.5 */ + ldp d2, d3, [x7, #16] /* c = 3.5, d = 0.5 */ + ldp d4, d5, [x7, #32] /* mul, small addend */ + ldp d6, d7, [x7, #48] /* 2.0, 1.0 */ + +.Lfp_loop: + /* four independent fused multiply-add chains */ + fmadd d0, d0, d4, d5 + fmadd d1, d1, d4, d5 + fmadd d2, d2, d4, d5 + fmadd d3, d3, d4, d5 + + /* keep the accumulators bounded so they never reach inf/NaN */ + fmin d0, d0, d6 + fmin d1, d1, d6 + fmin d2, d2, d6 + fmin d3, d3, d6 + + /* square root: long latency, low throughput */ + fsqrt d16, d0 + fsqrt d17, d1 + fadd d2, d2, d16 + fadd d3, d3, d17 + + /* divide: the other long-latency unit */ + fadd d18, d2, d7 /* divisor >= 1, never zero */ + fdiv d19, d7, d18 + fadd d0, d0, d19 + + fadd d20, d3, d7 + fdiv d21, d7, d20 + fadd d1, d1, d21 + + /* abs/neg/compare-select: cheap ops to balance the mix */ + fabs d2, d2 + fneg d22, d3 + fabs d3, d22 + fmax d3, d3, d7 + + subs x6, x6, #1 + b.ne .Lfp_loop + + fadd d0, d0, d1 + fadd d2, d2, d3 + fadd d0, d0, d2 + fmov x0, d0 + ret + +.Lfp_zero: + mov x0, xzr + ret +FN_END(fm_fp_math) + + .p2align 4 +.Lfp_consts: + .double 1.5, 2.5 + .double 3.5, 0.5 + .double 1.0625, 0.0009765625 + .double 2.0, 1.0 + + +/* =================================================================== + * uint64_t fm_primes(uint64_t limit, uint8_t *sieve) + * + * Sieve of Eratosthenes over [0, limit). The caller supplies `limit` bytes of + * scratch; this routine clears it itself, so the clearing pass is part of the + * measured work (as it would be in any real use). Returns the prime count. + * + * Strided stores over a buffer larger than L1 make this a memory-hierarchy + * test as much as an arithmetic one. + * =================================================================== */ +FN_BEGIN(fm_primes) + cmp x0, #2 + b.lo .Lpr_none + + mov x2, x1 /* sieve base */ + mov x3, x0 /* limit */ + + /* zero the sieve, 32 bytes per iteration */ + movi v0.16b, #0 + mov x4, xzr + and x5, x3, #~31 /* bulk portion */ +.Lpr_clear32: + cmp x4, x5 + b.hs .Lpr_clear1 + add x6, x2, x4 + stp q0, q0, [x6] + add x4, x4, #32 + b .Lpr_clear32 +.Lpr_clear1: + cmp x4, x3 + b.hs .Lpr_clear_done + strb wzr, [x2, x4] + add x4, x4, #1 + b .Lpr_clear1 +.Lpr_clear_done: + + /* mark 0 and 1 as composite */ + mov w6, #1 + strb w6, [x2] + strb w6, [x2, #1] + + /* outer loop: i = 2; i*i < limit; i++ */ + mov x7, #2 +.Lpr_outer: + mul x9, x7, x7 + cmp x9, x3 + b.hs .Lpr_count + + ldrb w10, [x2, x7] + cbnz w10, .Lpr_outer_next /* already composite, skip */ + + /* inner loop: mark multiples starting at i*i, stride i */ + mov x11, x9 +.Lpr_inner: + cmp x11, x3 + b.hs .Lpr_outer_next + strb w6, [x2, x11] + add x11, x11, x7 + b .Lpr_inner + +.Lpr_outer_next: + add x7, x7, #1 + b .Lpr_outer + + /* count the survivors */ +.Lpr_count: + mov x0, xzr /* count */ + mov x4, #2 +.Lpr_count_loop: + cmp x4, x3 + b.hs .Lpr_done + ldrb w10, [x2, x4] + cmp w10, #0 + cinc x0, x0, eq + add x4, x4, #1 + b .Lpr_count_loop + +.Lpr_done: + ret + +.Lpr_none: + mov x0, xzr + ret +FN_END(fm_primes) + + +/* =================================================================== + * uint64_t fm_simd(uint64_t iters, void *buf) + * + * "Extended instructions": the ASIMD/NEON unit, which is architecturally + * mandatory on AArch64 and therefore safe to use without runtime feature + * detection. (AES, SHA and DotProd are all *optional* extensions; using them + * unguarded would fault on cores that lack them, so they are deliberately + * avoided here.) + * + * buf must be at least 128 bytes and 16-byte aligned. Returns a checksum. + * =================================================================== */ +FN_BEGIN(fm_simd) + cbz x0, .Lsd_zero + mov x6, x0 + + /* seed eight vectors from the scratch buffer */ + ldp q0, q1, [x1] + ldp q2, q3, [x1, #32] + ldp q4, q5, [x1, #64] + ldp q6, q7, [x1, #96] + + movi v28.4s, #3 + movi v29.4s, #7 + movi v30.16b, #0x5A + + /* byte-permute table for tbl: reverses the 16 byte lanes */ + adr x9, .Lsd_perm + ldr q31, [x9] + + /* float accumulator and operands: derived from the integer seeds by + * conversion, so the FP pipeline sees real finite values rather than + * reinterpreted integer bit patterns (which would be NaNs/denormals + * and would measure the slow path, not the common one). */ + movi v22.4s, #0 + scvtf v26.4s, v0.4s + scvtf v27.4s, v1.4s + +.Lsd_loop: + /* integer SIMD: multiply-accumulate across four independent vectors */ + mla v0.4s, v1.4s, v28.4s + mla v1.4s, v2.4s, v29.4s + mla v2.4s, v3.4s, v28.4s + mla v3.4s, v0.4s, v29.4s + + /* saturating and halving arithmetic */ + sqadd v4.4s, v4.4s, v0.4s + uhadd v5.4s, v5.4s, v1.4s + srhadd v6.4s, v6.4s, v2.4s + sqsub v7.4s, v7.4s, v3.4s + + /* widening ops: 16->32 bit lane expansion */ + umull v16.4s, v0.4h, v1.4h + umull2 v17.4s, v0.8h, v1.8h + saddw v2.4s, v2.4s, v16.4h + ssubw2 v3.4s, v3.4s, v17.8h + + /* pairwise reduction */ + uaddlp v18.2d, v4.4s + addp v19.4s, v5.4s, v6.4s + add v4.4s, v4.4s, v19.4s + + /* table lookup: the byte-permute network */ + tbl v20.16b, {v0.16b}, v31.16b + eor v1.16b, v1.16b, v20.16b + + /* shifts, logic, min/max, reverse */ + shl v21.4s, v2.4s, #3 + usra v21.4s, v2.4s, #29 + orr v2.16b, v2.16b, v21.16b + bic v3.16b, v3.16b, v30.16b + umax v5.4s, v5.4s, v0.4s + umin v6.4s, v6.4s, v1.4s + rev32 v7.16b, v7.16b + + /* single-precision float SIMD: 4-wide fmla, plus the reciprocal and + * rsqrt estimate instructions that shader-style code leans on */ + /* v16/v17 are dead after the widening ops above, so they are reused + * here as scratch rather than touching the callee-saved v8-v15 bank. */ + fmla v22.4s, v26.4s, v27.4s + frecpe v23.4s, v26.4s + frsqrte v16.4s, v22.4s + fadd v22.4s, v22.4s, v23.4s + fmul v26.4s, v26.4s, v16.4s + + /* population count and leading-zero count */ + cnt v24.16b, v0.16b + uaddlv h25, v24.8b + clz v17.4s, v1.4s + add v0.4s, v0.4s, v17.4s + + /* fold the reduction results back in so nothing is dead code */ + add v4.2d, v4.2d, v18.2d + + subs x6, x6, #1 + b.ne .Lsd_loop + + /* horizontal fold to a single 64-bit checksum */ + eor v0.16b, v0.16b, v1.16b + eor v2.16b, v2.16b, v3.16b + eor v4.16b, v4.16b, v5.16b + eor v6.16b, v6.16b, v7.16b + eor v0.16b, v0.16b, v2.16b + eor v4.16b, v4.16b, v6.16b + eor v0.16b, v0.16b, v4.16b + addv s0, v0.4s + fmov w0, s0 + ret + +.Lsd_zero: + mov x0, xzr + ret +FN_END(fm_simd) + + .p2align 4 +.Lsd_perm: + .byte 15,14,13,12,11,10,9,8,7,6,5,4,3,2,1,0 + + +/* =================================================================== + * uint64_t fm_compress(const uint8_t *src, uint64_t len, uint32_t *ht) + * + * The match-finding inner loop of an LZ77 compressor (the LZ4 fast strategy): + * hash the next 4 bytes, probe a single-entry-per-bucket table, verify, then + * extend the match. This is where real compressors spend their time - it is + * branch-heavy with a data-dependent, cache-missing table probe. + * + * ht must hold 1<<16 uint32_t (256 KiB); this routine clears it itself. + * Returns the encoded size in bytes. + * =================================================================== */ +FN_BEGIN(fm_compress) + stp x29, x30, [sp, #-96]! + mov x29, sp + stp x19, x20, [sp, #16] + stp x21, x22, [sp, #32] + stp x23, x24, [sp, #48] + stp x25, x26, [sp, #64] + stp x27, x28, [sp, #80] + + mov x19, x0 /* src */ + mov x20, x1 /* len */ + mov x21, x2 /* ht */ + + /* clear the hash table: 1<<16 entries * 4 bytes = 262144 bytes */ + movi v0.16b, #0 + mov x9, xzr + MOV64(x10, 262144) +.Lcm_clear: + add x11, x21, x9 + stp q0, q0, [x11] + stp q0, q0, [x11, #32] + add x9, x9, #64 + cmp x9, x10 + b.lo .Lcm_clear + + cmp x20, #16 + b.lo .Lcm_tiny + + mov x22, x19 /* ip */ + mov x23, x19 /* anchor */ + add x24, x19, x20 /* end */ + sub x25, x24, #12 /* mflimit */ + mov x26, xzr /* outsize */ + + MOV64(x27, 2654435761) /* Knuth multiplicative hash */ + +.Lcm_loop: + cmp x22, x25 + b.hs .Lcm_flush + + ldr w9, [x22] /* seq = load32(ip) */ + mul w10, w9, w27 + lsr w10, w10, #16 /* h = (seq * prime) >> 16 */ + + ldr w11, [x21, x10, lsl #2] /* ref_off = ht[h] */ + sub x12, x22, x19 /* cur_off = ip - src */ + str w12, [x21, x10, lsl #2] /* ht[h] = cur_off */ + + add x13, x19, x11 /* ref = src + ref_off */ + cmp x13, x22 + b.hs .Lcm_no_match /* ref must be strictly behind ip */ + + sub x14, x22, x13 /* distance */ + MOV64(x15, 65536) + cmp x14, x15 + b.hs .Lcm_no_match /* 16-bit offset window */ + + ldr w16, [x13] + cmp w16, w9 + b.ne .Lcm_no_match + + /* match confirmed: extend it byte by byte */ + mov x28, #4 /* ml */ +.Lcm_extend: + add x9, x22, x28 + cmp x9, x24 + b.hs .Lcm_emit + ldrb w10, [x22, x28] + ldrb w11, [x13, x28] + cmp w10, w11 + b.ne .Lcm_emit + add x28, x28, #1 + b .Lcm_extend + +.Lcm_emit: + /* token(1) + offset(2) + literals + varint extensions */ + sub x9, x22, x23 /* literal run length */ + add x26, x26, x9 + add x26, x26, #3 + cmp x9, #15 + cinc x26, x26, hs /* literal-length extension byte */ + cmp x28, #19 + cinc x26, x26, hs /* match-length extension byte */ + + add x22, x22, x28 + mov x23, x22 + b .Lcm_loop + +.Lcm_no_match: + add x22, x22, #1 + b .Lcm_loop + +.Lcm_flush: + /* trailing literals */ + sub x9, x24, x23 + add x26, x26, x9 + add x26, x26, #1 + mov x0, x26 + b .Lcm_ret + +.Lcm_tiny: + add x0, x20, #1 + +.Lcm_ret: + ldp x27, x28, [sp, #80] + ldp x25, x26, [sp, #64] + ldp x23, x24, [sp, #48] + ldp x21, x22, [sp, #32] + ldp x19, x20, [sp, #16] + ldp x29, x30, [sp], #96 + ret +FN_END(fm_compress) + + +/* =================================================================== + * uint64_t fm_chacha20(uint8_t *buf, uint64_t len, const uint8_t key[32], + * uint64_t rounds) + * + * ChaCha20 stream cipher, NEON, four 128-bit state rows. Chosen over AES + * deliberately: the ARMv8 AES extension is optional, so an AES-instruction + * benchmark would fault on cores without it. ChaCha20 needs only baseline + * ASIMD and is a real, widely deployed cipher (TLS, WireGuard, SSH). + * + * len is rounded down to a multiple of 64. `rounds` = number of passes over + * the buffer. Returns a checksum of the keystream output. + * =================================================================== */ + +/* rotate each 32-bit lane left by n, via shl + shift-right-and-insert */ +#define VROTL(vd, vs, vt, n) \ + shl vt##.4s, vs##.4s, #(n) ;\ + sri vt##.4s, vs##.4s, #(32 - (n)) ;\ + mov vd##.16b, vt##.16b + +/* one ChaCha quarter-round over rows a,b,c,d using v24 as scratch */ +#define QROUND(a, b, c, d) \ + add a##.4s, a##.4s, b##.4s ;\ + eor d##.16b, d##.16b, a##.16b ;\ + rev32 d##.8h, d##.8h ;\ + add c##.4s, c##.4s, d##.4s ;\ + eor b##.16b, b##.16b, c##.16b ;\ + VROTL(b, b, v24, 12) ;\ + add a##.4s, a##.4s, b##.4s ;\ + eor d##.16b, d##.16b, a##.16b ;\ + VROTL(d, d, v24, 8) ;\ + add c##.4s, c##.4s, d##.4s ;\ + eor b##.16b, b##.16b, c##.16b ;\ + VROTL(b, b, v24, 7) + +FN_BEGIN(fm_chacha20) + and x1, x1, #~63 /* whole 64-byte blocks only */ + cbz x1, .Lcc_zero + cbz x3, .Lcc_zero + + stp x29, x30, [sp, #-32]! + mov x29, sp + stp x19, x20, [sp, #16] + + mov x19, x0 /* buf */ + mov x20, x1 /* len */ + + /* v4..v7 hold the base state */ + adr x9, .Lcc_sigma + ldr q4, [x9] /* "expand 32-byte k" */ + ldp q5, q6, [x2] /* key[0..31] */ + movi v7.4s, #0 /* counter || nonce */ + + movi v25.16b, #0 /* running checksum */ + mov x10, xzr /* counter value */ + +.Lcc_pass: + mov x11, xzr /* byte offset into buf */ + +.Lcc_block: + /* working state = base state, with the block counter in lane 0 of v7 */ + mov v0.16b, v4.16b + mov v1.16b, v5.16b + mov v2.16b, v6.16b + mov v3.16b, v7.16b + mov v3.s[0], w10 + + /* keep originals for the final feed-forward add */ + mov v16.16b, v0.16b + mov v17.16b, v1.16b + mov v18.16b, v2.16b + mov v19.16b, v3.16b + + mov w12, #10 /* 10 double rounds = 20 rounds */ +.Lcc_rounds: + /* column round */ + QROUND(v0, v1, v2, v3) + + /* rotate lanes to form the diagonals */ + ext v1.16b, v1.16b, v1.16b, #4 + ext v2.16b, v2.16b, v2.16b, #8 + ext v3.16b, v3.16b, v3.16b, #12 + + /* diagonal round */ + QROUND(v0, v1, v2, v3) + + /* undo the lane rotation */ + ext v1.16b, v1.16b, v1.16b, #12 + ext v2.16b, v2.16b, v2.16b, #8 + ext v3.16b, v3.16b, v3.16b, #4 + + subs w12, w12, #1 + b.ne .Lcc_rounds + + /* feed-forward: keystream = working + original */ + add v0.4s, v0.4s, v16.4s + add v1.4s, v1.4s, v17.4s + add v2.4s, v2.4s, v18.4s + add v3.4s, v3.4s, v19.4s + + /* XOR the keystream into the buffer */ + add x13, x19, x11 + ldp q20, q21, [x13] + ldp q22, q23, [x13, #32] + eor v20.16b, v20.16b, v0.16b + eor v21.16b, v21.16b, v1.16b + eor v22.16b, v22.16b, v2.16b + eor v23.16b, v23.16b, v3.16b + stp q20, q21, [x13] + stp q22, q23, [x13, #32] + + /* accumulate a checksum of the keystream */ + eor v25.16b, v25.16b, v0.16b + eor v25.16b, v25.16b, v3.16b + + add x10, x10, #1 + add x11, x11, #64 + cmp x11, x20 + b.lo .Lcc_block + + subs x3, x3, #1 + b.ne .Lcc_pass + + addv s25, v25.4s + fmov w0, s25 + + ldp x19, x20, [sp, #16] + ldp x29, x30, [sp], #32 + ret + +.Lcc_zero: + mov x0, xzr + ret +FN_END(fm_chacha20) + + .p2align 4 +.Lcc_sigma: + .word 0x61707865, 0x3320646e, 0x79622d32, 0x6b206574 + + +/* =================================================================== + * uint64_t fm_physics(double *bodies, uint64_t n, uint64_t steps) + * + * Direct-summation N-body gravity, O(n^2) per step, double precision. + * Layout per body, 8 doubles (64 bytes, one cache line): + * [0]=x [1]=y [2]=z [3]=mass [4]=vx [5]=vy [6]=vz [7]=pad + * + * The 1/sqrt is done with a real fsqrt+fdiv rather than the frsqrte estimate, + * so this exercises the divide/sqrt unit the way physics code actually does. + * Returns a checksum bit-cast from the final velocity sum. + * =================================================================== */ +FN_BEGIN(fm_physics) + cbz x1, .Lph_zero + cbz x2, .Lph_zero + + stp x29, x30, [sp, #-64]! + mov x29, sp + stp x19, x20, [sp, #16] + stp x21, x22, [sp, #32] + stp x23, x24, [sp, #48] + + mov x19, x0 /* bodies */ + mov x20, x1 /* n */ + mov x21, x2 /* steps */ + + adr x9, .Lph_consts + ldp d28, d29, [x9] /* dt, eps^2 */ + ldr d30, [x9, #16] /* 1.0 */ + +.Lph_step: + mov x22, xzr /* i */ + +.Lph_body_i: + lsl x9, x22, #6 /* i * 64 */ + add x23, x19, x9 /* &bodies[i] */ + + ldp d0, d1, [x23] /* xi, yi */ + ldr d2, [x23, #16] /* zi */ + + movi d16, #0 /* ax */ + movi d17, #0 /* ay */ + movi d18, #0 /* az */ + + mov x24, xzr /* j */ + mov x10, x19 /* &bodies[j] */ + +.Lph_body_j: + ldp d3, d4, [x10] /* xj, yj */ + ldp d5, d6, [x10, #16] /* zj, mj */ + + fsub d3, d3, d0 /* dx */ + fsub d4, d4, d1 /* dy */ + fsub d5, d5, d2 /* dz */ + + /* d2 = dx*dx + dy*dy + dz*dz + eps^2 (always >= eps^2, never zero, + * so the i==j self-term is finite and contributes exactly 0 below) */ + fmul d7, d3, d3 + fmadd d7, d4, d4, d7 + fmadd d7, d5, d5, d7 + fadd d7, d7, d29 + + fsqrt d19, d7 /* r */ + fdiv d20, d30, d19 /* 1/r */ + fmul d21, d20, d20 /* 1/r^2 */ + fmul d21, d21, d20 /* 1/r^3 */ + fmul d21, d21, d6 /* m/r^3 */ + + fmadd d16, d3, d21, d16 /* ax += dx * m/r^3 */ + fmadd d17, d4, d21, d17 + fmadd d18, d5, d21, d18 + + add x10, x10, #64 + add x24, x24, #1 + cmp x24, x20 + b.lo .Lph_body_j + + /* v += a * dt */ + ldp d22, d23, [x23, #32] + ldr d24, [x23, #48] + fmadd d22, d16, d28, d22 + fmadd d23, d17, d28, d23 + fmadd d24, d18, d28, d24 + stp d22, d23, [x23, #32] + str d24, [x23, #48] + + add x22, x22, #1 + cmp x22, x20 + b.lo .Lph_body_i + + /* second pass: x += v * dt (positions updated only after all forces) */ + mov x22, xzr + mov x10, x19 +.Lph_integrate: + ldp d0, d1, [x10] + ldr d2, [x10, #16] + ldp d22, d23, [x10, #32] + ldr d24, [x10, #48] + fmadd d0, d22, d28, d0 + fmadd d1, d23, d28, d1 + fmadd d2, d24, d28, d2 + stp d0, d1, [x10] + str d2, [x10, #16] + add x10, x10, #64 + add x22, x22, #1 + cmp x22, x20 + b.lo .Lph_integrate + + subs x21, x21, #1 + b.ne .Lph_step + + /* checksum: sum of all velocity components */ + movi d0, #0 + mov x22, xzr + mov x10, x19 +.Lph_sum: + ldp d22, d23, [x10, #32] + ldr d24, [x10, #48] + fadd d0, d0, d22 + fadd d0, d0, d23 + fadd d0, d0, d24 + add x10, x10, #64 + add x22, x22, #1 + cmp x22, x20 + b.lo .Lph_sum + + fmov x0, d0 + + ldp x23, x24, [sp, #48] + ldp x21, x22, [sp, #32] + ldp x19, x20, [sp, #16] + ldp x29, x30, [sp], #64 + ret + +.Lph_zero: + mov x0, xzr + ret +FN_END(fm_physics) + + .p2align 4 +.Lph_consts: + .double 0.0078125 /* dt */ + .double 0.0625 /* eps^2 */ + .double 1.0 + + +/* =================================================================== + * uint64_t fm_sort(uint32_t *a, uint64_t n) + * + * In-place heapsort. Chosen over quicksort because it needs no recursion or + * explicit stack, yet is aggressively branch-unpredictable and touches memory + * in a scattered pattern - it stresses the branch predictor and the cache + * hierarchy, which is what a sort benchmark should measure. + * + * Returns an order-sensitive checksum, which also verifies the sort. + * =================================================================== */ +FN_BEGIN(fm_sort) + cmp x1, #2 + b.lo .Lst_trivial + + stp x29, x30, [sp, #-48]! + mov x29, sp + stp x19, x20, [sp, #16] + stp x21, x22, [sp, #32] + + mov x19, x0 /* a */ + mov x20, x1 /* n */ + + /* ---- build the max-heap: for i = n/2 - 1 down to 0 ---- */ + lsr x21, x20, #1 /* i = n/2 */ +.Lst_build: + cbz x21, .Lst_extract + sub x21, x21, #1 /* i-- */ + mov x0, x21 /* root */ + mov x1, x20 /* end */ + bl .Lst_siftdown + cbnz x21, .Lst_build + + /* ---- extract: for end = n-1 down to 1 ---- */ +.Lst_extract: + sub x22, x20, #1 /* end = n-1 */ +.Lst_extract_loop: + cbz x22, .Lst_checksum + + /* swap a[0] and a[end] */ + ldr w9, [x19] + ldr w10, [x19, x22, lsl #2] + str w10, [x19] + str w9, [x19, x22, lsl #2] + + mov x0, xzr /* root = 0 */ + mov x1, x22 /* end = end */ + bl .Lst_siftdown + + sub x22, x22, #1 + b .Lst_extract_loop + + /* ---- order-sensitive checksum ---- */ +.Lst_checksum: + mov x0, xzr + mov x9, xzr +.Lst_cksum_loop: + ldr w10, [x19, x9, lsl #2] + eor x0, x0, x10 + ror x0, x0, #7 + add x0, x0, x10 + add x9, x9, #1 + cmp x9, x20 + b.lo .Lst_cksum_loop + + ldp x21, x22, [sp, #32] + ldp x19, x20, [sp, #16] + ldp x29, x30, [sp], #48 + ret + +.Lst_trivial: + mov x0, xzr + cbz x1, .Lst_trivial_ret + ldr w0, [x0] +.Lst_trivial_ret: + ret + + /* ---- local helper: siftdown(root = x0, end = x1) + * clobbers x9-x15 only; x19 (base) is live across the call. ---- */ +.Lst_siftdown: + mov x11, x0 /* root */ +.Lst_sift_loop: + lsl x12, x11, #1 + add x12, x12, #1 /* child = 2*root + 1 */ + cmp x12, x1 + b.hs .Lst_sift_done /* no children */ + + /* pick the larger of the two children */ + add x13, x12, #1 /* child + 1 */ + cmp x13, x1 + b.hs .Lst_sift_have_child + ldr w14, [x19, x12, lsl #2] + ldr w15, [x19, x13, lsl #2] + cmp w15, w14 + csel x12, x13, x12, hi + +.Lst_sift_have_child: + ldr w14, [x19, x11, lsl #2] /* a[root] */ + ldr w15, [x19, x12, lsl #2] /* a[child] */ + cmp w14, w15 + b.hs .Lst_sift_done /* heap property holds */ + + /* swap and descend */ + str w15, [x19, x11, lsl #2] + str w14, [x19, x12, lsl #2] + mov x11, x12 + b .Lst_sift_loop + +.Lst_sift_done: + ret +FN_END(fm_sort) + + +/* =================================================================== + * uint64_t fm_chase(void **ptrs, uint64_t steps) + * + * Pointer chase around a randomised cycle. Every load depends on the previous + * one, so nothing can be prefetched, overlapped or reordered - this measures + * the pure serial latency of the memory hierarchy, which is the single + * hardest thing for a wide out-of-order core to hide. It is the truest + * "single-threaded" test in the suite. + * =================================================================== */ +FN_BEGIN(fm_chase) + cbz x1, .Lch_zero + mov x2, x0 /* p = ptrs */ + mov x3, x1 + +.Lch_loop: + ldr x2, [x2] + subs x3, x3, #1 + b.ne .Lch_loop + + sub x0, x2, x0 /* final offset, keeps p live */ + ret + +.Lch_zero: + mov x0, xzr + ret +FN_END(fm_chase) + + +#if defined(__ELF__) + .section .note.GNU-stack, "", %progbits +#endif diff --git a/src/fossmark_x86_64.S b/src/fossmark_x86_64.S new file mode 100644 index 0000000..4332ecf --- /dev/null +++ b/src/fossmark_x86_64.S @@ -0,0 +1,987 @@ +/* + * fossmark_x86_64.S - x86-64 (AMD64) CPU benchmark kernels + * + * The AMD64 counterpart to fossmark.S. Same nine routines, same contract: each + * is a pure function of its arguments under the System V AMD64 ABI, contains no + * syscalls, no libc calls and no external data relocations, so it assembles and + * runs on Linux (ELF), macOS (Mach-O) and the BSDs. The portable C driver in + * main.c is shared unchanged between this file and the AArch64 one. + * + * Only baseline instructions are used: general-purpose AMD64 plus SSE2, which + * is architecturally mandatory on x86-64. The optional extensions (SSE4, AVX, + * FMA, AES-NI, POPCNT/BMI) are deliberately avoided - using them unguarded + * would fault (#UD) on cores that lack them - so, exactly as the NEON file + * sticks to mandatory ASIMD and shuns the optional AES/SHA/DotProd, this file + * sticks to mandatory SSE2 and shuns everything above it. FMA in particular is + * optional here, so every fused multiply-add is written as a separate multiply + * and add. + * + * Register conventions (System V AMD64): + * integer args rdi, rsi, rdx, rcx, r8, r9 (return in rax) + * callee-saved rbx, rbp, r12, r13, r14, r15 (saved when used) + * all of xmm0-15 are caller-saved, so no vector register need be preserved. + */ + + .intel_syntax noprefix + +#if defined(__APPLE__) +# define SYM(name) _##name +#else +# define SYM(name) name +#endif + +#if defined(__ELF__) +# define FN_BEGIN(name) .p2align 4 ; .globl SYM(name) ; .type SYM(name), @function ; SYM(name): +# define FN_END(name) .size SYM(name), . - SYM(name) +#else +# define FN_BEGIN(name) .p2align 4 ; .globl SYM(name) ; SYM(name): +# define FN_END(name) +#endif + + .text + +/* =================================================================== + * uint64_t fm_int_math(uint64_t iters) [rdi = iters] + * + * Four independent multiply-accumulate chains for instruction-level + * parallelism, mixed with the long-latency serialising ops (mul/div) and + * bit-manipulation. Returns a checksum so nothing can be elided. + * =================================================================== */ +FN_BEGIN(fm_int_math) + test rdi, rdi + jz .Lim_zero + + mov r8, 0x9E3779B97F4A7C15 /* a */ + mov r9, 0xBF58476D1CE4E5B9 /* b */ + mov r10, 0x94D049BB133111EB /* c */ + mov r11, 0x2545F4914F6CDD1D /* d */ + mov rsi, 0x00000000DEADBEEF /* odd multiplier, never zero */ + +.Lim_loop: + /* four independent multiply-accumulate chains */ + imul r8, rsi + add r8, r9 + imul r9, rsi + add r9, r10 + imul r10, rsi + add r10, r11 + imul r11, rsi + add r11, r8 + + /* cross-mix with shifts and logic ops */ + mov rax, r10 + shr rax, 29 + xor r8, rax + mov rax, r11 + shl rax, 17 + xor r9, rax + mov rax, r8 + ror rax, 31 + xor r10, rax + mov rax, r9 + sar rax, 7 + not rax + and r11, rax /* r11 &= ~(r9 >> 7 arith) */ + + /* wide multiplies: the long-latency 128-bit multiplier path */ + mov rax, r8 + mul r10 /* rdx:rax = a*c, high in rdx */ + add r8, rdx + mov rax, r9 + imul r11 /* rdx:rax = b*d signed, high rdx */ + add r9, rdx + + /* bit manipulation: byte reverse */ + mov rax, r10 + bswap rax + xor r8, rax + mov rax, r11 + bswap rax + xor r9, rax + + /* division: fully serialising, not pipelined */ + mov rcx, rsi + or rcx, 1 /* guarantee a non-zero divisor */ + mov rax, r8 + xor edx, edx + div rcx /* rax = a / rcx (unsigned) */ + imul rax, rcx + sub r10, rax /* c -= (a/div)*div */ + mov rax, r9 + cqo + idiv rcx /* rax = b / rcx (signed) */ + add r11, rax + + /* bitfield ops */ + ror r9, 11 + shld r8, r9, 23 + + dec rdi + jnz .Lim_loop + + mov rax, r8 + xor rax, r9 + xor rax, r10 + xor rax, r11 + ret + +.Lim_zero: + xor eax, eax + ret +FN_END(fm_int_math) + + +/* =================================================================== + * uint64_t fm_fp_math(uint64_t iters) [rdi = iters] + * + * Double-precision scalar FP. Four multiply-add chains for the pipelined + * path; sqrtsd and divsd for the non-pipelined divide/sqrt unit that usually + * separates cores. Returns the result bit-cast to u64. + * =================================================================== */ +FN_BEGIN(fm_fp_math) + test rdi, rdi + jz .Lfp_zero + + movsd xmm0, [rip + .Lfp_consts + 0] /* a = 1.5 */ + movsd xmm1, [rip + .Lfp_consts + 8] /* b = 2.5 */ + movsd xmm2, [rip + .Lfp_consts + 16] /* c = 3.5 */ + movsd xmm3, [rip + .Lfp_consts + 24] /* d = 0.5 */ + movsd xmm4, [rip + .Lfp_consts + 32] /* mul */ + movsd xmm5, [rip + .Lfp_consts + 40] /* addend */ + movsd xmm6, [rip + .Lfp_consts + 48] /* 2.0 */ + movsd xmm7, [rip + .Lfp_consts + 56] /* 1.0 */ + +.Lfp_loop: + /* four independent multiply-add chains (no baseline FMA) */ + mulsd xmm0, xmm4 + addsd xmm0, xmm5 + mulsd xmm1, xmm4 + addsd xmm1, xmm5 + mulsd xmm2, xmm4 + addsd xmm2, xmm5 + mulsd xmm3, xmm4 + addsd xmm3, xmm5 + + /* keep the accumulators bounded so they never reach inf/NaN */ + minsd xmm0, xmm6 + minsd xmm1, xmm6 + minsd xmm2, xmm6 + minsd xmm3, xmm6 + + /* square root: long latency, low throughput */ + sqrtsd xmm8, xmm0 + sqrtsd xmm9, xmm1 + addsd xmm2, xmm8 + addsd xmm3, xmm9 + + /* divide: 1/(c+1), divisor >= 1 so never zero */ + movapd xmm10, xmm2 + addsd xmm10, xmm7 + movapd xmm11, xmm7 + divsd xmm11, xmm10 + addsd xmm0, xmm11 + + movapd xmm10, xmm3 + addsd xmm10, xmm7 + movapd xmm11, xmm7 + divsd xmm11, xmm10 + addsd xmm1, xmm11 + + /* abs/neg/max: cheap ops to balance the mix */ + andpd xmm2, [rip + .Lfp_absmask] /* fabs(c) */ + xorpd xmm3, [rip + .Lfp_signmask] /* fneg(d) */ + andpd xmm3, [rip + .Lfp_absmask] /* fabs() */ + maxsd xmm3, xmm7 + + dec rdi + jnz .Lfp_loop + + addsd xmm0, xmm1 + addsd xmm2, xmm3 + addsd xmm0, xmm2 + movq rax, xmm0 + ret + +.Lfp_zero: + xor eax, eax + ret +FN_END(fm_fp_math) + + .p2align 4 +.Lfp_consts: + .double 1.5, 2.5 + .double 3.5, 0.5 + .double 1.0625, 0.0009765625 + .double 2.0, 1.0 + .p2align 4 +.Lfp_absmask: + .quad 0x7fffffffffffffff, 0x7fffffffffffffff +.Lfp_signmask: + .quad 0x8000000000000000, 0x8000000000000000 + + +/* =================================================================== + * uint64_t fm_primes(uint64_t limit, uint8_t *sieve) [rdi, rsi] + * + * Sieve of Eratosthenes over [0, limit). The routine clears the caller's + * scratch itself, so the clearing pass counts as measured work. Strided stores + * over a buffer larger than L1 make this a memory-hierarchy test too. Returns + * the prime count. + * =================================================================== */ +FN_BEGIN(fm_primes) + cmp rdi, 2 + jb .Lpr_none + + /* zero the sieve, 32 bytes per iteration */ + pxor xmm0, xmm0 + xor rax, rax /* index */ + mov rcx, rdi + and rcx, -32 /* bulk portion */ +.Lpr_clear32: + cmp rax, rcx + jae .Lpr_clear1 + movdqu [rsi + rax], xmm0 + movdqu [rsi + rax + 16], xmm0 + add rax, 32 + jmp .Lpr_clear32 +.Lpr_clear1: + cmp rax, rdi + jae .Lpr_clear_done + mov byte ptr [rsi + rax], 0 + inc rax + jmp .Lpr_clear1 +.Lpr_clear_done: + + /* mark 0 and 1 as composite */ + mov byte ptr [rsi], 1 + mov byte ptr [rsi + 1], 1 + + /* outer loop: i = 2; i*i < limit; i++ */ + mov r8, 2 +.Lpr_outer: + mov rax, r8 + imul rax, r8 /* i*i */ + cmp rax, rdi + jae .Lpr_count + + movzx edx, byte ptr [rsi + r8] + test dl, dl + jnz .Lpr_outer_next /* already composite, skip */ + + /* inner loop: mark multiples starting at i*i, stride i */ + mov r9, rax /* j = i*i */ +.Lpr_inner: + cmp r9, rdi + jae .Lpr_outer_next + mov byte ptr [rsi + r9], 1 + add r9, r8 + jmp .Lpr_inner + +.Lpr_outer_next: + inc r8 + jmp .Lpr_outer + + /* count the survivors */ +.Lpr_count: + xor eax, eax /* count */ + mov r9, 2 +.Lpr_count_loop: + cmp r9, rdi + jae .Lpr_done + cmp byte ptr [rsi + r9], 0 + jne .Lpr_count_next + inc rax +.Lpr_count_next: + inc r9 + jmp .Lpr_count_loop + +.Lpr_done: + ret + +.Lpr_none: + xor eax, eax + ret +FN_END(fm_primes) + + +/* =================================================================== + * uint64_t fm_simd(uint64_t iters, void *buf) [rdi = iters, rsi = buf] + * + * "Extended instructions": the SSE2 unit, which is architecturally mandatory + * on x86-64 and therefore safe without runtime feature detection. Packed + * 16-bit integer arithmetic is used (SSE2's widest integer multiply is 16-bit; + * 32-bit packed multiply, pmulld, is an SSE4.1 extension and is avoided), plus + * saturating/averaging ops, widening multiply-add, shuffles and the packed + * single-precision float path including the reciprocal/rsqrt estimates. + * + * buf must be at least 128 bytes. Returns a checksum. + * =================================================================== */ +FN_BEGIN(fm_simd) + test rdi, rdi + jz .Lsd_zero + + /* seed eight vectors from the scratch buffer */ + movdqu xmm0, [rsi] + movdqu xmm1, [rsi + 16] + movdqu xmm2, [rsi + 32] + movdqu xmm3, [rsi + 48] + movdqu xmm4, [rsi + 64] + movdqu xmm5, [rsi + 80] + movdqu xmm6, [rsi + 96] + movdqu xmm7, [rsi + 112] + + /* float operands: convert the integer seeds to finite floats rather + * than reinterpreting bit patterns (which would be NaNs/denormals) */ + cvtdq2ps xmm12, xmm0 + cvtdq2ps xmm13, xmm1 + pxor xmm14, xmm14 /* float accumulator */ + +.Lsd_loop: + /* 16-bit integer multiply-accumulate across independent vectors */ + pmullw xmm0, xmm1 + paddw xmm0, xmm2 + pmullw xmm1, xmm2 + paddw xmm1, xmm3 + pmullw xmm2, xmm3 + paddw xmm2, xmm0 + + /* saturating and averaging arithmetic */ + paddsw xmm4, xmm0 + paddusw xmm5, xmm1 + psubsw xmm6, xmm2 + psubusw xmm7, xmm3 + + /* widening multiply-add: 16->32 bit lanes */ + movdqa xmm8, xmm0 + pmaddwd xmm8, xmm1 + paddd xmm3, xmm8 + + /* high-half multiply */ + movdqa xmm9, xmm0 + pmulhw xmm9, xmm1 + pxor xmm2, xmm9 + + /* shifts and logic */ + movdqa xmm10, xmm2 + pslld xmm10, 3 + psrld xmm2, 29 + por xmm2, xmm10 + pand xmm3, xmm4 + + /* min/max */ + pmaxsw xmm5, xmm0 + pminsw xmm6, xmm1 + + /* byte average and sum-of-absolute-differences */ + pavgb xmm7, xmm4 + movdqa xmm11, xmm0 + psadbw xmm11, xmm5 + paddw xmm4, xmm11 + + /* lane shuffle: reverse the four 32-bit lanes */ + pshufd xmm0, xmm0, 0x1B + pxor xmm1, xmm0 + + /* single-precision float SIMD: multiply-add plus the reciprocal and + * rsqrt estimates that shader-style code leans on */ + movaps xmm15, xmm12 + mulps xmm15, xmm13 + addps xmm14, xmm15 + rcpps xmm8, xmm12 + rsqrtps xmm9, xmm14 + addps xmm14, xmm8 + mulps xmm12, xmm9 + + dec rdi + jnz .Lsd_loop + + /* fold the eight integer vectors together */ + pxor xmm0, xmm1 + pxor xmm2, xmm3 + pxor xmm4, xmm5 + pxor xmm6, xmm7 + pxor xmm0, xmm2 + pxor xmm4, xmm6 + pxor xmm0, xmm4 + + /* fold in the float accumulator (truncate to int lanes) */ + cvttps2dq xmm14, xmm14 + pxor xmm0, xmm14 + + /* horizontal add of the four 32-bit lanes -> single checksum */ + pshufd xmm1, xmm0, 0x4E + paddd xmm0, xmm1 + pshufd xmm1, xmm0, 0xB1 + paddd xmm0, xmm1 + movd eax, xmm0 + ret + +.Lsd_zero: + xor eax, eax + ret +FN_END(fm_simd) + + +/* =================================================================== + * uint64_t fm_compress(const uint8_t *src, uint64_t len, uint32_t *ht) + * [rdi, rsi, rdx] + * + * The match-finding inner loop of an LZ77 compressor (the LZ4 fast strategy): + * hash the next 4 bytes, probe a single-entry-per-bucket table, verify, then + * extend. Branch-heavy with a data-dependent, cache-missing table probe. + * + * ht must hold 1<<16 uint32_t (256 KiB); this routine clears it itself. + * Returns the encoded size in bytes. + * =================================================================== */ +FN_BEGIN(fm_compress) + push rbp + push rbx + push r12 + push r13 + push r14 + push r15 + + mov r12, rdi /* src */ + mov r13, rdx /* ht */ + + /* clear the hash table: 1<<16 entries * 4 bytes = 262144 bytes */ + pxor xmm0, xmm0 + xor rax, rax + mov ecx, 262144 +.Lcm_clear: + movdqu [r13 + rax], xmm0 + movdqu [r13 + rax + 16], xmm0 + movdqu [r13 + rax + 32], xmm0 + movdqu [r13 + rax + 48], xmm0 + add rax, 64 + cmp rax, rcx + jb .Lcm_clear + + cmp rsi, 16 + jb .Lcm_tiny + + mov r14, r12 /* ip */ + mov r15, r12 /* anchor */ + lea rbx, [r12 + rsi] /* end */ + lea r10, [rbx - 12] /* mflimit = end - 12 */ + xor ebp, ebp /* outsize */ + +.Lcm_loop: + cmp r14, r10 + jae .Lcm_flush + + mov eax, [r14] /* seq = load32(ip) */ + imul eax, eax, 0x9E3779B1 /* * Knuth prime 2654435761 */ + shr eax, 16 /* h = (seq*prime) >> 16 */ + + mov ecx, [r13 + rax*4] /* ref_off = ht[h] */ + mov rdx, r14 + sub rdx, r12 /* cur_off = ip - src */ + mov [r13 + rax*4], edx /* ht[h] = cur_off */ + + lea rsi, [r12 + rcx] /* ref = src + ref_off */ + cmp rsi, r14 + jae .Lcm_no_match /* ref must be strictly behind ip */ + + mov rax, r14 + sub rax, rsi /* distance */ + cmp rax, 65536 + jae .Lcm_no_match /* 16-bit offset window */ + + mov eax, [rsi] + cmp eax, [r14] /* verify the 4-byte match */ + jne .Lcm_no_match + + /* match confirmed: extend it byte by byte */ + mov r9, 4 /* ml */ +.Lcm_extend: + lea rax, [r14 + r9] + cmp rax, rbx /* ip + ml vs end */ + jae .Lcm_emit + mov cl, [r14 + r9] + cmp cl, [rsi + r9] + jne .Lcm_emit + inc r9 + jmp .Lcm_extend + +.Lcm_emit: + /* token(1) + offset(2) + literals + varint extensions */ + mov rax, r14 + sub rax, r15 /* literal run length */ + add rbp, rax + add rbp, 3 + cmp rax, 15 + jb .Lcm_no_lit_ext + inc rbp /* literal-length extension byte */ +.Lcm_no_lit_ext: + cmp r9, 19 + jb .Lcm_no_ml_ext + inc rbp /* match-length extension byte */ +.Lcm_no_ml_ext: + add r14, r9 /* ip += ml */ + mov r15, r14 /* anchor = ip */ + jmp .Lcm_loop + +.Lcm_no_match: + inc r14 + jmp .Lcm_loop + +.Lcm_flush: + /* trailing literals */ + mov rax, rbx + sub rax, r15 /* end - anchor */ + add rbp, rax + add rbp, 1 + mov rax, rbp + jmp .Lcm_ret + +.Lcm_tiny: + lea rax, [rsi + 1] /* len + 1 */ + +.Lcm_ret: + pop r15 + pop r14 + pop r13 + pop r12 + pop rbx + pop rbp + ret +FN_END(fm_compress) + + +/* =================================================================== + * uint64_t fm_chacha20(uint8_t *buf, uint64_t len, + * const uint8_t key[32], uint64_t rounds) + * [rdi, rsi, rdx, rcx] + * + * ChaCha20 stream cipher, SSE2, four 128-bit state rows. Chosen over AES for + * the same reason the NEON file chose it: the AES-NI extension is optional, so + * an AES-instruction benchmark would fault (#UD) on cores without it. ChaCha20 + * needs only baseline SSE2 and is a real, widely deployed cipher. + * + * len is rounded down to a multiple of 64. `rounds` = passes over the buffer. + * Returns a checksum of the keystream output. + * =================================================================== */ + +/* rotate each 32-bit lane left by n, via shift-left + shift-right + or */ +#define ROL32(v, n) \ + movdqa xmm14, v ;\ + pslld v, n ;\ + psrld xmm14, (32 - (n)) ;\ + por v, xmm14 + +/* one ChaCha quarter-round over rows a,b,c,d (xmm14 is scratch, via ROL32) */ +#define QROUND(a, b, c, d) \ + paddd a, b ;\ + pxor d, a ;\ + ROL32(d, 16) ;\ + paddd c, d ;\ + pxor b, c ;\ + ROL32(b, 12) ;\ + paddd a, b ;\ + pxor d, a ;\ + ROL32(d, 8) ;\ + paddd c, d ;\ + pxor b, c ;\ + ROL32(b, 7) + +FN_BEGIN(fm_chacha20) + and rsi, -64 /* whole 64-byte blocks only */ + jz .Lcc_zero + test rcx, rcx + jz .Lcc_zero + + /* xmm4..7 hold the base state */ + movdqa xmm4, [rip + .Lcc_sigma] /* "expand 32-byte k" */ + movdqu xmm5, [rdx] /* key[0..15] */ + movdqu xmm6, [rdx + 16] /* key[16..31] */ + pxor xmm7, xmm7 /* counter || nonce = 0 */ + + pxor xmm13, xmm13 /* running checksum */ + xor r8, r8 /* block counter value */ + +.Lcc_pass: + xor r9, r9 /* byte offset into buf */ + +.Lcc_block: + /* working state = base state, block counter in lane 0 of row 3. + * Row 3 is all-zero (nonce and counter), so a plain movd both sets + * the counter lane and clears the nonce lanes. */ + movdqa xmm0, xmm4 + movdqa xmm1, xmm5 + movdqa xmm2, xmm6 + movd xmm3, r8d + + /* keep originals for the final feed-forward add */ + movdqa xmm8, xmm0 + movdqa xmm9, xmm1 + movdqa xmm10, xmm2 + movdqa xmm11, xmm3 + + mov r10d, 10 /* 10 double rounds = 20 rounds */ +.Lcc_rounds: + /* column round */ + QROUND(xmm0, xmm1, xmm2, xmm3) + + /* rotate lanes to form the diagonals */ + pshufd xmm1, xmm1, 0x39 /* <<< 1 lane */ + pshufd xmm2, xmm2, 0x4E /* <<< 2 lanes */ + pshufd xmm3, xmm3, 0x93 /* <<< 3 lanes */ + + /* diagonal round */ + QROUND(xmm0, xmm1, xmm2, xmm3) + + /* undo the lane rotation */ + pshufd xmm1, xmm1, 0x93 + pshufd xmm2, xmm2, 0x4E + pshufd xmm3, xmm3, 0x39 + + dec r10d + jnz .Lcc_rounds + + /* feed-forward: keystream = working + original */ + paddd xmm0, xmm8 + paddd xmm1, xmm9 + paddd xmm2, xmm10 + paddd xmm3, xmm11 + + /* XOR the keystream into the buffer */ + lea rax, [rdi + r9] + movdqu xmm12, [rax] + pxor xmm12, xmm0 + movdqu [rax], xmm12 + movdqu xmm12, [rax + 16] + pxor xmm12, xmm1 + movdqu [rax + 16], xmm12 + movdqu xmm12, [rax + 32] + pxor xmm12, xmm2 + movdqu [rax + 32], xmm12 + movdqu xmm12, [rax + 48] + pxor xmm12, xmm3 + movdqu [rax + 48], xmm12 + + /* accumulate a checksum of the keystream */ + pxor xmm13, xmm0 + pxor xmm13, xmm3 + + inc r8 /* counter++ */ + add r9, 64 + cmp r9, rsi + jb .Lcc_block + + dec rcx + jnz .Lcc_pass + + /* horizontal add of the checksum lanes */ + pshufd xmm0, xmm13, 0x4E + paddd xmm13, xmm0 + pshufd xmm0, xmm13, 0xB1 + paddd xmm13, xmm0 + movd eax, xmm13 + ret + +.Lcc_zero: + xor eax, eax + ret +FN_END(fm_chacha20) + + .p2align 4 +.Lcc_sigma: + .long 0x61707865, 0x3320646e, 0x79622d32, 0x6b206574 + + +/* =================================================================== + * uint64_t fm_physics(double *bodies, uint64_t n, uint64_t steps) + * [rdi, rsi, rdx] + * + * Direct-summation N-body gravity, O(n^2) per step, double precision. + * Layout per body, 8 doubles (64 bytes): [x y z mass vx vy vz pad]. + * The 1/sqrt is a real sqrtsd+divsd (not the rsqrt estimate), exercising the + * divide/sqrt unit the way physics code does. Returns a velocity checksum. + * =================================================================== */ +FN_BEGIN(fm_physics) + test rsi, rsi + jz .Lph_zero + test rdx, rdx + jz .Lph_zero + + movsd xmm13, [rip + .Lph_dt] /* dt */ + movsd xmm14, [rip + .Lph_eps2] /* eps^2 */ + movsd xmm15, [rip + .Lph_one] /* 1.0 */ + +.Lph_step: + xor r8, r8 /* i */ + +.Lph_body_i: + mov rax, r8 + shl rax, 6 /* i * 64 */ + lea r9, [rdi + rax] /* &bodies[i] */ + + movsd xmm0, [r9] /* xi */ + movsd xmm1, [r9 + 8] /* yi */ + movsd xmm2, [r9 + 16] /* zi */ + + xorpd xmm3, xmm3 /* ax */ + xorpd xmm4, xmm4 /* ay */ + xorpd xmm5, xmm5 /* az */ + + xor r10, r10 /* j */ + mov r11, rdi /* &bodies[j] */ + +.Lph_body_j: + movsd xmm6, [r11] /* xj */ + movsd xmm7, [r11 + 8] /* yj */ + movsd xmm8, [r11 + 16] /* zj */ + movsd xmm9, [r11 + 24] /* mj */ + + subsd xmm6, xmm0 /* dx */ + subsd xmm7, xmm1 /* dy */ + subsd xmm8, xmm2 /* dz */ + + /* r2 = dx*dx + dy*dy + dz*dz + eps^2 (>= eps^2, so the i==j self-term + * is finite and contributes exactly 0 below) */ + movsd xmm10, xmm6 + mulsd xmm10, xmm6 + movsd xmm11, xmm7 + mulsd xmm11, xmm7 + addsd xmm10, xmm11 + movsd xmm11, xmm8 + mulsd xmm11, xmm8 + addsd xmm10, xmm11 + addsd xmm10, xmm14 /* + eps^2 */ + + sqrtsd xmm10, xmm10 /* r */ + movsd xmm11, xmm15 + divsd xmm11, xmm10 /* 1/r */ + movsd xmm12, xmm11 + mulsd xmm12, xmm11 /* 1/r^2 */ + mulsd xmm12, xmm11 /* 1/r^3 */ + mulsd xmm12, xmm9 /* m/r^3 */ + + mulsd xmm6, xmm12 /* dx * m/r^3 */ + addsd xmm3, xmm6 + mulsd xmm7, xmm12 + addsd xmm4, xmm7 + mulsd xmm8, xmm12 + addsd xmm5, xmm8 + + add r11, 64 + inc r10 + cmp r10, rsi + jb .Lph_body_j + + /* v += a * dt */ + movsd xmm6, [r9 + 32] /* vx */ + movsd xmm7, [r9 + 40] /* vy */ + movsd xmm8, [r9 + 48] /* vz */ + mulsd xmm3, xmm13 + addsd xmm6, xmm3 + mulsd xmm4, xmm13 + addsd xmm7, xmm4 + mulsd xmm5, xmm13 + addsd xmm8, xmm5 + movsd [r9 + 32], xmm6 + movsd [r9 + 40], xmm7 + movsd [r9 + 48], xmm8 + + inc r8 + cmp r8, rsi + jb .Lph_body_i + + /* second pass: x += v * dt (positions move only after all forces) */ + xor r8, r8 + mov r11, rdi +.Lph_integrate: + movsd xmm0, [r11] + movsd xmm1, [r11 + 8] + movsd xmm2, [r11 + 16] + movsd xmm6, [r11 + 32] + movsd xmm7, [r11 + 40] + movsd xmm8, [r11 + 48] + mulsd xmm6, xmm13 + addsd xmm0, xmm6 + mulsd xmm7, xmm13 + addsd xmm1, xmm7 + mulsd xmm8, xmm13 + addsd xmm2, xmm8 + movsd [r11], xmm0 + movsd [r11 + 8], xmm1 + movsd [r11 + 16], xmm2 + add r11, 64 + inc r8 + cmp r8, rsi + jb .Lph_integrate + + dec rdx + jnz .Lph_step + + /* checksum: sum of all velocity components */ + xorpd xmm0, xmm0 + xor r8, r8 + mov r11, rdi +.Lph_sum: + movsd xmm6, [r11 + 32] + movsd xmm7, [r11 + 40] + movsd xmm8, [r11 + 48] + addsd xmm0, xmm6 + addsd xmm0, xmm7 + addsd xmm0, xmm8 + add r11, 64 + inc r8 + cmp r8, rsi + jb .Lph_sum + + movq rax, xmm0 + ret + +.Lph_zero: + xor eax, eax + ret +FN_END(fm_physics) + + .p2align 4 +.Lph_dt: + .double 0.0078125 /* dt */ +.Lph_eps2: + .double 0.0625 /* eps^2 */ +.Lph_one: + .double 1.0 + + +/* =================================================================== + * uint64_t fm_sort(uint32_t *a, uint64_t n) [rdi = a, rsi = n] + * + * In-place heapsort: no recursion or explicit stack, aggressively + * branch-unpredictable, with scattered memory access - it stresses the branch + * predictor and the cache hierarchy. Returns an order-sensitive checksum, + * which also verifies the sort. + * + * Uses only caller-saved registers, so no prologue is needed; the internal + * siftdown is reached with `call` (contract below). + * =================================================================== */ +FN_BEGIN(fm_sort) + cmp rsi, 2 + jb .Lst_trivial + + /* ---- build the max-heap: for i = n/2 - 1 down to 0 ---- */ + mov r8, rsi + shr r8, 1 /* i = n/2 */ +.Lst_build: + test r8, r8 + jz .Lst_extract + dec r8 /* i-- */ + mov rcx, r8 /* root */ + mov rdx, rsi /* end */ + call .Lst_siftdown + test r8, r8 + jnz .Lst_build + + /* ---- extract: for end = n-1 down to 1 ---- */ +.Lst_extract: + mov r9, rsi + dec r9 /* end = n-1 */ +.Lst_extract_loop: + test r9, r9 + jz .Lst_checksum + + /* swap a[0] and a[end] */ + mov eax, [rdi] + mov r10d, [rdi + r9*4] + mov [rdi], r10d + mov [rdi + r9*4], eax + + xor ecx, ecx /* root = 0 */ + mov rdx, r9 /* end = end */ + call .Lst_siftdown + + dec r9 + jmp .Lst_extract_loop + + /* ---- order-sensitive checksum ---- */ +.Lst_checksum: + xor eax, eax + xor rcx, rcx /* index */ +.Lst_cksum_loop: + mov r10d, [rdi + rcx*4] + xor rax, r10 + ror rax, 7 + add rax, r10 + inc rcx + cmp rcx, rsi + jb .Lst_cksum_loop + ret + +.Lst_trivial: + xor eax, eax + test rsi, rsi + jz .Lst_trivial_ret + mov eax, [rdi] +.Lst_trivial_ret: + ret + + /* ---- local helper: siftdown(root = rcx, end = rdx) + * base a = rdi; clobbers rax, rcx, r10, r11 only. r8 (build i), + * r9 (extract end), rsi (n), rdx (end) all survive. ---- */ +.Lst_siftdown: + mov r11, rcx /* root */ +.Lst_sift_loop: + lea r10, [r11 + r11 + 1] /* child = 2*root + 1 */ + cmp r10, rdx + jae .Lst_sift_done /* no children */ + + lea rax, [r10 + 1] /* right = child + 1 */ + cmp rax, rdx + jae .Lst_sift_have_child /* no right child */ + mov ecx, [rdi + rax*4] /* a[right] */ + cmp ecx, [rdi + r10*4] /* a[right] vs a[child] */ + jbe .Lst_sift_have_child /* keep child if a[right] <= it */ + mov r10, rax /* else child = right */ + +.Lst_sift_have_child: + mov eax, [rdi + r11*4] /* a[root] */ + mov ecx, [rdi + r10*4] /* a[child] */ + cmp eax, ecx + jae .Lst_sift_done /* heap property holds */ + + /* swap and descend */ + mov [rdi + r11*4], ecx + mov [rdi + r10*4], eax + mov r11, r10 + jmp .Lst_sift_loop + +.Lst_sift_done: + ret +FN_END(fm_sort) + + +/* =================================================================== + * uint64_t fm_chase(void **ptrs, uint64_t steps) [rdi = ptrs, rsi = steps] + * + * Pointer chase around a randomised cycle. Every load depends on the previous + * one, so nothing can be prefetched, overlapped or reordered - this measures + * the pure serial latency of the memory hierarchy. The truest single-threaded + * test in the suite. + * =================================================================== */ +FN_BEGIN(fm_chase) + test rsi, rsi + jz .Lch_zero + mov rax, rdi /* p = ptrs */ + mov rcx, rsi + +.Lch_loop: + mov rax, [rax] + dec rcx + jnz .Lch_loop + + sub rax, rdi /* final offset, keeps p live */ + ret + +.Lch_zero: + xor eax, eax + ret +FN_END(fm_chase) + + +#if defined(__ELF__) + .section .note.GNU-stack, "", @progbits +#endif diff --git a/src/main.c b/src/main.c new file mode 100644 index 0000000..ee8124c --- /dev/null +++ b/src/main.c @@ -0,0 +1,709 @@ +/* + * fossmark - a multi-core AArch64 CPU benchmark + * + * This file is the portable driver: it owns everything the assembly kernels + * deliberately do not (timing, memory, I/O, scoring). The kernels in + * fossmark.S are pure computation and identical on every OS; only this file + * knows what an operating system is. + * + * Every workload is run twice: once on a single core, and once on all available + * cores at once - one identical copy of the kernel per core, each with its own + * private buffers, so the machine is driven to 100%% and the rate is whole-machine + * throughput. From these two passes fossmark reports two composite scores, a + * SINGLECORE and a MULTICORE, from the same tests and the same weights. + * + * Build: cc -O2 -pthread main.c fossmark.S -o fossmark -lm + */ + +#include +#include +#include +#include +#include +#include +#include + +/* ---------- platform identification (for the banner only) ---------- */ + +#if defined(_WIN32) +# define FM_OS "Windows" +#elif defined(__APPLE__) +# define FM_OS "macOS" +#elif defined(__linux__) +# define FM_OS "Linux" +#else +# define FM_OS "POSIX" +#endif + +#if defined(__aarch64__) || defined(_M_ARM64) +# define FM_ARCH "ARM64" +# define D_INT "64-bit ALU: madd, umulh, udiv, bitops" +# define D_FP "double: fmadd, fdiv, fsqrt" +# define D_SIMD "NEON ASIMD: 128-bit integer + float" +#elif defined(__x86_64__) || defined(_M_X64) +# define FM_ARCH "x86-64" +# define D_INT "64-bit ALU: imul, mul, div, bitops" +# define D_FP "double: mulsd/addsd, divsd, sqrtsd" +# define D_SIMD "SSE2: 128-bit integer + float" +#else +# define FM_ARCH "unknown" +# define D_INT "64-bit integer ALU" +# define D_FP "double-precision FP" +# define D_SIMD "128-bit SIMD: integer + float" +#endif + +/* ---------- portable monotonic clock ---------- */ + +#if defined(_WIN32) +# define WIN32_LEAN_AND_MEAN +# include +static double now_seconds(void) +{ + LARGE_INTEGER f, t; + QueryPerformanceFrequency(&f); + QueryPerformanceCounter(&t); + return (double)t.QuadPart / (double)f.QuadPart; +} +#else +# include +static double now_seconds(void) +{ + struct timespec ts; + clock_gettime(CLOCK_MONOTONIC, &ts); + return (double)ts.tv_sec + (double)ts.tv_nsec * 1e-9; +} +#endif + +/* ---------- the assembly kernels ---------- */ + +extern uint64_t fm_int_math(uint64_t iters); +extern uint64_t fm_fp_math(uint64_t iters); +extern uint64_t fm_primes(uint64_t limit, uint8_t *sieve); +extern uint64_t fm_simd(uint64_t iters, void *buf); +extern uint64_t fm_compress(const uint8_t *src, uint64_t len, uint32_t *ht); +extern uint64_t fm_chacha20(uint8_t *buf, uint64_t len, + const uint8_t key[32], uint64_t rounds); +extern uint64_t fm_physics(double *bodies, uint64_t n, uint64_t steps); +extern uint64_t fm_sort(uint32_t *a, uint64_t n); +extern uint64_t fm_chase(void **ptrs, uint64_t steps); + +/* ---------- tuning ---------- */ + +#define PRIME_LIMIT (2u * 1000u * 1000u) /* sieve span */ +#define COMPRESS_LEN (4u * 1024u * 1024u) /* corpus size */ +#define HT_ENTRIES (1u << 16) /* LZ77 hash buckets */ +#define CIPHER_LEN (1u * 1024u * 1024u) /* plaintext size */ +#define SIMD_BUF 256 /* NEON scratch */ +#define NBODY_N 512 /* bodies */ +#define SORT_N (1u << 20) /* elements to sort */ +#define CHASE_NODES (1u << 21) /* 16 MiB cycle, > any L2 */ + +#define MIN_SECONDS 2.0 /* per-test measured floor */ +#define REPEATS 3 /* best-of, to reject noise */ + +/* ---------- scoring configuration ---------- + * + * The overall score is a WEIGHTED geometric mean of each test's rate expressed + * relative to a reference machine. Two knobs per test: + * + * FM_REF_* the reference rate (this machine's measured rate). A machine + * matching the reference scores FM_TARGET_SCORE on that test. + * FM_WEIGHT_* how much that test counts toward the overall, by its + * influence on everyday user experience. Weights are relative: + * only their ratios matter, so they need not sum to anything - + * the code normalises by their sum. (They happen to sum to 100 + * here, so each reads as a percent.) + * + * Per-test score: S_i = FM_TARGET_SCORE * (rate_i / FM_REF_i) + * Overall score: Overall = FM_TARGET_SCORE * + * exp( Sum(w_i * ln(rate_i/FM_REF_i)) / Sum(w_i) ) + * + * On the reference machine every ratio is 1, so every S_i and the overall come + * out to exactly FM_TARGET_SCORE, regardless of the weights. Scaling is linear + * in performance, so far slower machines fall well below (half as fast -> half + * the score) and faster future machines rise above. + */ + +#define FM_TARGET_SCORE 10000.0 /* reference-machine overall */ + +/* Reference rates: this machine, in each test's native unit (see tests[]). */ +#define FM_REF_INT 3086.0 /* Mops/s */ +#define FM_REF_FP 1682.0 /* Mops/s */ +#define FM_REF_PRIMES 812.0 /* Mcand/s */ +#define FM_REF_SIMD 6576.0 /* Mops/s */ +#define FM_REF_COMPRESS 674.0 /* MB/s */ +#define FM_REF_CRYPTO 406.0 /* MB/s */ +#define FM_REF_PHYSICS 631.0 /* Mpair/s */ +#define FM_REF_SORT 363.0 /* Mkey-cmp/s*/ +#define FM_REF_CHASE 79.0 /* Mhop/s (scoring); shown as ns/access */ + +/* Weights: influence on day-to-day, common-workload user experience. + * Rationale: integer/general-purpose code and memory-latency-bound + * responsiveness dominate everyday use; specialised FP/physics matter least. + * Roughly an 80/20 integer-vs-FP split, in the spirit of Geekbench 6's + * weighted, integer-dominant methodology. Retune freely. */ +#define FM_WEIGHT_INT 20.0 /* general-purpose ALU: everything */ +#define FM_WEIGHT_CHASE 16.0 /* memory latency: responsiveness */ +#define FM_WEIGHT_COMPRESS 14.0 /* web, storage, RAM compression */ +#define FM_WEIGHT_SORT 12.0 /* general data-structure work */ +#define FM_WEIGHT_SIMD 11.0 /* codecs, mem/string ops, parsing */ +#define FM_WEIGHT_FP 9.0 /* spreadsheets, app/media math */ +#define FM_WEIGHT_CRYPTO 8.0 /* TLS, disk encryption (small frac) */ +#define FM_WEIGHT_PRIMES 6.0 /* synthetic ALU+memory proxy */ +#define FM_WEIGHT_PHYSICS 4.0 /* niche simulation/games */ + +/* ---------- deterministic PRNG (splitmix64) ---------- */ + +static uint64_t rng_state = 0x853c49e6748fea9bULL; + +static uint64_t rng_next(void) +{ + uint64_t z = (rng_state += 0x9e3779b97f4a7c15ULL); + z = (z ^ (z >> 30)) * 0xbf58476d1ce4e5b9ULL; + z = (z ^ (z >> 27)) * 0x94d049bb133111ebULL; + return z ^ (z >> 31); +} + +static void rng_reset(void) { rng_state = 0x853c49e6748fea9bULL; } + +/* ---------- aligned allocation ---------- */ + +static void *xalloc(size_t n) +{ + void *p = NULL; +#if defined(_WIN32) + p = _aligned_malloc(n, 64); +#else + if (posix_memalign(&p, 64, n) != 0) + p = NULL; +#endif + if (!p) { + fprintf(stderr, "fossmark: out of memory (%zu bytes)\n", n); + exit(1); + } + return p; +} + +static void xfree(void *p) +{ +#if defined(_WIN32) + _aligned_free(p); +#else + free(p); +#endif +} + +/* ---------- workload state ---------- + * + * Because every core runs the same kernel simultaneously, each core needs its + * OWN mutable buffers - sharing them would be a data race and would corrupt + * both the results and the determinism check. Per-core scratch lives in a + * `workspace`, one per thread. Read-only inputs (the corpus, the pristine + * physics/sort seeds, the key, the chase graph) are genuinely shared. + */ +struct workspace { + uint8_t *sieve; /* prime sieve scratch */ + uint32_t *ht; /* LZ77 hash table scratch */ + uint8_t *cipher_buf; /* ChaCha20 buffer, encrypted in place */ + uint8_t *simd_buf; /* NEON scratch */ + double *bodies; /* n-body integration buffer */ + uint32_t *sort_work; /* the buffer we actually sort */ + void **chase; /* private 16 MiB pointer-chase cycle */ +}; + +static long g_ncores = 1; /* active online cores */ +static struct workspace *g_ws; /* g_ncores per-thread workspaces */ + +static uint8_t *g_corpus; /* shared, read-only compression input */ +static uint8_t g_key[32]; /* shared, read-only cipher key */ +static uint8_t *g_cipher_src; /* pristine plaintext, copied per-core */ +static uint8_t *g_simd_src; /* pristine NEON seed, copied per-core */ +static double *g_bodies_src; /* pristine initial conditions */ +static uint32_t *g_sort_src; /* pristine unsorted data */ + +/* + * Synthesise a compressible corpus. Random bytes would be incompressible and + * would make the match-finder trivially miss every probe, measuring nothing + * interesting. This builds text-like data with realistic repetition instead. + */ +static void build_corpus(uint8_t *buf, size_t len) +{ + static const char *words[] = { + "the", "quick", "brown", "fox", "jumps", "over", "lazy", + "dog", "benchmark", "processor", "assembly", "vector", + "memory", "cache", "pipeline", "instruction", "compress", + "data", "system", "performance", "register", "kernel" + }; + const size_t nwords = sizeof(words) / sizeof(words[0]); + size_t pos = 0; + + while (pos < len) { + const char *w = words[rng_next() % nwords]; + size_t wl = strlen(w); + + if (pos + wl + 1 > len) + break; + memcpy(buf + pos, w, wl); + pos += wl; + buf[pos++] = (rng_next() % 8 == 0) ? '\n' : ' '; + } + while (pos < len) + buf[pos++] = ' '; +} + +/* Build a single random cycle through the node array (Sattolo's algorithm), + * guaranteeing one cycle of exactly CHASE_NODES steps with no early closure. */ +static void build_chase(void **nodes, size_t n) +{ + size_t *perm = xalloc(n * sizeof(size_t)); + size_t i; + + for (i = 0; i < n; i++) + perm[i] = i; + for (i = n - 1; i > 0; i--) { + size_t j = (size_t)(rng_next() % i); /* strictly j < i */ + size_t t = perm[i]; + perm[i] = perm[j]; + perm[j] = t; + } + for (i = 0; i < n; i++) + nodes[perm[i]] = (void *)&nodes[perm[(i + 1) % n]]; + + xfree(perm); +} + +static void setup(void) +{ + size_t i; + long t; + + rng_reset(); + + /* shared read-only inputs and pristine per-core seeds */ + g_corpus = xalloc(COMPRESS_LEN); + g_cipher_src = xalloc(CIPHER_LEN); + g_simd_src = xalloc(SIMD_BUF); + g_bodies_src = xalloc(NBODY_N * 8 * sizeof(double)); + g_sort_src = xalloc(SORT_N * sizeof(uint32_t)); + + build_corpus(g_corpus, COMPRESS_LEN); + + for (i = 0; i < CIPHER_LEN; i++) + g_cipher_src[i] = (uint8_t)rng_next(); + for (i = 0; i < 32; i++) + g_key[i] = (uint8_t)rng_next(); + for (i = 0; i < SIMD_BUF; i++) + g_simd_src[i] = (uint8_t)rng_next(); + for (i = 0; i < SORT_N; i++) + g_sort_src[i] = (uint32_t)rng_next(); + + /* bodies: [x y z mass vx vy vz pad], positions in a unit-ish cube */ + for (i = 0; i < NBODY_N; i++) { + double *b = &g_bodies_src[i * 8]; + b[0] = (double)(rng_next() % 2000) / 1000.0 - 1.0; + b[1] = (double)(rng_next() % 2000) / 1000.0 - 1.0; + b[2] = (double)(rng_next() % 2000) / 1000.0 - 1.0; + b[3] = (double)(rng_next() % 900) / 1000.0 + 0.1; /* mass > 0 */ + b[4] = b[5] = b[6] = 0.0; + b[7] = 0.0; + } + + /* one private workspace per core, every copy seeded identically so all + * cores compute the same deterministic result. Each core also gets its + * own pointer-chase cycle: sharing one would collapse the multi-core + * latency test into a shared-cache test instead of a memory test. */ + g_ws = xalloc((size_t)g_ncores * sizeof *g_ws); + for (t = 0; t < g_ncores; t++) { + struct workspace *w = &g_ws[t]; + + w->sieve = xalloc(PRIME_LIMIT); + w->ht = xalloc(HT_ENTRIES * sizeof(uint32_t)); + w->cipher_buf = xalloc(CIPHER_LEN); + w->simd_buf = xalloc(SIMD_BUF); + w->bodies = xalloc(NBODY_N * 8 * sizeof(double)); + w->sort_work = xalloc(SORT_N * sizeof(uint32_t)); + w->chase = xalloc(CHASE_NODES * sizeof(void *)); + + memcpy(w->cipher_buf, g_cipher_src, CIPHER_LEN); + memcpy(w->simd_buf, g_simd_src, SIMD_BUF); + build_chase(w->chase, CHASE_NODES); + } +} + +static void teardown(void) +{ + long t; + + for (t = 0; t < g_ncores; t++) { + struct workspace *w = &g_ws[t]; + + xfree(w->sieve); xfree(w->ht); xfree(w->cipher_buf); + xfree(w->simd_buf); xfree(w->bodies); xfree(w->sort_work); + xfree(w->chase); + } + xfree(g_ws); + + xfree(g_corpus); xfree(g_cipher_src); xfree(g_simd_src); + xfree(g_bodies_src); xfree(g_sort_src); +} + +/* ---------- the test harness ---------- */ + +/* + * Each test runs a kernel `n` times against a per-core workspace and returns a + * checksum. The harness auto-calibrates `n` upward until the run exceeds + * MIN_SECONDS, so the result is insensitive to clock granularity and to how + * fast the machine is. + */ +typedef uint64_t (*run_fn)(uint64_t n, struct workspace *ws); + +struct test { + const char *name; + const char *detail; + run_fn run; + uint64_t start_n; + double work_per_n; /* abstract work units, for scoring */ + const char *unit; + double ref_rate; /* reference-machine rate, in `unit` */ + double weight; /* relative weight in the overall score */ +}; + +static uint64_t run_int(uint64_t n, struct workspace *ws) +{ + (void)ws; + return fm_int_math(n * 100000); +} +static uint64_t run_fp(uint64_t n, struct workspace *ws) +{ + (void)ws; + return fm_fp_math(n * 100000); +} +static uint64_t run_primes(uint64_t n, struct workspace *ws) +{ + uint64_t c = 0; + for (uint64_t i = 0; i < n; i++) + c += fm_primes(PRIME_LIMIT, ws->sieve); + return c; +} +static uint64_t run_simd(uint64_t n, struct workspace *ws) +{ + return fm_simd(n * 100000, ws->simd_buf); +} +static uint64_t run_compress(uint64_t n, struct workspace *ws) +{ + uint64_t c = 0; + for (uint64_t i = 0; i < n; i++) + c += fm_compress(g_corpus, COMPRESS_LEN, ws->ht); + return c; +} +static uint64_t run_crypto(uint64_t n, struct workspace *ws) +{ + return fm_chacha20(ws->cipher_buf, CIPHER_LEN, g_key, n); +} +static uint64_t run_physics(uint64_t n, struct workspace *ws) +{ + /* restore initial conditions: the integrator mutates the bodies, so + * a re-run must start from the same state to be reproducible */ + memcpy(ws->bodies, g_bodies_src, NBODY_N * 8 * sizeof(double)); + return fm_physics(ws->bodies, NBODY_N, n); +} +static uint64_t run_sort(uint64_t n, struct workspace *ws) +{ + uint64_t c = 0; + for (uint64_t i = 0; i < n; i++) { + /* restore the pristine data: sorting an already-sorted array + * would measure the best case, not the real one */ + memcpy(ws->sort_work, g_sort_src, SORT_N * sizeof(uint32_t)); + c ^= fm_sort(ws->sort_work, SORT_N); + } + return c; +} +static uint64_t run_chase(uint64_t n, struct workspace *ws) +{ + return fm_chase(ws->chase, n * 1000000); +} + +static const struct test tests[] = { + { "Integer Math", D_INT, + run_int, 20, 100000.0 * 24, "Mops/s", + FM_REF_INT, FM_WEIGHT_INT }, + { "Floating Point Math", D_FP, + run_fp, 20, 100000.0 * 20, "Mops/s", + FM_REF_FP, FM_WEIGHT_FP }, + { "Prime Numbers", "sieve of Eratosthenes to 2M", + run_primes, 1, (double)PRIME_LIMIT, "Mcand/s", + FM_REF_PRIMES, FM_WEIGHT_PRIMES }, + { "Extended Instructions",D_SIMD, + run_simd, 10, 100000.0 * 32, "Mops/s", + FM_REF_SIMD, FM_WEIGHT_SIMD }, + { "Compression", "LZ77 match finder, 4 MiB corpus", + run_compress, 1, (double)COMPRESS_LEN, "MB/s", + FM_REF_COMPRESS, FM_WEIGHT_COMPRESS }, + { "Encryption", "ChaCha20, 20 rounds, 1 MiB", + run_crypto, 4, (double)CIPHER_LEN, "MB/s", + FM_REF_CRYPTO, FM_WEIGHT_CRYPTO }, + { "Physics", "512-body direct-sum gravity", + run_physics, 4, (double)NBODY_N * NBODY_N, "Mpair/s", + FM_REF_PHYSICS, FM_WEIGHT_PHYSICS }, + { "Sorting", "heapsort, 1M uint32", + run_sort, 1, (double)SORT_N * 20, "Mkey-cmp/s", + FM_REF_SORT, FM_WEIGHT_SORT }, + { "Memory Latency", "dependent-load pointer chase, 16 MiB", + run_chase, 1, 1000000.0, "ns/access", + FM_REF_CHASE, FM_WEIGHT_CHASE }, +}; + +#define NTESTS (sizeof(tests) / sizeof(tests[0])) + +struct result { + double rate; /* work units per second */ + double score; + uint64_t checksum; + double seconds; + uint64_t iters; + int threads; /* cores this test was spread across */ +}; + +/* + * One unit of parallel work: run `run(n, ws)` on a private workspace. Every + * core executes the identical kernel on identically-seeded data, so all cores + * return the same checksum; the harness sums them into one aggregate that stays + * deterministic across repeats. + */ +struct job { + run_fn run; + uint64_t n; + struct workspace *ws; + uint64_t result; +}; + +static void *job_entry(void *arg) +{ + struct job *j = arg; + j->result = j->run(j->n, j->ws); + return NULL; +} + +/* + * Run the kernel on `threads` cores at once and return the summed checksum. + * The calling thread runs job 0 itself; threads 1..N-1 run on spawned workers. + * A thread that fails to spawn simply runs inline, so the benchmark still + * completes (with less parallelism) rather than aborting. + */ +static uint64_t dispatch(run_fn run, uint64_t n, int threads) +{ + struct job *jobs = xalloc((size_t)threads * sizeof *jobs); + pthread_t *tids = threads > 1 + ? xalloc((size_t)(threads - 1) * sizeof *tids) : NULL; + int i, spawned = 0; + uint64_t agg = 0; + + for (i = 0; i < threads; i++) { + jobs[i].run = run; + jobs[i].n = n; + jobs[i].ws = &g_ws[i]; + } + for (i = 1; i < threads; i++) { + if (pthread_create(&tids[spawned], NULL, job_entry, &jobs[i]) == 0) + spawned++; + else + job_entry(&jobs[i]); /* fall back to inline */ + } + + job_entry(&jobs[0]); /* this thread runs job 0 */ + + for (i = 0; i < spawned; i++) + pthread_join(tids[i], NULL); + for (i = 0; i < threads; i++) + agg += jobs[i].result; + + xfree(jobs); + xfree(tids); + return agg; +} + +static struct result run_test(const struct test *t, int threads) +{ + struct result r; + uint64_t n = t->start_n; + double elapsed = 0.0, best = 0.0; + uint64_t checksum = 0; + int i; + + /* calibrate: grow n until a single run clears the noise floor */ + for (;;) { + double t0 = now_seconds(); + checksum = dispatch(t->run, n, threads); + elapsed = now_seconds() - t0; + + if (elapsed >= MIN_SECONDS) + break; + if (elapsed < 0.001) { + n *= 8; /* far too fast to measure */ + } else { + double scale = (MIN_SECONDS * 1.3) / elapsed; + if (scale < 1.5) + scale = 1.5; + if (scale > 8.0) + scale = 8.0; + n = (uint64_t)((double)n * scale) + 1; + } + } + + /* best-of: the fastest run is the one least disturbed by the OS */ + best = elapsed; + for (i = 1; i < REPEATS; i++) { + double t0 = now_seconds(); + uint64_t c = dispatch(t->run, n, threads); + double e = now_seconds() - t0; + + if (c != checksum) { + fprintf(stderr, + "fossmark: %s is non-deterministic " + "(checksum %llu != %llu)\n", t->name, + (unsigned long long)c, + (unsigned long long)checksum); + exit(2); + } + if (e < best) + best = e; + } + + r.seconds = best; + r.iters = n; + r.checksum = checksum; + r.threads = threads; + /* aggregate throughput: `threads` cores each did n*work_per_n of work in + * the same wall-clock window, so the machine's rate is their sum */ + r.rate = ((double)threads * (double)n * t->work_per_n) / best / 1e6; + /* normalise against the reference machine: this is the per-test score */ + r.score = FM_TARGET_SCORE * (r.rate / t->ref_rate); + return r; +} + +/* + * The number shown in the RATE column. Most tests report throughput in their + * `unit`. The memory-latency test is different: throughput (hops/s) is not what + * anyone reasons about for memory, so we report the actual per-access latency + * in nanoseconds instead. That is a per-core property -- the time for one + * dependent load in the chain -- so it is derived from a single core's hop + * count and is independent of how many cores ran, unlike the aggregate `rate`. + */ +static double display_metric(const struct test *t, const struct result *r) +{ + if (t->run == run_chase) { + double hops = (double)r->iters * t->work_per_n; /* per core */ + return r->seconds / hops * 1e9; /* ns/access */ + } + return r->rate; +} + +/* ---------- output ---------- */ + +static void print_header(void) +{ + printf("\n"); + printf(" fossmark 1.0 - multi-core CPU benchmark\n"); + printf(" ------------------------------------------------------------------\n"); + printf(" platform: %s/%s\n", FM_OS, FM_ARCH); + printf(" cores: %ld (each test is run once on 1 core, once on all %ld)\n", + g_ncores, g_ncores); + printf("\n"); + printf(" RATE/TIME/SCORE below are the all-core (multi-core) pass.\n"); + printf("\n"); + printf(" %-24s %12s %-11s %8s %9s\n", + "TEST", "RATE", "UNIT", "TIME", "SCORE"); + printf(" --------------------------------------------------------------------------\n"); + fflush(stdout); +} + +int main(int argc, char **argv) +{ + struct result multi[NTESTS], single[NTESTS]; + double multi_log_sum = 0.0, single_log_sum = 0.0; + double weight_sum = 0.0; + int verbose = 0; + size_t i; + + for (i = 1; i < (size_t)argc; i++) { + if (strcmp(argv[i], "-v") == 0 || + strcmp(argv[i], "--verbose") == 0) { + verbose = 1; + } else if (strcmp(argv[i], "-h") == 0 || + strcmp(argv[i], "--help") == 0) { + printf("usage: %s [-v|--verbose]\n", argv[0]); + return 0; + } else { + fprintf(stderr, "fossmark: unknown option '%s'\n", + argv[i]); + return 1; + } + } + + { + long n = sysconf(_SC_NPROCESSORS_ONLN); + g_ncores = n > 0 ? n : 1; + } + + printf("\n preparing workloads..."); + fflush(stdout); + setup(); + printf(" done\n"); + + print_header(); + + for (i = 0; i < NTESTS; i++) { + double sm, ss; + + printf(" %-24s", tests[i].name); + fflush(stdout); + + /* each test runs twice: the all-core pass (shown) and the + * single-core pass (folded into the SINGLECORE score) */ + multi[i] = run_test(&tests[i], (int)g_ncores); + single[i] = run_test(&tests[i], 1); + + printf(" %12.1f %-11s %7.2fs %9.0f\n", + display_metric(&tests[i], &multi[i]), tests[i].unit, + multi[i].seconds, multi[i].score); + if (verbose) + printf(" %-24s %s\n" + " %-24s weight=%.0f%% 1-core: %.1f %s / %.0f %ld-core: %.1f %s / %.0f\n", + "", tests[i].detail, "", tests[i].weight, + display_metric(&tests[i], &single[i]), tests[i].unit, + single[i].score, g_ncores, + display_metric(&tests[i], &multi[i]), tests[i].unit, + multi[i].score); + fflush(stdout); + + /* accumulate the weighted geometric mean of BOTH passes, same + * weights, so the two composite scores are directly comparable */ + sm = multi[i].score > 0.0 ? multi[i].score : 1e-9; + ss = single[i].score > 0.0 ? single[i].score : 1e-9; + multi_log_sum += tests[i].weight * log(sm); + single_log_sum += tests[i].weight * log(ss); + weight_sum += tests[i].weight; + } + + printf(" --------------------------------------------------------------------------\n"); + + /* + * Two composite scores, each the WEIGHTED geometric mean of the per-test + * scores from one pass. Per-test scores are already normalised so the + * single-thread reference machine reads FM_TARGET_SCORE. Geometric rather + * than arithmetic so no single test dominates; weighted so tests count in + * proportion to their influence on everyday use (the FM_WEIGHT_* config). + * The two passes share tests and weights, so MULTICORE / SINGLECORE is a + * clean read of how much the machine gains from all its cores. + */ + printf(" %-24s %44.0f\n", "MULTICORE SCORE", + exp(multi_log_sum / weight_sum)); + printf(" %-24s %44.0f\n", "SINGLECORE SCORE", + exp(single_log_sum / weight_sum)); + printf(" %-24s (weighted geometric means, reference machine = %.0f)\n", + "", FM_TARGET_SCORE); + printf("\n"); + + teardown(); + return 0; +} diff --git a/src/test_kernels.c b/src/test_kernels.c new file mode 100644 index 0000000..9ae3595 --- /dev/null +++ b/src/test_kernels.c @@ -0,0 +1,462 @@ +/* + * test_kernels.c - correctness checks for the fossmark assembly kernels + * + * The benchmark's own best-of-N run guards against non-determinism, but a + * kernel can be perfectly deterministic and still wrong. This file is the + * "single C file to poke at and test with": it validates each kernel against + * an independent reference or an invariant, so a mistake in the assembly is + * caught here rather than silently skewing a score. + * + * Every check (except the single-threaded pointer-chase) is run concurrently + * on all available cores. The kernels take their buffers as arguments and hold + * no shared state, so a correct kernel must give identical, correct results no + * matter how many copies run at once; a hidden global or a reentrancy bug would + * survive a single-threaded run but fail here. + * + * Build: cc -O2 -pthread test_kernels.c fossmark.S -o test_kernels -lm + * Exit status is 0 iff every check passes. + */ + +#include +#include +#include +#include +#include +#include +#include +#include + +extern uint64_t fm_int_math(uint64_t iters); +extern uint64_t fm_fp_math(uint64_t iters); +extern uint64_t fm_primes(uint64_t limit, uint8_t *sieve); +extern uint64_t fm_simd(uint64_t iters, void *buf); +extern uint64_t fm_compress(const uint8_t *src, uint64_t len, uint32_t *ht); +extern uint64_t fm_chacha20(uint8_t *buf, uint64_t len, + const uint8_t key[32], uint64_t rounds); +extern uint64_t fm_physics(double *bodies, uint64_t n, uint64_t steps); +extern uint64_t fm_sort(uint32_t *a, uint64_t n); +extern uint64_t fm_chase(void **ptrs, uint64_t steps); + +static int failures = 0; +static int checks = 0; + +/* + * Concurrency plumbing. Each check runs on every core at once; the counters and + * stdout are shared, so ok()/note() serialise on this lock. `fm_primary` is set + * on exactly one thread per check (the one running on the main thread): it owns + * the human-readable output so the "[ ok ]" lines and diagnostics appear once, + * not once per core. Every thread still evaluates every assertion, so a failure + * on any core - even a silent secondary - is reported and counted. + */ +static pthread_mutex_t io_lock = PTHREAD_MUTEX_INITIALIZER; +static __thread int fm_primary = 1; +static long fm_ncores = 1; + +static void ok(const char *what, int cond) +{ + pthread_mutex_lock(&io_lock); + if (fm_primary) { + checks++; + if (cond) { + printf(" [ ok ] %s\n", what); + } else { + printf(" [FAIL] %s\n", what); + failures++; + } + } else if (!cond) { + /* a secondary core disagrees: surface it explicitly */ + printf(" [FAIL] %s (concurrent core)\n", what); + failures++; + } + pthread_mutex_unlock(&io_lock); +} + +/* Diagnostic output that should appear once per check, not once per core. */ +static void note(const char *fmt, ...) +{ + va_list ap; + + if (!fm_primary) + return; + pthread_mutex_lock(&io_lock); + va_start(ap, fmt); + vprintf(fmt, ap); + va_end(ap); + pthread_mutex_unlock(&io_lock); +} + +/* Run `check` on every core simultaneously. The main thread is the primary; + * fm_ncores-1 workers run the same check as silent secondaries. */ +static void *fm_worker(void *arg) +{ + void (*check)(void) = *(void (**)(void))arg; + + fm_primary = 0; + check(); + return NULL; +} + +static void parallel(void (*check)(void)) +{ + long extra = fm_ncores - 1; + pthread_t *th = NULL; + long i, spawned = 0; + + if (extra > 0) { + th = calloc((size_t)extra, sizeof *th); + if (th) { + for (i = 0; i < extra; i++) + if (pthread_create(&th[spawned], NULL, + fm_worker, &check) == 0) + spawned++; + } + } + + check(); /* primary runs on this thread */ + + for (i = 0; i < spawned; i++) + pthread_join(th[i], NULL); + free(th); +} + +/* ---------- reference implementations ---------- */ + +static uint64_t ref_prime_count(uint64_t limit) +{ + uint8_t *s = calloc(limit, 1); + uint64_t count = 0, i, j; + + for (i = 2; i * i < limit; i++) + if (!s[i]) + for (j = i * i; j < limit; j += i) + s[j] = 1; + for (i = 2; i < limit; i++) + if (!s[i]) + count++; + free(s); + return count; +} + +/* A textbook scalar ChaCha20 block function, used both to anchor against the + * RFC 8439 known-answer vector and to validate the NEON kernel block-for-block. + * `out` receives 64 keystream bytes for the given counter and 12-byte nonce. */ +#define ROTL32(x, n) (((x) << (n)) | ((x) >> (32 - (n)))) + +static void ref_chacha_block(uint32_t out_words[16], const uint8_t key[32], + uint32_t counter, const uint8_t nonce[12]) +{ + static const uint32_t c[4] = { + 0x61707865, 0x3320646e, 0x79622d32, 0x6b206574 + }; + uint32_t s[16], x[16]; + int i; + + for (i = 0; i < 4; i++) + s[i] = c[i]; + for (i = 0; i < 8; i++) + s[4 + i] = (uint32_t)key[4 * i] | (uint32_t)key[4 * i + 1] << 8 | + (uint32_t)key[4 * i + 2] << 16 | + (uint32_t)key[4 * i + 3] << 24; + s[12] = counter; + for (i = 0; i < 3; i++) + s[13 + i] = (uint32_t)nonce[4 * i] | (uint32_t)nonce[4 * i + 1] << 8 | + (uint32_t)nonce[4 * i + 2] << 16 | + (uint32_t)nonce[4 * i + 3] << 24; + + memcpy(x, s, sizeof x); +#define QR(a, b, cc, d) \ + x[a] += x[b]; x[d] ^= x[a]; x[d] = ROTL32(x[d], 16); \ + x[cc] += x[d]; x[b] ^= x[cc]; x[b] = ROTL32(x[b], 12); \ + x[a] += x[b]; x[d] ^= x[a]; x[d] = ROTL32(x[d], 8); \ + x[cc] += x[d]; x[b] ^= x[cc]; x[b] = ROTL32(x[b], 7) + for (i = 0; i < 10; i++) { + QR(0, 4, 8, 12); QR(1, 5, 9, 13); + QR(2, 6, 10, 14); QR(3, 7, 11, 15); + QR(0, 5, 10, 15); QR(1, 6, 11, 12); + QR(2, 7, 8, 13); QR(3, 4, 9, 14); + } +#undef QR + for (i = 0; i < 16; i++) + out_words[i] = x[i] + s[i]; +} + +/* ---------- checks ---------- */ + +static void check_int(void) +{ + /* determinism and non-triviality: the checksum must be stable and + * must actually change with the iteration count */ + uint64_t a = fm_int_math(1000); + uint64_t b = fm_int_math(1000); + uint64_t c = fm_int_math(2000); + + ok("int_math is deterministic", a == b); + ok("int_math depends on iters", a != c); + ok("int_math(0) is zero", fm_int_math(0) == 0); +} + +static void check_fp(void) +{ + uint64_t a = fm_fp_math(1000); + uint64_t b = fm_fp_math(1000); + double da; + + memcpy(&da, &a, sizeof da); + ok("fp_math is deterministic", a == b); + ok("fp_math result is finite", isfinite(da)); + ok("fp_math(0) is zero", fm_fp_math(0) == 0); +} + +static void check_primes(void) +{ + enum { LIM = 1000000 }; + uint8_t *sieve = malloc(LIM); + uint64_t got = fm_primes(LIM, sieve); + uint64_t ref = ref_prime_count(LIM); + + note(" primes < %d: got %llu, expected %llu\n", + LIM, (unsigned long long)got, (unsigned long long)ref); + ok("primes matches reference sieve", got == ref); + ok("primes < 10 == 4", fm_primes(10, sieve) == 4); /* 2,3,5,7 */ + ok("primes < 2 == 0", fm_primes(2, sieve) == 0); + free(sieve); +} + +static void check_simd(void) +{ + uint8_t *buf = aligned_alloc(16, 256); + uint64_t a, b; + + memset(buf, 0xA5, 256); + a = fm_simd(500, buf); + memset(buf, 0xA5, 256); + b = fm_simd(500, buf); + ok("simd is deterministic", a == b); + ok("simd(0) is zero", fm_simd(0, buf) == 0); + free(buf); +} + +static void check_compress(void) +{ + enum { N = 65536 }; + uint8_t *src = malloc(N); + uint32_t *ht = malloc((1 << 16) * sizeof(uint32_t)); + uint64_t incompressible, compressible; + size_t i; + + /* genuinely incompressible data (splitmix64 output): with no matches + * to exploit, an LZ coder's output must be at least the input size */ + { + uint64_t st = 0x1234567890abcdefULL; + for (i = 0; i < N; i++) { + uint64_t z = (st += 0x9e3779b97f4a7c15ULL); + z = (z ^ (z >> 30)) * 0xbf58476d1ce4e5b9ULL; + z = (z ^ (z >> 27)) * 0x94d049bb133111ebULL; + src[i] = (uint8_t)(z ^ (z >> 31)); + } + } + incompressible = fm_compress(src, N, ht); + + /* all-zero data is maximally compressible: it must shrink hugely */ + memset(src, 0, N); + compressible = fm_compress(src, N, ht); + + note(" 64KiB random -> %llu bytes, 64KiB zeros -> %llu bytes\n", + (unsigned long long)incompressible, + (unsigned long long)compressible); + ok("compress expands random data", incompressible >= N); + ok("compress shrinks constant data", compressible < N / 10); + ok("compress is deterministic", fm_compress(src, N, ht) == compressible); + free(src); + free(ht); +} + +static void check_crypto(void) +{ + uint8_t key[32]; + size_t i; + + /* (1) anchor the scalar reference to the RFC 8439 s.2.3.2 vector: + * key = 00,01,...,1f; counter = 1; nonce = 00,00,00,09,...,4a,... + * serialised keystream begins 10 f1 e7 e4. */ + { + uint32_t w[16]; + uint8_t rnonce[12] = {0,0,0,9, 0,0,0,0x4a, 0,0,0,0}; + uint8_t ks0[4]; + for (i = 0; i < 32; i++) + key[i] = (uint8_t)i; + ref_chacha_block(w, key, 1, rnonce); + for (i = 0; i < 4; i++) + ks0[i] = (uint8_t)(w[0] >> (8 * i)); + note(" ref keystream[0..3] = %02x %02x %02x %02x " + "(RFC 8439 expects 10 f1 e7 e4)\n", + ks0[0], ks0[1], ks0[2], ks0[3]); + ok("scalar ChaCha20 matches RFC 8439 vector", + ks0[0] == 0x10 && ks0[1] == 0xf1 && + ks0[2] == 0xe7 && ks0[3] == 0xe4); + } + + /* (2) validate the NEON kernel against that reference. The kernel + * hardwires nonce = 0 and starts the block counter at 0, so we + * compare its keystream to the reference block-for-block. */ + { + uint8_t buf[128]; + uint8_t zero_nonce[12] = {0}; + uint32_t ref0[16], ref1[16]; + int match = 1; + + for (i = 0; i < 32; i++) + key[i] = (uint8_t)(i * 5 + 1); + memset(buf, 0, sizeof buf); /* zeros -> raw keystream */ + fm_chacha20(buf, sizeof buf, key, 1); + + ref_chacha_block(ref0, key, 0, zero_nonce); + ref_chacha_block(ref1, key, 1, zero_nonce); + for (i = 0; i < 16; i++) { + uint32_t k0 = (uint32_t)buf[4 * i] | + (uint32_t)buf[4 * i + 1] << 8 | + (uint32_t)buf[4 * i + 2] << 16 | + (uint32_t)buf[4 * i + 3] << 24; + uint32_t k1 = (uint32_t)buf[64 + 4 * i] | + (uint32_t)buf[64 + 4 * i + 1] << 8 | + (uint32_t)buf[64 + 4 * i + 2] << 16 | + (uint32_t)buf[64 + 4 * i + 3] << 24; + if (k0 != ref0[i] || k1 != ref1[i]) + match = 0; + } + ok("NEON ChaCha20 matches scalar reference (2 blocks)", match); + } + + /* (3) the cipher is a real XOR stream: applying it twice is identity */ + { + uint8_t plain[128], work[128], k2[32]; + for (i = 0; i < 128; i++) + plain[i] = (uint8_t)(i * 7 + 1); + for (i = 0; i < 32; i++) + k2[i] = (uint8_t)(i * 3); + memcpy(work, plain, 128); + fm_chacha20(work, 128, k2, 1); + ok("chacha20 actually changes data", memcmp(work, plain, 128) != 0); + fm_chacha20(work, 128, k2, 1); + ok("chacha20 round-trips (XOR is involutive)", + memcmp(work, plain, 128) == 0); + } +} + +static void check_physics(void) +{ + /* two equal masses released from rest must accelerate toward each + * other: symmetric, momentum-conserving, and bounded. */ + double bodies[2 * 8] = {0}; + double total_p; + + bodies[0] = -1.0; bodies[3] = 1.0; /* body 0 at x=-1, mass 1 */ + bodies[8] = 1.0; bodies[11] = 1.0; /* body 1 at x=+1, mass 1 */ + + fm_physics(bodies, 2, 200); + + /* velocities must be equal and opposite (Newton's third law) */ + total_p = bodies[4] + bodies[12]; /* vx0 + vx1 */ + note(" 2-body: vx0=%.6f vx1=%.6f (sum should be ~0)\n", + bodies[4], bodies[12]); + ok("physics conserves momentum", fabs(total_p) < 1e-9); + ok("physics: bodies attract", bodies[4] > 0.0 && bodies[12] < 0.0); + ok("physics values stay finite", isfinite(bodies[0]) && isfinite(bodies[8])); +} + +static int cmp_u32(const void *p, const void *q) +{ + uint32_t x = *(const uint32_t *)p, y = *(const uint32_t *)q; + return (x > y) - (x < y); +} + +static int is_sorted(const uint32_t *a, size_t n) +{ + for (size_t i = 1; i < n; i++) + if (a[i - 1] > a[i]) + return 0; + return 1; +} + +static void check_sort(void) +{ + enum { N = 10000 }; + uint32_t *a = malloc(N * sizeof(uint32_t)); + uint32_t *b = malloc(N * sizeof(uint32_t)); + uint64_t s; + size_t i; + uint32_t r = 12345; + + for (i = 0; i < N; i++) { + r = r * 1103515245u + 12345u; + a[i] = r; + } + memcpy(b, a, N * sizeof(uint32_t)); + + s = fm_sort(a, N); + ok("sort produces sorted output", is_sorted(a, N)); + + /* multiset is preserved: sort the reference with the C library and + * compare element by element */ + qsort(b, N, sizeof(uint32_t), cmp_u32); + ok("sort is a permutation of the input", + memcmp(a, b, N * sizeof(uint32_t)) == 0); + + /* already-sorted input stays sorted and gives the same checksum */ + { + uint64_t s2 = fm_sort(a, N); + ok("sort is idempotent on sorted data", + is_sorted(a, N) && s2 == s); + } + + ok("sort of empty array is zero", fm_sort(a, 0) == 0); + free(a); + free(b); +} + +static void check_chase(void) +{ + /* build a tiny 4-node cycle by hand and confirm the walk returns to + * the start after exactly `n` steps (offset 0 relative to entry) */ + void *nodes[4]; + + nodes[0] = &nodes[1]; + nodes[1] = &nodes[2]; + nodes[2] = &nodes[3]; + nodes[3] = &nodes[0]; + + /* 4 hops from &nodes[0] returns to &nodes[0]; fm_chase returns the + * final pointer minus the starting pointer, so a full loop gives 0 */ + ok("chase completes a full cycle", fm_chase(nodes, 4) == 0); + ok("chase(0) is zero", fm_chase(nodes, 0) == 0); + /* one hop lands on &nodes[1], i.e. one pointer-width past the start */ + ok("chase single hop offset", + fm_chase(nodes, 1) == (uint64_t)((char *)&nodes[1] - (char *)&nodes[0])); +} + +int main(void) +{ + long n = sysconf(_SC_NPROCESSORS_ONLN); + + fm_ncores = n > 0 ? n : 1; + + printf("\nfossmark kernel correctness tests\n"); + printf("=================================\n"); + printf("running each check on %ld core%s in parallel\n\n", + fm_ncores, fm_ncores == 1 ? "" : "s"); + + printf("Integer Math:\n"); parallel(check_int); + printf("Floating Point Math:\n"); parallel(check_fp); + printf("Prime Numbers:\n"); parallel(check_primes); + printf("Extended Instructions:\n"); parallel(check_simd); + printf("Compression:\n"); parallel(check_compress); + printf("Encryption:\n"); parallel(check_crypto); + printf("Physics:\n"); parallel(check_physics); + printf("Sorting:\n"); parallel(check_sort); + /* the pointer chase is the single-threaded test: run it on one core */ + printf("Single-Threaded (chase):\n"); check_chase(); + + printf("\n=================================\n"); + printf("%d checks, %d failures\n\n", checks, failures); + return failures ? 1 : 0; +}