From 61a3ea22e4afe9db89ed555afa4c503c97bbe9a3 Mon Sep 17 00:00:00 2001 From: tannevaled Date: Sat, 7 Nov 2026 22:20:12 +0100 Subject: [PATCH 1/2] kernels: import go-simd/floats for sum, elementwise, extrema, clip, exp, log (K1) K1 of the go-scikit performance campaign (go-scikit/roadmap#27): ndarray now imports go-simd/floats v0.9.0 and deletes its copies. - exp, log, and the correctly rounded float32 exp/log (with the AVX2 kernels of #94/#96 and the NEON kernel of #95, now go-simd/floats#29): deleted here (Go sources, tests, generator kernels, assembly); Exp, Log, Exp32, Log32 and their P forms call floats. - amd64 and arm64: the sum/sqrt/max/min/add/sub/mul/div kernels are floats'; Clip is floats.ClampScalar. Max/Min resolve a zero result with the serial oracle (floats leaves the sign of zero unspecified). - Sum is floats.SumPairwise on every architecture: NumPy's order, bit-identical to np.sum; serial at every size. A changed last bit is the fix; the golden values of TestSumMatchesNumPy are NumPy 2.5.4's. - Not moved: the float32 elementwise/sqrt/extremum/AddBlocks8F32 kernels, the dot inside GEMM/Dot1D, and the ppc64le, loong64, riscv64 and s390x kernels (floats has none for them yet). Co-Authored-By: Claude Sonnet 5.5 --- README.md | 13 + go.mod | 7 +- go.sum | 4 + internal/kernels/asmgen/amd64/gen.go | 17 +- internal/kernels/asmgen/arm64/gen.go | 11 +- internal/kernels/exp.go | 270 --------- internal/kernels/exp_test.go | 178 ------ internal/kernels/expcr.go | 94 --- internal/kernels/expcr_amd64.go | 50 -- internal/kernels/expcr_amd64_test.go | 45 -- internal/kernels/expcr_fence_amd64_test.go | 20 - internal/kernels/expcr_other.go | 10 - internal/kernels/expcr_test.go | 53 -- internal/kernels/explog.go | 83 +++ internal/kernels/generate.go | 10 + internal/kernels/kernels.go | 16 +- internal/kernels/log.go | 305 ---------- internal/kernels/log_amd64.go | 54 -- internal/kernels/log_amd64_test.go | 84 --- internal/kernels/log_fence_amd64_test.go | 20 - internal/kernels/log_other.go | 6 - internal/kernels/log_test.go | 154 ----- internal/kernels/logcr.go | 138 ----- internal/kernels/logcr_amd64.go | 53 -- internal/kernels/logcr_amd64_test.go | 91 --- internal/kernels/logcr_fence_amd64_test.go | 20 - internal/kernels/logcr_test.go | 150 ----- internal/kernels/parallel.go | 14 +- internal/kernels/reduce_amd64.go | 124 ---- internal/kernels/reduce_arm64.go | 118 ---- internal/kernels/reduce_floats.go | 51 ++ internal/kernels/reduce_generic.go | 1 - internal/kernels/reduce_loong64.go | 8 +- internal/kernels/reduce_ppc64le.go | 10 +- internal/kernels/reduce_riscv64.go | 17 +- internal/kernels/reduce_s390x.go | 8 +- internal/kernels/reduce_sum.go | 12 + internal/kernels/sum_amd64.s | 630 --------------------- internal/kernels/sum_arm64.s | 243 -------- internal/kernels/sum_numpy_test.go | 37 ++ internal/kernels/via64.go | 37 -- internal/kernels/via64_exhaustive_test.go | 267 --------- internal/kernels/via64_oracle_test.go | 68 --- internal/kernels/via64_test.go | 146 ----- 44 files changed, 248 insertions(+), 3499 deletions(-) delete mode 100644 internal/kernels/exp.go delete mode 100644 internal/kernels/exp_test.go delete mode 100644 internal/kernels/expcr.go delete mode 100644 internal/kernels/expcr_amd64.go delete mode 100644 internal/kernels/expcr_amd64_test.go delete mode 100644 internal/kernels/expcr_fence_amd64_test.go delete mode 100644 internal/kernels/expcr_other.go delete mode 100644 internal/kernels/expcr_test.go create mode 100644 internal/kernels/explog.go create mode 100644 internal/kernels/generate.go delete mode 100644 internal/kernels/log.go delete mode 100644 internal/kernels/log_amd64.go delete mode 100644 internal/kernels/log_amd64_test.go delete mode 100644 internal/kernels/log_fence_amd64_test.go delete mode 100644 internal/kernels/log_other.go delete mode 100644 internal/kernels/log_test.go delete mode 100644 internal/kernels/logcr.go delete mode 100644 internal/kernels/logcr_amd64.go delete mode 100644 internal/kernels/logcr_amd64_test.go delete mode 100644 internal/kernels/logcr_fence_amd64_test.go delete mode 100644 internal/kernels/logcr_test.go delete mode 100644 internal/kernels/reduce_amd64.go delete mode 100644 internal/kernels/reduce_arm64.go create mode 100644 internal/kernels/reduce_floats.go create mode 100644 internal/kernels/reduce_sum.go create mode 100644 internal/kernels/sum_numpy_test.go delete mode 100644 internal/kernels/via64.go delete mode 100644 internal/kernels/via64_exhaustive_test.go delete mode 100644 internal/kernels/via64_oracle_test.go delete mode 100644 internal/kernels/via64_test.go diff --git a/README.md b/README.md index e7d4332..5b6613d 100644 --- a/README.md +++ b/README.md @@ -175,6 +175,19 @@ arch is therefore a durable foundation. The numeric loops live in on every 64-bit target Go supports (amd64, arm64, ppc64le, loong64, riscv64, s390x), behind the same API and held to the same tests. +The one-dimensional float64 kernels that other libraries need too live once, in +[go-simd/floats](https://github.com/go-simd/floats) (K1 of the go-scikit +[performance campaign](https://github.com/go-scikit/roadmap/blob/main/docs/performance-campaign.md)), +which this module imports: on amd64 (AVX2) and arm64 (NEON) the elementwise +`Add`/`Sub`/`Mul`/`Div`, `Sqrt`, `Max`/`Min` and `Clip` kernels, and everywhere +`Exp`, `Log` and the correctly rounded float32 `Exp`/`Log` (the AVX2 and NEON +kernels were written here and moved down). **`Sum` is NumPy's pairwise sum** +(`floats.SumPairwise`): bit-identical to `np.sum` of a contiguous 1-D float64 +array, on every architecture, at every size. It replaces lane-parallel sums whose +last bits differed from NumPy's and from one architecture to the next, and it is +serial (a multicore split would regroup NumPy's order). The other architectures +keep their own elementwise kernels here until floats has them. + ## Goroutines Operations on large arrays are spread over `GOMAXPROCS` goroutines. After the diff --git a/go.mod b/go.mod index c1e9113..ac8b53d 100644 --- a/go.mod +++ b/go.mod @@ -2,4 +2,9 @@ module github.com/go-ndarray/ndarray go 1.27.2 -require github.com/go-fft/fft v0.28.11 +require ( + github.com/go-fft/fft v0.28.11 + github.com/go-simd/floats v0.9.0 +) + +require golang.org/x/sys v0.49.0 // indirect diff --git a/go.sum b/go.sum index 65bd968..6671d94 100644 --- a/go.sum +++ b/go.sum @@ -1,2 +1,6 @@ github.com/go-fft/fft v0.28.11 h1:ZfT78fSmPjZU2E9PgkSI/0f0E3sesYpl13GHmIFzRGc= github.com/go-fft/fft v0.28.11/go.mod h1:GN5V4TWCjV05ECqH1ONfjFM6Z1m2shqCdssxKGqQ41s= +github.com/go-simd/floats v0.9.0 h1:P/sMlwYWwaqCah/cMcBRfffK1JcUxXZVj76f4AwQRN0= +github.com/go-simd/floats v0.9.0/go.mod h1:Aq5a6oY5pUqm9Hpi+iV2nrNmaaUb0rwHvUHIerXWSvA= +golang.org/x/sys v0.49.0 h1:XbzkgYJHdqh/8m2Uu0W/dQv8nktxx4BFHp1M0gROTXA= +golang.org/x/sys v0.49.0/go.mod h1:hNLxWAXmnKAxqDtdwIYC4bM9oQPEecfsnNMuSxOs3og= diff --git a/internal/kernels/asmgen/amd64/gen.go b/internal/kernels/asmgen/amd64/gen.go index 21df7f3..c698fd6 100644 --- a/internal/kernels/asmgen/amd64/gen.go +++ b/internal/kernels/asmgen/amd64/gen.go @@ -37,14 +37,11 @@ import ( func main() { f := emit.NewFile("amd64") - f.Add(sumKernel()) - f.Add(sqrtKernel()) - f.Add(extremeKernel("maxSSE2", "MAXPD", "MAXSD")) - f.Add(extremeKernel("minSSE2", "MINPD", "MINSD")) - f.Add(binKernel("addSSE2", "ADDPD", "ADDSD")) - f.Add(binKernel("subSSE2", "SUBPD", "SUBSD")) - f.Add(binKernel("mulSSE2", "MULPD", "MULSD")) - f.Add(binKernel("divSSE2", "DIVPD", "DIVSD")) + // The float64 sum, sqrt, max/min and add/sub/mul/div kernels, and the float32 + // exp/log and float64 log kernels, moved to go-simd/floats (K1 of the + // go-scikit performance campaign): sumKernel, sqrtKernel, extremeKernel, + // binKernel, exp32CRKernel, log32CRKernel and logKernel below are no longer + // emitted. f.Add(gemmFMAKernel()) f.Add(gemm32FMAKernel()) for mr := 1; mr <= gemmMR; mr++ { @@ -54,10 +51,6 @@ func main() { f.Add(dotFMAKernel()) f.Add(gemmSSE2HalfKernel()) f.Add(amd64.FeatureProbe("hasFMA", amd64.FMA)) - f.Add(amd64.FeatureProbe("hasAVX2", amd64.AVX2)) - f.Add(exp32CRKernel()) - f.Add(log32CRKernel()) - f.Add(logKernel()) f.Add(blocks8Kernel()) f.Add(bin32Kernel("add32SSE2", "ADDPS", "ADDSS")) f.Add(bin32Kernel("sub32SSE2", "SUBPS", "SUBSS")) diff --git a/internal/kernels/asmgen/arm64/gen.go b/internal/kernels/asmgen/arm64/gen.go index d07eb3f..ebd7d73 100644 --- a/internal/kernels/asmgen/arm64/gen.go +++ b/internal/kernels/asmgen/arm64/gen.go @@ -41,20 +41,15 @@ func reduceSig() arm64.Signature { func main() { f := emit.NewFile("arm64") - f.Add(sumKernel()) + // The float64 sum, sqrt, add/sub/mul/div and max/min kernels moved to + // go-simd/floats (K1 of the go-scikit performance campaign): sumKernel, + // sqrtKernel, binKernel and extremumKernel below are no longer emitted. f.Add(dotKernel()) - f.Add(sqrtKernel()) - f.Add(binKernel("addNEON", "VFADD", "FADDD")) - f.Add(binKernel("subNEON", "VFSUB", "FSUBD")) - f.Add(binKernel("mulNEON", "VFMUL", "FMULD")) - f.Add(binKernel("divNEON", "VFDIV", "FDIVD")) f.Add(bin32Kernel("add32NEON", "VFADD", "FADDS")) f.Add(bin32Kernel("sub32NEON", "VFSUB", "FSUBS")) f.Add(bin32Kernel("mul32NEON", "VFMUL", "FMULS")) f.Add(bin32Kernel("div32NEON", "VFDIV", "FDIVS")) f.Add(sqrt32Kernel()) - f.Add(extremumKernel("maxNEON", "VFMAX", "FMAXD")) - f.Add(extremumKernel("minNEON", "VFMIN", "FMIND")) f.Add(extremum32Kernel("max32NEON", "VFMAX", "FMAXS")) f.Add(extremum32Kernel("min32NEON", "VFMIN", "FMINS")) f.Add(gemmKernel()) diff --git a/internal/kernels/exp.go b/internal/kernels/exp.go deleted file mode 100644 index e2511ec..0000000 --- a/internal/kernels/exp.go +++ /dev/null @@ -1,270 +0,0 @@ -package kernels - -import "math" - -// Exp is a port of the double-precision exp of Arm's optimized-routines -// (math/exp.c and math/exp_data.c, by Szabolcs Nagy; also glibc's exp since -// 2.28), used under its MIT licence: -// -// Copyright (c) 2018-2025, Arm Limited. -// SPDX-License-Identifier: MIT OR Apache-2.0 WITH LLVM-exception -// -// exp(x) = 2^(k/N) * exp(r) with N = 128, x = k*ln2/N + r, |r| <= ln2/2N. -// 2^(k/N) comes from a 128-entry table as scale*(1+tail), and exp(r)-1 from a -// degree-5 polynomial; the result is scale + scale*(tail + exp(r)-1). Arm -// documents a worst case of 0.509 ULP (0.511 without FMA); TestExpULP measures -// this port against a 256-bit reference. It is branch-light and table-driven, -// unlike math.Exp, which is what makes it several times faster per element. - -const ( - expTableBits = 7 - expN = 1 << expTableBits - expInvLn2N = 0x1.71547652b82fep0 * expN - expNegLn2hiN = -0x1.62e42fefa0000p-8 - expNegLn2loN = -0x1.cf79abc9e3b3ap-47 - expShift = 0x1.8p52 - expC2 = 0x1.ffffffffffdbdp-2 - expC3 = 0x1.555555555543cp-3 - expC4 = 0x1.55555cf172b91p-5 - expC5 = 0x1.1111167a4d017p-7 -) - -// expTab holds, for k in [0, N): expTab[2k] = the bits of tail[k], and -// expTab[2k+1] = the bits of H[k] minus k<<45, where 2^(k/N) ~= H[k]*(1+tail[k]). -// Copied from exp_data.c (N == 128); TestExpTable re-derives every entry. -var expTab = [2 * expN]uint64{ - 0x0, 0x3ff0000000000000, - 0x3c9b3b4f1a88bf6e, 0x3feff63da9fb3335, - 0xbc7160139cd8dc5d, 0x3fefec9a3e778061, - 0xbc905e7a108766d1, 0x3fefe315e86e7f85, - 0x3c8cd2523567f613, 0x3fefd9b0d3158574, - 0xbc8bce8023f98efa, 0x3fefd06b29ddf6de, - 0x3c60f74e61e6c861, 0x3fefc74518759bc8, - 0x3c90a3e45b33d399, 0x3fefbe3ecac6f383, - 0x3c979aa65d837b6d, 0x3fefb5586cf9890f, - 0x3c8eb51a92fdeffc, 0x3fefac922b7247f7, - 0x3c3ebe3d702f9cd1, 0x3fefa3ec32d3d1a2, - 0xbc6a033489906e0b, 0x3fef9b66affed31b, - 0xbc9556522a2fbd0e, 0x3fef9301d0125b51, - 0xbc5080ef8c4eea55, 0x3fef8abdc06c31cc, - 0xbc91c923b9d5f416, 0x3fef829aaea92de0, - 0x3c80d3e3e95c55af, 0x3fef7a98c8a58e51, - 0xbc801b15eaa59348, 0x3fef72b83c7d517b, - 0xbc8f1ff055de323d, 0x3fef6af9388c8dea, - 0x3c8b898c3f1353bf, 0x3fef635beb6fcb75, - 0xbc96d99c7611eb26, 0x3fef5be084045cd4, - 0x3c9aecf73e3a2f60, 0x3fef54873168b9aa, - 0xbc8fe782cb86389d, 0x3fef4d5022fcd91d, - 0x3c8a6f4144a6c38d, 0x3fef463b88628cd6, - 0x3c807a05b0e4047d, 0x3fef3f49917ddc96, - 0x3c968efde3a8a894, 0x3fef387a6e756238, - 0x3c875e18f274487d, 0x3fef31ce4fb2a63f, - 0x3c80472b981fe7f2, 0x3fef2b4565e27cdd, - 0xbc96b87b3f71085e, 0x3fef24dfe1f56381, - 0x3c82f7e16d09ab31, 0x3fef1e9df51fdee1, - 0xbc3d219b1a6fbffa, 0x3fef187fd0dad990, - 0x3c8b3782720c0ab4, 0x3fef1285a6e4030b, - 0x3c6e149289cecb8f, 0x3fef0cafa93e2f56, - 0x3c834d754db0abb6, 0x3fef06fe0a31b715, - 0x3c864201e2ac744c, 0x3fef0170fc4cd831, - 0x3c8fdd395dd3f84a, 0x3feefc08b26416ff, - 0xbc86a3803b8e5b04, 0x3feef6c55f929ff1, - 0xbc924aedcc4b5068, 0x3feef1a7373aa9cb, - 0xbc9907f81b512d8e, 0x3feeecae6d05d866, - 0xbc71d1e83e9436d2, 0x3feee7db34e59ff7, - 0xbc991919b3ce1b15, 0x3feee32dc313a8e5, - 0x3c859f48a72a4c6d, 0x3feedea64c123422, - 0xbc9312607a28698a, 0x3feeda4504ac801c, - 0xbc58a78f4817895b, 0x3feed60a21f72e2a, - 0xbc7c2c9b67499a1b, 0x3feed1f5d950a897, - 0x3c4363ed60c2ac11, 0x3feece086061892d, - 0x3c9666093b0664ef, 0x3feeca41ed1d0057, - 0x3c6ecce1daa10379, 0x3feec6a2b5c13cd0, - 0x3c93ff8e3f0f1230, 0x3feec32af0d7d3de, - 0x3c7690cebb7aafb0, 0x3feebfdad5362a27, - 0x3c931dbdeb54e077, 0x3feebcb299fddd0d, - 0xbc8f94340071a38e, 0x3feeb9b2769d2ca7, - 0xbc87deccdc93a349, 0x3feeb6daa2cf6642, - 0xbc78dec6bd0f385f, 0x3feeb42b569d4f82, - 0xbc861246ec7b5cf6, 0x3feeb1a4ca5d920f, - 0x3c93350518fdd78e, 0x3feeaf4736b527da, - 0x3c7b98b72f8a9b05, 0x3feead12d497c7fd, - 0x3c9063e1e21c5409, 0x3feeab07dd485429, - 0x3c34c7855019c6ea, 0x3feea9268a5946b7, - 0x3c9432e62b64c035, 0x3feea76f15ad2148, - 0xbc8ce44a6199769f, 0x3feea5e1b976dc09, - 0xbc8c33c53bef4da8, 0x3feea47eb03a5585, - 0xbc845378892be9ae, 0x3feea34634ccc320, - 0xbc93cedd78565858, 0x3feea23882552225, - 0x3c5710aa807e1964, 0x3feea155d44ca973, - 0xbc93b3efbf5e2228, 0x3feea09e667f3bcd, - 0xbc6a12ad8734b982, 0x3feea012750bdabf, - 0xbc6367efb86da9ee, 0x3fee9fb23c651a2f, - 0xbc80dc3d54e08851, 0x3fee9f7df9519484, - 0xbc781f647e5a3ecf, 0x3fee9f75e8ec5f74, - 0xbc86ee4ac08b7db0, 0x3fee9f9a48a58174, - 0xbc8619321e55e68a, 0x3fee9feb564267c9, - 0x3c909ccb5e09d4d3, 0x3feea0694fde5d3f, - 0xbc7b32dcb94da51d, 0x3feea11473eb0187, - 0x3c94ecfd5467c06b, 0x3feea1ed0130c132, - 0x3c65ebe1abd66c55, 0x3feea2f336cf4e62, - 0xbc88a1c52fb3cf42, 0x3feea427543e1a12, - 0xbc9369b6f13b3734, 0x3feea589994cce13, - 0xbc805e843a19ff1e, 0x3feea71a4623c7ad, - 0xbc94d450d872576e, 0x3feea8d99b4492ed, - 0x3c90ad675b0e8a00, 0x3feeaac7d98a6699, - 0x3c8db72fc1f0eab4, 0x3feeace5422aa0db, - 0xbc65b6609cc5e7ff, 0x3feeaf3216b5448c, - 0x3c7bf68359f35f44, 0x3feeb1ae99157736, - 0xbc93091fa71e3d83, 0x3feeb45b0b91ffc6, - 0xbc5da9b88b6c1e29, 0x3feeb737b0cdc5e5, - 0xbc6c23f97c90b959, 0x3feeba44cbc8520f, - 0xbc92434322f4f9aa, 0x3feebd829fde4e50, - 0xbc85ca6cd7668e4b, 0x3feec0f170ca07ba, - 0x3c71affc2b91ce27, 0x3feec49182a3f090, - 0x3c6dd235e10a73bb, 0x3feec86319e32323, - 0xbc87c50422622263, 0x3feecc667b5de565, - 0x3c8b1c86e3e231d5, 0x3feed09bec4a2d33, - 0xbc91bbd1d3bcbb15, 0x3feed503b23e255d, - 0x3c90cc319cee31d2, 0x3feed99e1330b358, - 0x3c8469846e735ab3, 0x3feede6b5579fdbf, - 0xbc82dfcd978e9db4, 0x3feee36bbfd3f37a, - 0x3c8c1a7792cb3387, 0x3feee89f995ad3ad, - 0xbc907b8f4ad1d9fa, 0x3feeee07298db666, - 0xbc55c3d956dcaeba, 0x3feef3a2b84f15fb, - 0xbc90a40e3da6f640, 0x3feef9728de5593a, - 0xbc68d6f438ad9334, 0x3feeff76f2fb5e47, - 0xbc91eee26b588a35, 0x3fef05b030a1064a, - 0x3c74ffd70a5fddcd, 0x3fef0c1e904bc1d2, - 0xbc91bdfbfa9298ac, 0x3fef12c25bd71e09, - 0x3c736eae30af0cb3, 0x3fef199bdd85529c, - 0x3c8ee3325c9ffd94, 0x3fef20ab5fffd07a, - 0x3c84e08fd10959ac, 0x3fef27f12e57d14b, - 0x3c63cdaf384e1a67, 0x3fef2f6d9406e7b5, - 0x3c676b2c6c921968, 0x3fef3720dcef9069, - 0xbc808a1883ccb5d2, 0x3fef3f0b555dc3fa, - 0xbc8fad5d3ffffa6f, 0x3fef472d4a07897c, - 0xbc900dae3875a949, 0x3fef4f87080d89f2, - 0x3c74a385a63d07a7, 0x3fef5818dcfba487, - 0xbc82919e2040220f, 0x3fef60e316c98398, - 0x3c8e5a50d5c192ac, 0x3fef69e603db3285, - 0x3c843a59ac016b4b, 0x3fef7321f301b460, - 0xbc82d52107b43e1f, 0x3fef7c97337b9b5f, - 0xbc892ab93b470dc9, 0x3fef864614f5a129, - 0x3c74b604603a88d3, 0x3fef902ee78b3ff6, - 0x3c83c5ec519d7271, 0x3fef9a51fbc74c83, - 0xbc8ff7128fd391f0, 0x3fefa4afa2a490da, - 0xbc8dae98e223747d, 0x3fefaf482d8e67f1, - 0x3c8ec3bc41aa2008, 0x3fefba1bee615a27, - 0x3c842b94c3a9eb32, 0x3fefc52b376bba97, - 0x3c8a64a931d185ee, 0x3fefd0765b6e4540, - 0xbc8e37bae43be3ed, 0x3fefdbfdad9cbe14, - 0x3c77893b4d91cd9d, 0x3fefe7c1819e90d8, - 0x3c5305c14160cc89, 0x3feff3c22b8f71f1, -} - -// top12 is the sign and exponent bits of x. -func top12(x float64) uint32 { return uint32(math.Float64bits(x) >> 52) } - -// expSpecial finishes exp when scale*(1+tmp) may overflow or underflow: the -// exponent in sbits may have wrapped, k's sign (bit 31 of ki) tells which. -// The subnormal branch rounds y once before scaling it down, avoiding the -// double rounding that would cost up to 0.5 ULP more. Conversions to float64 -// keep the compiler from fusing the compensated sums into an FMA. -func expSpecial(tmp float64, sbits, ki uint64) float64 { - if ki&0x80000000 == 0 { - // k > 0: the exponent of scale might have overflowed by <= 460. - sbits -= 1009 << 52 - scale := math.Float64frombits(sbits) - return 0x1p1009 * float64(scale+float64(scale*tmp)) - } - // k < 0: care is needed in the subnormal range. - sbits += 1022 << 52 - scale := math.Float64frombits(sbits) - y := scale + float64(scale*tmp) - if y < 1.0 { - lo := scale - y + float64(scale*tmp) - hi := 1.0 + y - lo = 1.0 - hi + y + lo - y = float64(hi+lo) - 1.0 - if y == 0 { - y = 0 // never -0 - } - } - return 0x1p-1022 * y -} - -// exp1 is exp(x) for one element (exp_inline of exp.c with xtail = 0). -func exp1(x float64) float64 { - abstop := top12(x) & 0x7ff - if abstop-0x3c9 >= 0x408-0x3c9 { // |x| < 2^-54, or |x| >= 512, or not finite - if abstop-0x3c9 >= 0x80000000 { - return 1.0 + x // tiny x (and 0, a common input) - } - if abstop >= 0x409 { // |x| >= 1024, Inf or NaN - switch { - case x == math.Inf(-1): - return 0 - case abstop >= 0x7ff: - return 1.0 + x // NaN, +Inf - case math.Signbit(x): - return 0 // underflow - default: - return math.Inf(1) // overflow - } - } - abstop = 0 // large |x|: finished by expSpecial - } - // z - kd is in [-1, 1] in non-nearest rounding modes. - z := expInvLn2N * x - kd := float64(z + expShift) - ki := math.Float64bits(kd) - kd -= expShift - r := x + kd*expNegLn2hiN + kd*expNegLn2loN - idx := 2 * (ki % expN) - top := ki << (52 - expTableBits) - tail := math.Float64frombits(expTab[idx]) - sbits := expTab[idx+1] + top - r2 := r * r - tmp := tail + r + r2*(expC2+r*expC3) + r2*r2*(expC4+r*expC5) - if abstop == 0 { - return expSpecial(tmp, sbits, ki) - } - scale := math.Float64frombits(sbits) - return scale + scale*tmp -} - -// Exp writes exp(src[i]) into dst[i]. exp1 is too large to inline, and a -// call per element cost half the time, so the common case (2^-54 <= |x| < -// 512) is exp1's main path written out in the loop; everything else calls -// exp1. TestExpLoopMatchesExp1 holds the two to the same bits. -func Exp(dst, src []float64) { - dst = dst[:len(src)] - for i, x := range src { - if top12(x)&0x7ff-0x3c9 >= 0x408-0x3c9 { - dst[i] = exp1(x) - continue - } - z := expInvLn2N * x - kd := float64(z + expShift) - ki := math.Float64bits(kd) - kd -= expShift - r := x + kd*expNegLn2hiN + kd*expNegLn2loN - idx := 2 * (ki % expN) - tail := math.Float64frombits(expTab[idx]) - scale := math.Float64frombits(expTab[idx+1] + ki<<(52-expTableBits)) - r2 := r * r - tmp := tail + r + r2*(expC2+r*expC3) + r2*r2*(expC4+r*expC5) - dst[i] = scale + scale*tmp - } -} - -// ExpP is Exp, split across cores above ParThreshold. -func ExpP(dst, src []float64) { - n := len(dst) - if n < ParThreshold { - Exp(dst, src) - return - } - parallelFor(n, numWorkers(n), func(lo, hi int) { Exp(dst[lo:hi], src[lo:hi]) }) -} diff --git a/internal/kernels/exp_test.go b/internal/kernels/exp_test.go deleted file mode 100644 index 3c06be7..0000000 --- a/internal/kernels/exp_test.go +++ /dev/null @@ -1,178 +0,0 @@ -package kernels - -import ( - "math" - "math/big" - "math/rand" - "testing" -) - -const refPrec = 300 - -// bigLn2 is ln 2 to refPrec bits: 2*atanh(1/3) = 2*sum (1/3)^(2n+1)/(2n+1). -func bigLn2() *big.Float { - third := new(big.Float).SetPrec(refPrec).Quo(big.NewFloat(1).SetPrec(refPrec), big.NewFloat(3)) - ninth := new(big.Float).SetPrec(refPrec).Mul(third, third) - sum := new(big.Float).SetPrec(refPrec) - pow := new(big.Float).SetPrec(refPrec).Set(third) - for n := 0; n < 300; n++ { - term := new(big.Float).SetPrec(refPrec).Quo(pow, big.NewFloat(float64(2*n+1))) - sum.Add(sum, term) - pow.Mul(pow, ninth) - } - return sum.Mul(sum, big.NewFloat(2)) -} - -var ln2Big = bigLn2() - -// refExp is exp(x) to refPrec bits, independent of the code under test: -// x = k*ln2 + r with |r| <= ln2/2, exp(r) by its Taylor series, times 2^k. -func refExp(x float64) *big.Float { - bx := new(big.Float).SetPrec(refPrec).SetFloat64(x) - q := new(big.Float).SetPrec(refPrec).Quo(bx, ln2Big) - qf, _ := q.Float64() - k := math.Round(qf) - r := new(big.Float).SetPrec(refPrec).Mul(ln2Big, big.NewFloat(k)) - r.Sub(bx, r) - sum := new(big.Float).SetPrec(refPrec).SetInt64(1) - term := new(big.Float).SetPrec(refPrec).SetInt64(1) - for n := 1; n < 80; n++ { - term.Mul(term, r) - term.Quo(term, big.NewFloat(float64(n))) - sum.Add(sum, term) - } - return sum.SetMantExp(sum, int(k)) -} - -// ulpErr is |got - exact| in units of the float64 spacing at exact (2^-1074 -// in the subnormal range). -func ulpErr(got float64, exact *big.Float) float64 { - e := exact.MantExp(nil) - 1 // exact = m * 2^e, 1 <= m < 2 - ulpExp := max(e-52, -1074) - d := new(big.Float).SetPrec(refPrec).SetFloat64(got) - d.Sub(d, exact) - d.SetMantExp(d, -ulpExp) - f, _ := d.Float64() - return math.Abs(f) -} - -// TestExpTable re-derives every table entry from 2^(k/128) computed to 300 -// bits (2 square-rooted seven times, then raised to k): H[k] must be the double -// nearest 2^(k/N) and tail[k] the double nearest (2^(k/N) - H)/H. -func TestExpTable(t *testing.T) { - root := new(big.Float).SetPrec(refPrec).SetInt64(2) - for i := 0; i < expTableBits; i++ { - root.Sqrt(root) - } - e := new(big.Float).SetPrec(refPrec).SetInt64(1) - for k := 0; k < expN; k++ { - h := math.Float64frombits(expTab[2*k+1] + uint64(k)<<45) - if want, _ := e.Float64(); h != want { - t.Fatalf("H[%d] = %v, nearest double to 2^(k/N) is %v", k, h, want) - } - tail := new(big.Float).SetPrec(refPrec).SetFloat64(h) - tail.Sub(e, tail) - tail.Quo(tail, new(big.Float).SetFloat64(h)) - got := math.Float64frombits(expTab[2*k]) - if want, _ := tail.Float64(); got != want { - t.Fatalf("tail[%d] = %v, want %v", k, got, want) - } - e.Mul(e, root) - } -} - -// TestExpULP measures the worst error against the 256-bit reference over the -// whole finite domain, the region near 0, the overflow edge and the subnormal -// range, and requires it to stay within Arm's documented 0.511 ULP (+ margin). -func TestExpULP(t *testing.T) { - rng := rand.New(rand.NewSource(1)) - var xs []float64 - n := 5000 - if testing.Short() { - n = 2000 - } - for i := 0; i < n; i++ { - xs = append(xs, - rng.Float64()*1454.78-745, // the whole finite range - (rng.Float64()*2-1)*math.Ldexp(1, -rng.Intn(60)), // near 0, log-uniform - 709+rng.Float64()*0.7827, // just below overflow - -745.13+rng.Float64()*37, // subnormal and near-subnormal results - ) - } - worst, worstX, goWorst := 0.0, 0.0, 0.0 - for _, x := range xs { - ref := refExp(x) - if e := ulpErr(exp1(x), ref); e > worst { - worst, worstX = e, x - } - goWorst = max(goWorst, ulpErr(math.Exp(x), ref)) - } - t.Logf("exp: worst %.4f ULP at x=%v over %d inputs (math.Exp: %.4f)", worst, worstX, len(xs), goWorst) - if worst > 0.52 { - t.Errorf("exp error %.4f ULP at x=%v exceeds 0.52", worst, worstX) - } -} - -// TestExpLoopMatchesExp1: Exp's inlined main path and exp1 give the same bits -// for every input, special or not. -func TestExpLoopMatchesExp1(t *testing.T) { - rng := rand.New(rand.NewSource(2)) - xs := []float64{0, 1e-300, math.Ldexp(1, -54), math.Ldexp(1, -55), 511.99, 512, -512, - 709.78, -745.13, 1000, math.Inf(1), math.Inf(-1), math.NaN()} - for i := 0; i < 20000; i++ { - xs = append(xs, rng.Float64()*1460-750, (rng.Float64()*2-1)*math.Ldexp(1, -rng.Intn(70))) - } - got := make([]float64, len(xs)) - Exp(got, xs) - for i, x := range xs { - if w := exp1(x); math.Float64bits(got[i]) != math.Float64bits(w) && !(math.IsNaN(w) && math.IsNaN(got[i])) { - t.Fatalf("Exp(%v) = %v, exp1 = %v", x, got[i], w) - } - } -} - -func TestExpSpecial(t *testing.T) { - for _, c := range []struct{ x, want float64 }{ - {0, 1}, {math.Copysign(0, -1), 1}, {1e-300, 1}, {-1e-20, 1}, - {math.Inf(1), math.Inf(1)}, {math.Inf(-1), 0}, - {709.8, math.Inf(1)}, {1000, math.Inf(1)}, {2000, math.Inf(1)}, - {-745.2, 0}, {-1000, 0}, {-2000, 0}, - {1, math.E}, - } { - if got := exp1(c.x); got != c.want { - t.Errorf("exp(%v) = %v, want %v", c.x, got, c.want) - } - } - // Against the reference: the subnormal range (just above ln(2^-1075) exp - // rounds up to the smallest subnormal; arm64's math.Exp gives 0), and the - // top of the finite range, where amd64's math.Exp (golang/go#81995) returns +Inf from - // x >= 1023.5*ln2 ~ 709.436 (its k = round(x*log2e) reaches 1024 and the - // biased exponent check calls that overflow, though fr*2^1024 is finite). - for _, x := range []float64{-745.13, -745.1332, -740, -720.5, -708.5, -708.39, - 709.436, 709.5, 709.78, 709.782712893383} { - want, _ := refExp(x).Float64() - if got := exp1(x); got != want { - t.Errorf("exp(%v) = %v, want %v", x, got, want) - } - } - if got := exp1(math.NaN()); !math.IsNaN(got) { - t.Errorf("exp(NaN) = %v", got) - } - dst := make([]float64, 3) - Exp(dst, []float64{0, 1, math.Inf(-1)}) - eqSlice(t, dst, []float64{1, math.E, 0}) -} - -// TestExpP: the parallel form equals the serial one element for element. -func TestExpP(t *testing.T) { - src := randVec(5000, 3) - for i := range src { - src[i] *= 40 - } - want, got := make([]float64, len(src)), make([]float64, len(src)) - Exp(want, src) - withThresholds(64, 1<<14, func() { ExpP(got, src) }) - eqSlice(t, got, want) - ExpP(got[:10], src[:10]) - eqSlice(t, got[:10], want[:10]) -} diff --git a/internal/kernels/expcr.go b/internal/kernels/expcr.go deleted file mode 100644 index dee70cf..0000000 --- a/internal/kernels/expcr.go +++ /dev/null @@ -1,94 +0,0 @@ -package kernels - -import "math" - -// expCR is the double-precision exp of Arm's correctly rounded vector expf -// (optimized-routines, math/aarch64/advsimd/cr_expf.c and cr_expf_utils.h, -// MIT), transcribed one operation per line so that a SIMD kernel can repeat -// it bit for bit: every vfmaq/vfmsq is math.FMA. -// -// x/ln2 is rounded to a multiple of 1/64 by adding a shift whose ULP is 2^-6; -// the low 6 bits of the sum index a table of 2^(i/64), and the bits above -// (with the bias 0x3ff folded into the shift) are the exponent of the scale. -// exp(r) - 1 on |r| <= ln2/128 is a degree-5 polynomial. Valid for -// |x| <= 512, which covers every float32 whose exp is neither 0 nor +Inf. -func expCR(x float64) float64 { - z := math.FMA(x, crInvLn2, crShift) - n := z - crShift - r := math.FMA(-n, crLn2Hi, x) - r = math.FMA(-n, crLn2Lo, r) - r2 := r * r - p01 := math.FMA(r, crC1, crC0) - p23 := math.FMA(r, crC3, crC2) - p04 := math.FMA(r2, p23, p01) - y := math.FMA(r2, p04, r) - u := math.Float64bits(z) - s := math.Float64frombits(u<<46&0xfff0000000000000 | crMant[u&63]) - return math.FMA(s, y, s) -} - -const ( - crShift = 0x1.800000000ffc0p+46 - crInvLn2 = 0x1.71547652b82fep+0 - crLn2Hi = 0x1.62e42fefa39efp-1 - crLn2Lo = 0x1.abc9e3b39803fp-56 - crC0 = 0x1.fffffffffdbcep-2 - crC1 = 0x1.55555555543c2p-3 - crC2 = 0x1.555573c64f2e3p-5 - crC3 = 0x1.111126b4eff73p-7 -) - -// crMant are the mantissas of 2^(i/64), i = 0..63 (cr_expf_utils.h). -var crMant = [64]uint64{ - 0x0000000000000, 0x02c9a3e778061, 0x059b0d3158574, 0x0874518759bc8, - 0x0b5586cf9890f, 0x0e3ec32d3d1a2, 0x11301d0125b51, 0x1429aaea92de0, - 0x172b83c7d517b, 0x1a35beb6fcb75, 0x1d4873168b9aa, 0x2063b88628cd6, - 0x2387a6e756238, 0x26b4565e27cdd, 0x29e9df51fdee1, 0x2d285a6e4030b, - 0x306fe0a31b715, 0x33c08b26416ff, 0x371a7373aa9cb, 0x3a7db34e59ff7, - 0x3dea64c123422, 0x4160a21f72e2a, 0x44e086061892d, 0x486a2b5c13cd0, - 0x4bfdad5362a27, 0x4f9b2769d2ca7, 0x5342b569d4f82, 0x56f4736b527da, - 0x5ab07dd485429, 0x5e76f15ad2148, 0x6247eb03a5585, 0x6623882552225, - 0x6a09e667f3bcd, 0x6dfb23c651a2f, 0x71f75e8ec5f74, 0x75feb564267c9, - 0x7a11473eb0187, 0x7e2f336cf4e62, 0x82589994cce13, 0x868d99b4492ed, - 0x8ace5422aa0db, 0x8f1ae99157736, 0x93737b0cdc5e5, 0x97d829fde4e50, - 0x9c49182a3f090, 0xa0c667b5de565, 0xa5503b23e255d, 0xa9e6b5579fdbf, - 0xae89f995ad3ad, 0xb33a2b84f15fb, 0xb7f76f2fb5e47, 0xbcc1e904bc1d2, - 0xc199bdd85529c, 0xc67f12e57d14b, 0xcb720dcef9069, 0xd072d4a07897c, - 0xd5818dcfba487, 0xda9e603db3285, 0xdfc97337b9b5f, 0xe502ee78b3ff6, - 0xea4afa2a490da, 0xefa1bee615a27, 0xf50765b6e4540, 0xfa7c1819e90d8, -} - -// expCRSlice writes expCR(src[i]) into dst[i] for |src[i]| <= 512, and the -// limits outside: +Inf above, 0 below, NaN for NaN (the float32 results -// there are the same). -func expCRSlice(dst, src []float64) { - dst = dst[:len(src)] - for i, x := range src { - switch { - case math.Abs(x) <= 512: - dst[i] = expCR(x) - case x > 0: - dst[i] = math.Inf(1) - case x < 0: - dst[i] = 0 - default: - dst[i] = x // NaN - } - } -} - -// exp32CR is float32 exp on expCR, rounded once. -func exp32CR(dst, src []float32) { - src = src[:len(dst)] - var buf [via64Chunk]float64 - for i := 0; i < len(dst); i += via64Chunk { - b := buf[:min(via64Chunk, len(dst)-i)] - for j := range b { - b[j] = float64(src[i+j]) - } - expCRSlice(b, b) - for j, v := range b { - dst[i+j] = float32(v) - } - } -} diff --git a/internal/kernels/expcr_amd64.go b/internal/kernels/expcr_amd64.go deleted file mode 100644 index 66f8342..0000000 --- a/internal/kernels/expcr_amd64.go +++ /dev/null @@ -1,50 +0,0 @@ -package kernels - -// expCRConsts is the constant block exp32CRAVX2 reads (asmgen/amd64/gen.go, -// exp32CRKernel): each constant as a 4-lane vector at the offset the kernel -// uses, then the clamp bounds as 8 floats, then the 64 table mantissas. -type expCRConsts struct { - shift, invLn2, ln2Hi, ln2Lo [4]float64 // offsets 0, 32, 64, 96 - c0, c1, c2, c3 [4]float64 // 128, 160, 192, 224 - idxMask, expMask [4]uint64 // 256, 288 - lo, hi [8]float32 // 320, 352 - table [64]uint64 // 384 -} - -var expCRK = func() *expCRConsts { - k := new(expCRConsts) - for i := range 4 { - k.shift[i], k.invLn2[i], k.ln2Hi[i], k.ln2Lo[i] = crShift, crInvLn2, crLn2Hi, crLn2Lo - k.c0[i], k.c1[i], k.c2[i], k.c3[i] = crC0, crC1, crC2, crC3 - k.idxMask[i], k.expMask[i] = 63, 0xfff0000000000000 - } - for i := range 8 { - k.lo[i], k.hi[i] = -512, 512 - } - k.table = crMant - return k -}() - -// useExp32AVX2 selects the AVX2 kernel; a var so tests can force the Go path. -var useExp32AVX2 = hasAVX2() && useFMA - -// exp32Fast is the correctly rounded float32 exp: blocks of 8 on the AVX2 -// kernel, the rest (and every element without AVX2+FMA) on exp32CR. Both -// compute expCR's operations in the same order, so they agree bit for bit. -func exp32Fast(dst, src []float32) { - src = src[:len(dst)] - n := 0 - if useExp32AVX2 { - n = len(dst) &^ 7 - if n > 0 { - exp32CRAVX2(&dst[0], &src[0], n, expCRK) - } - } - exp32CR(dst[n:], src[n:]) -} - -// hasAVX2 reports CPU and OS support for AVX2 (go-asmgen FeatureProbe). -func hasAVX2() bool - -//go:noescape -func exp32CRAVX2(dst, src *float32, n int, k *expCRConsts) diff --git a/internal/kernels/expcr_amd64_test.go b/internal/kernels/expcr_amd64_test.go deleted file mode 100644 index bf8c2e3..0000000 --- a/internal/kernels/expcr_amd64_test.go +++ /dev/null @@ -1,45 +0,0 @@ -package kernels - -import ( - "math" - "math/rand" - "testing" -) - -// TestExp32AVX2MatchesGo holds the AVX2 kernel to exp32CR bit for bit (NaN -// to NaN) on random inputs across the float32 range, the special values and -// the clamp and overflow edges, at every length 0..70 (blocks of 8 and the Go -// tail). -func TestExp32AVX2MatchesGo(t *testing.T) { - if !useExp32AVX2 { - t.Skip("no AVX2+FMA here") - } - r := rand.New(rand.NewSource(3)) - src := make([]float32, 0, 100000) - for _, v := range []float32{0, float32(math.Copysign(0, -1)), 1, -1, 88.72283, 88.72284, 89, -103.972, -103.973, -104, - 511.9, 512, 512.1, -512, -512.1, 1e30, -1e30, float32(math.Inf(1)), float32(math.Inf(-1)), - float32(math.NaN()), math.Float32frombits(0xffc00000), math.SmallestNonzeroFloat32} { - src = append(src, v) - } - for len(src) < cap(src) { - src = append(src, math.Float32frombits(r.Uint32())) - } - want := make([]float32, len(src)) - exp32CR(want, src) - got := make([]float32, len(src)) - exp32Fast(got, src) - for i := range src { - if math.Float32bits(got[i]) != math.Float32bits(want[i]) && !(got[i] != got[i] && want[i] != want[i]) { - t.Fatalf("exp(0x%08x): AVX2 0x%08x, Go 0x%08x", math.Float32bits(src[i]), math.Float32bits(got[i]), math.Float32bits(want[i])) - } - } - for n := 0; n <= 70; n++ { - d := make([]float32, n) - exp32Fast(d, src[:n]) - for i := range d { - if math.Float32bits(d[i]) != math.Float32bits(want[i]) && !(d[i] != d[i] && want[i] != want[i]) { - t.Fatalf("length %d [%d]: 0x%08x, want 0x%08x", n, i, math.Float32bits(d[i]), math.Float32bits(want[i])) - } - } - } -} diff --git a/internal/kernels/expcr_fence_amd64_test.go b/internal/kernels/expcr_fence_amd64_test.go deleted file mode 100644 index 6d6b959..0000000 --- a/internal/kernels/expcr_fence_amd64_test.go +++ /dev/null @@ -1,20 +0,0 @@ -//go:build linux || darwin - -package kernels - -import "testing" - -// TestExp32AVX2StaysInside fences src and dst of the AVX2 exp kernel -// (fences32 needs mmap/mprotect, hence the build tag of fence_test.go). -func TestExp32AVX2StaysInside(t *testing.T) { - if !useExp32AVX2 { - t.Skip("no AVX2+FMA here") - } - for n := 8; n <= 64; n += 8 { - for _, s := range fences32(t, n) { - for _, d := range fences32(t, n) { - exp32CRAVX2(&d[0], &s[0], n, expCRK) - } - } - } -} diff --git a/internal/kernels/expcr_other.go b/internal/kernels/expcr_other.go deleted file mode 100644 index 80b6267..0000000 --- a/internal/kernels/expcr_other.go +++ /dev/null @@ -1,10 +0,0 @@ -//go:build !amd64 - -package kernels - -// exp32Fast is the correctly rounded float32 exp: expCR in Go on these -// architectures (no SIMD kernel yet). -func exp32Fast(dst, src []float32) { exp32CR(dst, src) } - -// log32Fast is the correctly rounded float32 log: logCR in Go here. -func log32Fast(dst, src []float32) { logCRSlice(dst, src) } diff --git a/internal/kernels/expcr_test.go b/internal/kernels/expcr_test.go deleted file mode 100644 index dea44ae..0000000 --- a/internal/kernels/expcr_test.go +++ /dev/null @@ -1,53 +0,0 @@ -package kernels - -import ( - "math" - "math/big" - "math/rand" - "testing" -) - -// TestExpCRULP bounds expCR's float64 error against the 300-bit reference -// over the inputs float32 exp needs (|x| <= 104, the float32 overflow and -// underflow edges included). Measured worst: 1.006 ULP. The exhaustive -// float32 test certifies an expCR result only when it is more than 4 float64 -// ULP from a rounding midpoint, so this bound must stay below 2 (half the -// margin). -func TestExpCRULP(t *testing.T) { - r := rand.New(rand.NewSource(11)) - var worst float64 - var at float64 - for i := 0; i < 200000; i++ { - x := float64(float32(r.Float64()*193 - 104)) // a float32 input, widened - if e := ulpErr(expCR(x), refExp(x)); e > worst { - worst, at = e, x - } - } - t.Logf("expCR: worst %.3f float64 ULP at x = %v", worst, at) - if worst >= 2 { - t.Errorf("expCR error %.3f ULP at %v: the certifier's 4-ULP margin needs it below 2", worst, at) - } -} - -// TestExp32CRCorrectlyRounded: the prototype, through the same exhaustive -// certifier as Exp32 (NDARRAY_EXHAUSTIVE=1 for all 2^32 inputs). -func TestExp32CRCorrectlyRounded(t *testing.T) { - runExhaustive(t, exhaustiveCase{"expCR", exp32CR, expCRSlice, refExp, exp32Special, 4}, exhaustiveStride()) -} - -// TestExpCRTable re-derives the 64 mantissas of 2^(i/64) from the 300-bit -// reference: the table was copied, not trusted. -func TestExpCRTable(t *testing.T) { - root := new(big.Float).SetPrec(refPrec).SetInt64(2) - for range 6 { // 2^(1/64) - root.Sqrt(root) - } - e := new(big.Float).SetPrec(refPrec).SetInt64(1) - for i, m := range crMant { - want, _ := e.Float64() - if got := math.Float64frombits(0x3ff0000000000000 | m); got != want { - t.Errorf("crMant[%d] gives %v, the nearest double to 2^(%d/64) is %v", i, got, i, want) - } - e.Mul(e, root) - } -} diff --git a/internal/kernels/explog.go b/internal/kernels/explog.go new file mode 100644 index 0000000..e8f54b4 --- /dev/null +++ b/internal/kernels/explog.go @@ -0,0 +1,83 @@ +package kernels + +import ( + "math" + + "github.com/go-simd/floats" +) + +// Exp, Log and their float32 forms are go-simd/floats': K1 of the go-scikit +// performance campaign moved the kernels there (they were written here: +// go-ndarray#94 and #96 for AVX2, #95 for NEON), so that a kernel lives once. +// What stays here is the parallel split above ParThreshold and Log10. +// +// - floats.Exp / floats.Log: float64, 0.511 / 0.519 ULP (Arm optimized-routines); +// - floats.Float32Exp / floats.Float32Log: correctly rounded float32, proven +// over every float32 input in go-simd/floats (TestExp32CorrectlyRounded, +// TestLog32CorrectlyRounded), on AVX2 (amd64) and NEON (arm64) kernels. + +// Exp writes exp(src[i]) into dst[i]. +func Exp(dst, src []float64) { floats.Exp(dst, src) } + +// Log writes log(src[i]) into dst[i]. +func Log(dst, src []float64) { floats.Log(dst, src) } + +// Exp32 writes the correctly rounded exp(src[i]) into dst[i]. +func Exp32(dst, src []float32) { floats.Float32Exp(dst, src) } + +// Log32 writes the correctly rounded log(src[i]) into dst[i]. +func Log32(dst, src []float32) { floats.Float32Log(dst, src) } + +// Log10 writes log10(src[i]) into dst[i] as log(x)*(1/ln10), the formula of +// math.Log10, on floats.Log. +func Log10(dst, src []float64) { + Log(dst, src) + for i := range dst { + dst[i] *= 1 / math.Ln10 + } +} + +// ExpP is Exp, split across cores above ParThreshold. +func ExpP(dst, src []float64) { + n := len(dst) + if n < ParThreshold { + Exp(dst, src) + return + } + parallelFor(n, numWorkers(n), func(lo, hi int) { Exp(dst[lo:hi], src[lo:hi]) }) +} + +// LogP is Log, split across cores above ParThreshold. +func LogP(dst, src []float64) { + n := len(dst) + if n < ParThreshold { + Log(dst, src) + return + } + parallelFor(n, numWorkers(n), func(lo, hi int) { Log(dst[lo:hi], src[lo:hi]) }) +} + +// Log10P is Log10, split across cores above ParThreshold. +func Log10P(dst, src []float64) { + n := len(dst) + if n < ParThreshold { + Log10(dst, src) + return + } + parallelFor(n, numWorkers(n), func(lo, hi int) { Log10(dst[lo:hi], src[lo:hi]) }) +} + +// Exp32P is Exp32, split across cores above ParThreshold. +func Exp32P(dst, src []float32) { split32(Exp32, dst, src) } + +// Log32P is Log32, split across cores above ParThreshold. +func Log32P(dst, src []float32) { split32(Log32, dst, src) } + +func split32(f func(dst, src []float32), dst, src []float32) { + n := len(dst) + if n < ParThreshold { + f(dst, src) + return + } + parallelFor(n, numWorkers(n), func(lo, hi int) { f(dst[lo:hi], src[lo:hi]) }) +} diff --git a/internal/kernels/generate.go b/internal/kernels/generate.go new file mode 100644 index 0000000..80fe3b8 --- /dev/null +++ b/internal/kernels/generate.go @@ -0,0 +1,10 @@ +package kernels + +// The amd64 and arm64 assembly of this package is generated by go-asmgen. The +// float64 sum, sqrt, max/min and elementwise kernels, and the exp/log kernels, +// live in go-simd/floats since K1 of the go-scikit performance campaign; what is +// generated here is GEMM, the float64 dot, and the float32 elementwise, sqrt, +// extremum and pairwise-block kernels. +// +//go:generate sh -c "cd asmgen/amd64 && GOFLAGS=-mod=mod go run gen.go && mv sum_amd64.s ../../sum_amd64.s" +//go:generate sh -c "cd asmgen/arm64 && GOFLAGS=-mod=mod go run gen.go && mv sum_arm64.s ../../sum_arm64.s" diff --git a/internal/kernels/kernels.go b/internal/kernels/kernels.go index cb54fee..67dc981 100644 --- a/internal/kernels/kernels.go +++ b/internal/kernels/kernels.go @@ -11,7 +11,11 @@ // a kernel it is already materialised in row-major order. package kernels -import "math" +import ( + "math" + + "github.com/go-simd/floats" +) // Add writes a[i]+b[i] into dst[i] for every element. func Add(dst, a, b []float64) { @@ -326,15 +330,7 @@ func CumProdAxis(dst, src []float64, outer, axisLen, inner int) { // numpy.clip. A NaN bound or element follows Go's min/max comparison order; // callers pass lo <= hi. func Clip(dst, src []float64, lo, hi float64) { - for i, v := range src { - if v < lo { - v = lo - } - if v > hi { - v = hi - } - dst[i] = v - } + floats.ClampScalar(dst, src, lo, hi) } // Where writes into dst the value t[i] where cond[i] is non-zero (truthy), diff --git a/internal/kernels/log.go b/internal/kernels/log.go deleted file mode 100644 index c3a5ae9..0000000 --- a/internal/kernels/log.go +++ /dev/null @@ -1,305 +0,0 @@ -package kernels - -import "math" - -// Log is a port of the double-precision log of Arm's optimized-routines -// (math/log.c and math/log_data.c, by Szabolcs Nagy; also glibc's log since -// 2.28), used under its MIT licence: -// -// Copyright (c) 2018-2025, Arm Limited. -// SPDX-License-Identifier: MIT OR Apache-2.0 WITH LLVM-exception -// -// x = 2^k z with z in [0x1.6p-1, 0x1.6p0), split into N = 128 subintervals; -// log(x) = k*ln2 + log(c) + log1p(z/c - 1) for c near the centre of z's -// subinterval, with 1/c and log(c) from a table and log1p from a degree-6 -// polynomial. Inputs within 2^-4 of 1 take a separate degree-12 polynomial. -// r = z/c - 1 is computed with one fused multiply-add (math.FMA, exact on every -// target), the "fast fma" variant of log.c, so the second table its non-FMA -// variant needs is not used. Arm documents 0.519 ULP (0.507 near 1); -// TestLogULP measures this port against a 300-bit reference. - -const ( - logTableBits = 7 - logN = 1 << logTableBits - logOff = 0x3fe6000000000000 - logLn2hi = 0x1.62e42fefa3800p-1 - logLn2lo = 0x1.ef35793c76730p-45 - logA0 = -0x1.0000000000001p-1 - logA1 = 0x1.555555551305bp-2 - logA2 = -0x1.fffffffeb459p-3 - logA3 = 0x1.999b324f10111p-3 - logA4 = -0x1.55575e506c89fp-3 - logB0 = -0x1p-1 - logB1 = 0x1.5555555555577p-2 - logB2 = -0x1.ffffffffffdcbp-3 - logB3 = 0x1.999999995dd0cp-3 - logB4 = -0x1.55555556745a7p-3 - logB5 = 0x1.24924a344de3p-3 - logB6 = -0x1.fffffa4423d65p-4 - logB7 = 0x1.c7184282ad6cap-4 - logB8 = -0x1.999eb43b068ffp-4 - logB9 = 0x1.78182f7afd085p-4 - logB10 = -0x1.5521375d145cdp-4 -) - -// logTab holds {1/c, log(c)} per subinterval, copied from log_data.c -// (N == 128); TestLogTable checks every entry against its defining properties. -// logEntry is one row of logTab: 1/c and log(c), 16 bytes (logAVX2 gathers -// invc at 16i and logc at 16i + 8). -type logEntry struct{ invc, logc float64 } - -var logTab = [logN]logEntry{ - {0x1.734f0c3e0de9fp+0, -0x1.7cc7f79e69000p-2}, - {0x1.713786a2ce91fp+0, -0x1.76feec20d0000p-2}, - {0x1.6f26008fab5a0p+0, -0x1.713e31351e000p-2}, - {0x1.6d1a61f138c7dp+0, -0x1.6b85b38287800p-2}, - {0x1.6b1490bc5b4d1p+0, -0x1.65d5590807800p-2}, - {0x1.69147332f0cbap+0, -0x1.602d076180000p-2}, - {0x1.6719f18224223p+0, -0x1.5a8ca86909000p-2}, - {0x1.6524f99a51ed9p+0, -0x1.54f4356035000p-2}, - {0x1.63356aa8f24c4p+0, -0x1.4f637c36b4000p-2}, - {0x1.614b36b9ddc14p+0, -0x1.49da7fda85000p-2}, - {0x1.5f66452c65c4cp+0, -0x1.445923989a800p-2}, - {0x1.5d867b5912c4fp+0, -0x1.3edf439b0b800p-2}, - {0x1.5babccb5b90dep+0, -0x1.396ce448f7000p-2}, - {0x1.59d61f2d91a78p+0, -0x1.3401e17bda000p-2}, - {0x1.5805612465687p+0, -0x1.2e9e2ef468000p-2}, - {0x1.56397cee76bd3p+0, -0x1.2941b3830e000p-2}, - {0x1.54725e2a77f93p+0, -0x1.23ec58cda8800p-2}, - {0x1.52aff42064583p+0, -0x1.1e9e129279000p-2}, - {0x1.50f22dbb2bddfp+0, -0x1.1956d2b48f800p-2}, - {0x1.4f38f4734ded7p+0, -0x1.141679ab9f800p-2}, - {0x1.4d843cfde2840p+0, -0x1.0edd094ef9800p-2}, - {0x1.4bd3ec078a3c8p+0, -0x1.09aa518db1000p-2}, - {0x1.4a27fc3e0258ap+0, -0x1.047e65263b800p-2}, - {0x1.4880524d48434p+0, -0x1.feb224586f000p-3}, - {0x1.46dce1b192d0bp+0, -0x1.f474a7517b000p-3}, - {0x1.453d9d3391854p+0, -0x1.ea4443d103000p-3}, - {0x1.43a2744b4845ap+0, -0x1.e020d44e9b000p-3}, - {0x1.420b54115f8fbp+0, -0x1.d60a22977f000p-3}, - {0x1.40782da3ef4b1p+0, -0x1.cc00104959000p-3}, - {0x1.3ee8f5d57fe8fp+0, -0x1.c202956891000p-3}, - {0x1.3d5d9a00b4ce9p+0, -0x1.b81178d811000p-3}, - {0x1.3bd60c010c12bp+0, -0x1.ae2c9ccd3d000p-3}, - {0x1.3a5242b75dab8p+0, -0x1.a45402e129000p-3}, - {0x1.38d22cd9fd002p+0, -0x1.9a877681df000p-3}, - {0x1.3755bc5847a1cp+0, -0x1.90c6d69483000p-3}, - {0x1.35dce49ad36e2p+0, -0x1.87120a645c000p-3}, - {0x1.34679984dd440p+0, -0x1.7d68fb4143000p-3}, - {0x1.32f5cceffcb24p+0, -0x1.73cb83c627000p-3}, - {0x1.3187775a10d49p+0, -0x1.6a39a9b376000p-3}, - {0x1.301c8373e3990p+0, -0x1.60b3154b7a000p-3}, - {0x1.2eb4ebb95f841p+0, -0x1.5737d76243000p-3}, - {0x1.2d50a0219a9d1p+0, -0x1.4dc7b8fc23000p-3}, - {0x1.2bef9a8b7fd2ap+0, -0x1.4462c51d20000p-3}, - {0x1.2a91c7a0c1babp+0, -0x1.3b08abc830000p-3}, - {0x1.293726014b530p+0, -0x1.31b996b490000p-3}, - {0x1.27dfa5757a1f5p+0, -0x1.2875490a44000p-3}, - {0x1.268b39b1d3bbfp+0, -0x1.1f3b9f879a000p-3}, - {0x1.2539d838ff5bdp+0, -0x1.160c8252ca000p-3}, - {0x1.23eb7aac9083bp+0, -0x1.0ce7f57f72000p-3}, - {0x1.22a012ba940b6p+0, -0x1.03cdc49fea000p-3}, - {0x1.2157996cc4132p+0, -0x1.f57bdbc4b8000p-4}, - {0x1.201201dd2fc9bp+0, -0x1.e370896404000p-4}, - {0x1.1ecf4494d480bp+0, -0x1.d17983ef94000p-4}, - {0x1.1d8f5528f6569p+0, -0x1.bf9674ed8a000p-4}, - {0x1.1c52311577e7cp+0, -0x1.adc79202f6000p-4}, - {0x1.1b17c74cb26e9p+0, -0x1.9c0c3e7288000p-4}, - {0x1.19e010c2c1ab6p+0, -0x1.8a646b372c000p-4}, - {0x1.18ab07bb670bdp+0, -0x1.78d01b3ac0000p-4}, - {0x1.1778a25efbcb6p+0, -0x1.674f145380000p-4}, - {0x1.1648d354c31dap+0, -0x1.55e0e6d878000p-4}, - {0x1.151b990275fddp+0, -0x1.4485cdea1e000p-4}, - {0x1.13f0ea432d24cp+0, -0x1.333d94d6aa000p-4}, - {0x1.12c8b7210f9dap+0, -0x1.22079f8c56000p-4}, - {0x1.11a3028ecb531p+0, -0x1.10e4698622000p-4}, - {0x1.107fbda8434afp+0, -0x1.ffa6c6ad20000p-5}, - {0x1.0f5ee0f4e6bb3p+0, -0x1.dda8d4a774000p-5}, - {0x1.0e4065d2a9fcep+0, -0x1.bbcece4850000p-5}, - {0x1.0d244632ca521p+0, -0x1.9a1894012c000p-5}, - {0x1.0c0a77ce2981ap+0, -0x1.788583302c000p-5}, - {0x1.0af2f83c636d1p+0, -0x1.5715e67d68000p-5}, - {0x1.09ddb98a01339p+0, -0x1.35c8a49658000p-5}, - {0x1.08cabaf52e7dfp+0, -0x1.149e364154000p-5}, - {0x1.07b9f2f4e28fbp+0, -0x1.e72c082eb8000p-6}, - {0x1.06ab58c358f19p+0, -0x1.a55f152528000p-6}, - {0x1.059eea5ecf92cp+0, -0x1.63d62cf818000p-6}, - {0x1.04949cdd12c90p+0, -0x1.228fb8caa0000p-6}, - {0x1.038c6c6f0ada9p+0, -0x1.c317b20f90000p-7}, - {0x1.02865137932a9p+0, -0x1.419355daa0000p-7}, - {0x1.0182427ea7348p+0, -0x1.81203c2ec0000p-8}, - {0x1.008040614b195p+0, -0x1.0040979240000p-9}, - {0x1.fe01ff726fa1ap-1, 0x1.feff384900000p-9}, - {0x1.fa11cc261ea74p-1, 0x1.7dc41353d0000p-7}, - {0x1.f6310b081992ep-1, 0x1.3cea3c4c28000p-6}, - {0x1.f25f63ceeadcdp-1, 0x1.b9fc114890000p-6}, - {0x1.ee9c8039113e7p-1, 0x1.1b0d8ce110000p-5}, - {0x1.eae8078cbb1abp-1, 0x1.58a5bd001c000p-5}, - {0x1.e741aa29d0c9bp-1, 0x1.95c8340d88000p-5}, - {0x1.e3a91830a99b5p-1, 0x1.d276aef578000p-5}, - {0x1.e01e009609a56p-1, 0x1.07598e598c000p-4}, - {0x1.dca01e577bb98p-1, 0x1.253f5e30d2000p-4}, - {0x1.d92f20b7c9103p-1, 0x1.42edd8b380000p-4}, - {0x1.d5cac66fb5ccep-1, 0x1.606598757c000p-4}, - {0x1.d272caa5ede9dp-1, 0x1.7da76356a0000p-4}, - {0x1.cf26e3e6b2ccdp-1, 0x1.9ab434e1c6000p-4}, - {0x1.cbe6da2a77902p-1, 0x1.b78c7bb0d6000p-4}, - {0x1.c8b266d37086dp-1, 0x1.d431332e72000p-4}, - {0x1.c5894bd5d5804p-1, 0x1.f0a3171de6000p-4}, - {0x1.c26b533bb9f8cp-1, 0x1.067152b914000p-3}, - {0x1.bf583eeece73fp-1, 0x1.147858292b000p-3}, - {0x1.bc4fd75db96c1p-1, 0x1.2266ecdca3000p-3}, - {0x1.b951e0c864a28p-1, 0x1.303d7a6c55000p-3}, - {0x1.b65e2c5ef3e2cp-1, 0x1.3dfc33c331000p-3}, - {0x1.b374867c9888bp-1, 0x1.4ba366b7a8000p-3}, - {0x1.b094b211d304ap-1, 0x1.5933928d1f000p-3}, - {0x1.adbe885f2ef7ep-1, 0x1.66acd2418f000p-3}, - {0x1.aaf1d31603da2p-1, 0x1.740f8ec669000p-3}, - {0x1.a82e63fd358a7p-1, 0x1.815c0f51af000p-3}, - {0x1.a5740ef09738bp-1, 0x1.8e92954f68000p-3}, - {0x1.a2c2a90ab4b27p-1, 0x1.9bb3602f84000p-3}, - {0x1.a01a01393f2d1p-1, 0x1.a8bed1c2c0000p-3}, - {0x1.9d79f24db3c1bp-1, 0x1.b5b515c01d000p-3}, - {0x1.9ae2505c7b190p-1, 0x1.c2967ccbcc000p-3}, - {0x1.9852ef297ce2fp-1, 0x1.cf635d5486000p-3}, - {0x1.95cbaeea44b75p-1, 0x1.dc1bd3446c000p-3}, - {0x1.934c69de74838p-1, 0x1.e8c01b8cfe000p-3}, - {0x1.90d4f2f6752e6p-1, 0x1.f5509c0179000p-3}, - {0x1.8e6528effd79dp-1, 0x1.00e6c121fb800p-2}, - {0x1.8bfce9fcc007cp-1, 0x1.071b80e93d000p-2}, - {0x1.899c0dabec30ep-1, 0x1.0d46b9e867000p-2}, - {0x1.87427aa2317fbp-1, 0x1.13687334bd000p-2}, - {0x1.84f00acb39a08p-1, 0x1.1980d67234800p-2}, - {0x1.82a49e8653e55p-1, 0x1.1f8ffe0cc8000p-2}, - {0x1.8060195f40260p-1, 0x1.2595fd7636800p-2}, - {0x1.7e22563e0a329p-1, 0x1.2b9300914a800p-2}, - {0x1.7beb377dcb5adp-1, 0x1.3187210436000p-2}, - {0x1.79baa679725c2p-1, 0x1.377266dec1800p-2}, - {0x1.77907f2170657p-1, 0x1.3d54ffbaf3000p-2}, - {0x1.756cadbd6130cp-1, 0x1.432eee32fe000p-2}, -} - -// logNear1 is log(x) for x within 2^-4 of 1, where the table path would lose -// accuracy: a degree-12 polynomial with r*r split into exact halves. The -// conversions to float64 keep the compiler from fusing the compensated sums. -func logNear1(x float64) float64 { - if x == 1 { - return 0 - } - r := x - 1.0 - r2 := r * r - r3 := r * r2 - y := r3 * (logB1 + r*logB2 + r2*logB3 + - r3*(logB4+r*logB5+r2*logB6+r3*(logB7+r*logB8+r2*logB9+r3*logB10))) - w := float64(r * 0x1p27) - rhi := float64(r+w) - w - rlo := r - rhi - w = float64(float64(rhi*rhi) * logB0) // exact: rhi has 26 bits, B0 = -0.5 - hi := r + w - lo := r - hi + w - lo += float64(float64(logB0*rlo) * float64(rhi+r)) - y += lo - y += hi - return y -} - -// log1 is log(x) for one element, every case (log.c's log). -func log1(x float64) float64 { - ix := math.Float64bits(x) - if ix-math.Float64bits(1.0-0x1p-4) < math.Float64bits(1.0+0x1.09p-4)-math.Float64bits(1.0-0x1p-4) { - return logNear1(x) - } - if top := uint32(ix >> 48); top-0x0010 >= 0x7ff0-0x0010 { - // x < 0x1p-1022, or Inf or NaN. - switch { - case ix*2 == 0: - return math.Inf(-1) // log(±0) - case ix == math.Float64bits(math.Inf(1)): - return x - case top&0x8000 != 0 || top&0x7ff0 == 0x7ff0: - return math.NaN() // x < 0, or NaN - } - // Subnormal: normalise it. - ix = math.Float64bits(x * 0x1p52) - ix -= 52 << 52 - } - return logReduced(ix) -} - -// logReduced is the table path for a positive normal x given by its bits. -func logReduced(ix uint64) float64 { - tmp := ix - logOff - i := (tmp >> (52 - logTableBits)) % logN - k := int64(tmp) >> 52 - iz := ix - tmp&(0xfff<<52) - e := &logTab[i] - r := math.FMA(math.Float64frombits(iz), e.invc, -1.0) - kd := float64(k) - w := kd*logLn2hi + e.logc // exact (table property) - hi := w + r - lo := w - hi + r + kd*logLn2lo - r2 := r * r - return lo + r2*logA0 + r*r2*(logA1+r*logA2+r2*(logA3+r*logA4)) + hi -} - -// Log writes log(src[i]) into dst[i]: on the AVX2 kernel where there is one -// (log_amd64.go), on logScalar otherwise and for the blocks the kernel hands -// back. -func Log(dst, src []float64) { logFast(dst[:len(src)], src) } - -// logScalar is Log in Go. Positive normal inputs away from 1 take the table -// path straight from the loop; everything else goes through log1. -// TestLogLoopMatchesLog1 holds the two to the same bits. -func logScalar(dst, src []float64) { - dst = dst[:len(src)] - const lo, span = 0x0010 << 48, (0x7ff0 - 0x0010) << 48 - near1lo := math.Float64bits(1.0 - 0x1p-4) - near1span := math.Float64bits(1.0+0x1.09p-4) - near1lo - for i, x := range src { - ix := math.Float64bits(x) - if ix-lo >= span || ix-near1lo < near1span { - dst[i] = log1(x) - continue - } - // logReduced, written out: it does not inline, and a call per - // element is a large share of the time (as for Exp). - tmp := ix - logOff - e := &logTab[(tmp>>(52-logTableBits))%logN] - r := math.FMA(math.Float64frombits(ix-tmp&(0xfff<<52)), e.invc, -1.0) - kd := float64(int64(tmp) >> 52) - w := kd*logLn2hi + e.logc - hi := w + r - lo := w - hi + r + kd*logLn2lo - r2 := r * r - dst[i] = lo + r2*logA0 + r*r2*(logA1+r*logA2+r2*(logA3+r*logA4)) + hi - } -} - -// Log10 writes log10(src[i]) into dst[i] as log(x)*(1/ln10), the formula of -// math.Log10, on this package's log. -func Log10(dst, src []float64) { - Log(dst, src) - for i := range dst { - dst[i] *= 1 / math.Ln10 - } -} - -// Log10P is Log10, split across cores above ParThreshold. -func Log10P(dst, src []float64) { - n := len(dst) - if n < ParThreshold { - Log10(dst, src) - return - } - parallelFor(n, numWorkers(n), func(lo, hi int) { Log10(dst[lo:hi], src[lo:hi]) }) -} - -// LogP is Log, split across cores above ParThreshold. -func LogP(dst, src []float64) { - n := len(dst) - if n < ParThreshold { - Log(dst, src) - return - } - parallelFor(n, numWorkers(n), func(lo, hi int) { Log(dst[lo:hi], src[lo:hi]) }) -} diff --git a/internal/kernels/log_amd64.go b/internal/kernels/log_amd64.go deleted file mode 100644 index 44f87ca..0000000 --- a/internal/kernels/log_amd64.go +++ /dev/null @@ -1,54 +0,0 @@ -package kernels - -import "math" - -// logConsts is the constant block logAVX2 reads (asmgen/amd64/gen.go, -// logKernel), each value as a 4-lane vector at the offset the kernel uses. -type logConsts struct { - lo, hiMax, n1loM1, n1hi [4]uint64 // offsets 0, 32, 64, 96 - off, idxMask, topMask [4]uint64 // 128, 160, 192 - magic [4]uint64 // 224 - magicSub [4]float64 // 256 - x2048 [4]uint64 // 288 - negOne, ln2hi, ln2lo [4]float64 // 320, 352, 384 - a0, a1, a2, a3, a4 [4]float64 // 416, 448, 480, 512, 544 -} - -var logK = func() *logConsts { - k := new(logConsts) - for i := range 4 { - k.lo[i], k.hiMax[i] = 0x0010<<48, 0x7fefffffffffffff - k.n1loM1[i] = math.Float64bits(1.0-0x1p-4) - 1 - k.n1hi[i] = math.Float64bits(1.0 + 0x1.09p-4) - k.off[i], k.idxMask[i], k.topMask[i] = logOff, logN-1, 0xfff<<52 - k.magic[i], k.magicSub[i], k.x2048[i] = 0x4330000000000000, 0x1p52+2048, 0x800 - k.negOne[i], k.ln2hi[i], k.ln2lo[i] = -1, logLn2hi, logLn2lo - k.a0[i], k.a1[i], k.a2[i], k.a3[i], k.a4[i] = logA0, logA1, logA2, logA3, logA4 - } - return k -}() - -// useLogAVX2 selects the AVX2 kernel; a var so tests can force the Go path. -var useLogAVX2 = hasAVX2() && useFMA - -// logFast runs blocks of 4 on the AVX2 kernel, which hands back each block -// holding an input Log does not take through its table path (not a positive -// normal, or within 2^-4 of 1) to logScalar, as is the tail. Both run the -// same operations in the same order, so they agree bit for bit. -func logFast(dst, src []float64) { - i := 0 - if useLogAVX2 { - n := len(src) &^ 3 - for i < n { - i += logAVX2(&dst[i], &src[i], n-i, logK, &logTab[0]) - if i < n { - logScalar(dst[i:i+4], src[i:i+4]) - i += 4 - } - } - } - logScalar(dst[i:], src[i:]) -} - -//go:noescape -func logAVX2(dst, src *float64, n int, k *logConsts, tab *logEntry) (done int) diff --git a/internal/kernels/log_amd64_test.go b/internal/kernels/log_amd64_test.go deleted file mode 100644 index f5c955c..0000000 --- a/internal/kernels/log_amd64_test.go +++ /dev/null @@ -1,84 +0,0 @@ -package kernels - -import ( - "math" - "math/rand" - "testing" - "unsafe" -) - -// TestLogConstsLayout holds logConsts to the offsets logKernel uses. -func TestLogConstsLayout(t *testing.T) { - var k logConsts - for _, c := range []struct { - name string - got, want uintptr - }{ - {"lo", unsafe.Offsetof(k.lo), 0}, {"n1hi", unsafe.Offsetof(k.n1hi), 96}, - {"off", unsafe.Offsetof(k.off), 128}, {"topMask", unsafe.Offsetof(k.topMask), 192}, - {"magic", unsafe.Offsetof(k.magic), 224}, {"magicSub", unsafe.Offsetof(k.magicSub), 256}, - {"x2048", unsafe.Offsetof(k.x2048), 288}, {"negOne", unsafe.Offsetof(k.negOne), 320}, - {"ln2lo", unsafe.Offsetof(k.ln2lo), 384}, {"a0", unsafe.Offsetof(k.a0), 416}, - {"a4", unsafe.Offsetof(k.a4), 544}, - } { - if c.got != c.want { - t.Errorf("%s at %d, the kernel reads %d", c.name, c.got, c.want) - } - } - if unsafe.Sizeof(logEntry{}) != 16 { - t.Errorf("logEntry is %d bytes, the kernel gathers at 16i", unsafe.Sizeof(logEntry{})) - } -} - -// TestLogAVX2MatchesGo holds logFast (AVX2 blocks, Go for the blocks it hands -// back and the tail) to logScalar bit for bit (NaN to NaN): random bit -// patterns (most blocks hand back), positive normals across the whole -// exponent range (few do, and every k from -1022 to 1023 and every table -// index occurs), values near 1, the special values, every length 0..70, and -// in place. -func TestLogAVX2MatchesGo(t *testing.T) { - if !useLogAVX2 { - t.Skip("no AVX2+FMA here") - } - r := rand.New(rand.NewSource(9)) - src := make([]float64, 0, 300000) - for _, v := range []float64{1, 1, 1, 1, 0, math.Copysign(0, -1), -1, math.Inf(1), math.Inf(-1), math.NaN(), - 5e-324, 0x1p-1022, math.MaxFloat64, 1 - 0x1p-4, 1 + 0x1.09p-4, math.Nextafter(1-0x1p-4, 0), 0.25, 24.25} { - src = append(src, v) - } - for len(src) < 100000 { - src = append(src, math.Float64frombits(r.Uint64())) - } - for len(src) < 250000 { - src = append(src, math.Float64frombits(0x0010<<48+uint64(r.Int63n(0x7fe0<<48)))) - } - for len(src) < cap(src) { - src = append(src, 0.9+0.2*r.Float64()) - } - want := make([]float64, len(src)) - logScalar(want, src) - same := func(a, b float64) bool { return math.Float64bits(a) == math.Float64bits(b) || (a != a && b != b) } - got := make([]float64, len(src)) - logFast(got, src) - for i := range src { - if !same(got[i], want[i]) { - t.Fatalf("log(0x%016x): AVX2 0x%016x, Go 0x%016x", math.Float64bits(src[i]), math.Float64bits(got[i]), math.Float64bits(want[i])) - } - } - for n := 0; n <= 70; n++ { - d := make([]float64, n) - logFast(d, src[:n]) - for i := range d { - if !same(d[i], want[i]) { - t.Fatalf("length %d [%d]: 0x%016x, want 0x%016x", n, i, math.Float64bits(d[i]), math.Float64bits(want[i])) - } - } - } - in := append([]float64(nil), src[100000:104099]...) - logFast(in, in) - for i := range in { - if !same(in[i], want[100000+i]) { - t.Fatalf("in place [%d]: 0x%016x, want 0x%016x", i, math.Float64bits(in[i]), math.Float64bits(want[100000+i])) - } - } -} diff --git a/internal/kernels/log_fence_amd64_test.go b/internal/kernels/log_fence_amd64_test.go deleted file mode 100644 index fdc3d0c..0000000 --- a/internal/kernels/log_fence_amd64_test.go +++ /dev/null @@ -1,20 +0,0 @@ -//go:build linux || darwin - -package kernels - -import "testing" - -// TestLogAVX2StaysInside fences src and dst of the AVX2 float64 log kernel -// (fences needs mmap/mprotect, hence the build tag of fence_test.go). -func TestLogAVX2StaysInside(t *testing.T) { - if !useLogAVX2 { - t.Skip("no AVX2+FMA here") - } - for n := 4; n <= 64; n += 4 { - for _, s := range fences(t, n) { - for _, d := range fences(t, n) { - logAVX2(&d[0], &s[0], n, logK, &logTab[0]) - } - } - } -} diff --git a/internal/kernels/log_other.go b/internal/kernels/log_other.go deleted file mode 100644 index 6ceb6b6..0000000 --- a/internal/kernels/log_other.go +++ /dev/null @@ -1,6 +0,0 @@ -//go:build !amd64 - -package kernels - -// logFast is Log's Go loop on these architectures (no SIMD kernel yet). -func logFast(dst, src []float64) { logScalar(dst, src) } diff --git a/internal/kernels/log_test.go b/internal/kernels/log_test.go deleted file mode 100644 index 4d601ae..0000000 --- a/internal/kernels/log_test.go +++ /dev/null @@ -1,154 +0,0 @@ -package kernels - -import ( - "math" - "math/big" - "math/rand" - "testing" -) - -const logRefPrec = 300 - -func bigF(x float64) *big.Float { return new(big.Float).SetPrec(logRefPrec).SetFloat64(x) } - -// refAtanhSeries is 2*atanh(t) = 2*sum t^(2n+1)/(2n+1), for |t| <= 1/3. -func refAtanhSeries(t *big.Float) *big.Float { - t2 := new(big.Float).SetPrec(logRefPrec).Mul(t, t) - sum := new(big.Float).SetPrec(logRefPrec) - pow := new(big.Float).SetPrec(logRefPrec).Set(t) - for n := 0; n < 400; n++ { - sum.Add(sum, new(big.Float).SetPrec(logRefPrec).Quo(pow, big.NewFloat(float64(2*n+1)))) - pow.Mul(pow, t2) - } - return sum.Mul(sum, big.NewFloat(2)) -} - -var logRefLn2 = refAtanhSeries(new(big.Float).SetPrec(logRefPrec).Quo(bigF(1), bigF(3))) - -// refLog is log(x) for x > 0 to 300 bits, independent of the code under test: -// x = m * 2^e with m in [1,2), log m = 2*atanh((m-1)/(m+1)). -func refLog(x float64) *big.Float { - m := bigF(x) - e := m.MantExp(m) - 1 // x = m' * 2^(e+1), m' in [0.5,1) -> m = 2m' in [1,2) - m.SetMantExp(m, 1) - t := new(big.Float).SetPrec(logRefPrec).Quo( - new(big.Float).SetPrec(logRefPrec).Sub(m, bigF(1)), - new(big.Float).SetPrec(logRefPrec).Add(m, bigF(1))) - r := refAtanhSeries(t) - return r.Add(r, new(big.Float).SetPrec(logRefPrec).Mul(logRefLn2, big.NewFloat(float64(e)))) -} - -// logULP is |got - exact| in units of the spacing of doubles at exact. -func logULP(got float64, exact *big.Float) float64 { - if exact.Sign() == 0 { - return math.Abs(got) * math.Ldexp(1, 1074) - } - e := exact.MantExp(nil) - 1 - d := bigF(got) - d.Sub(d, exact) - d.SetMantExp(d, -(max(e-52, -1074))) - f, _ := d.Float64() - return math.Abs(f) -} - -// TestLogTable checks every table entry against the properties log_data.c -// says it was built for: c = 1/invc lies in its subinterval of [0x1.6p-1, -// 0x1.6p0), logc is log(c) to within the rounding of invc, and 0x1.8p9 + logc -// is exact (so k*ln2hi + logc is computed without rounding). -func TestLogTable(t *testing.T) { - for i, e := range logTab { - lo := math.Float64frombits(logOff + uint64(i)<<45) - hi := math.Float64frombits(logOff + uint64(i+1)<<45) - if c := 1 / e.invc; c < lo || c >= hi { - t.Errorf("entry %d: 1/invc = %v outside [%v, %v)", i, c, lo, hi) - } - neg := refLog(e.invc) - d := bigF(e.logc) - d.Add(d, neg) // logc - log(c) ~= logc + log(invc) - if f, _ := d.Float64(); math.Abs(f) > 0x1p-52 { - t.Errorf("entry %d: logc %v differs from -log(invc) by %v", i, e.logc, f) - } - if (0x1.8p9+e.logc)-0x1.8p9 != e.logc { - t.Errorf("entry %d: 0x1.8p9 + logc is not exact", i) - } - } -} - -// TestLogULP measures the worst error over random inputs across the whole -// positive range (log-uniform exponent), the region near 1 and subnormals. -func TestLogULP(t *testing.T) { - rng := rand.New(rand.NewSource(3)) - var xs []float64 - for i := 0; i < 2500; i++ { - xs = append(xs, - math.Ldexp(1+rng.Float64(), rng.Intn(2046)-1022), // whole normal range - 1+(rng.Float64()*2-1)*0x1.2p-4, // near 1 (both paths) - math.Ldexp(rng.Float64(), -1022), // subnormal - 0.5+rng.Float64()*2, // the table's own range - ) - } - worst, worstX, goWorst := 0.0, 0.0, 0.0 - for _, x := range xs { - if x <= 0 { - continue - } - ref := refLog(x) - if e := logULP(log1(x), ref); e > worst { - worst, worstX = e, x - } - goWorst = max(goWorst, logULP(math.Log(x), ref)) - } - t.Logf("log: worst %.4f ULP at x=%v over %d inputs (math.Log: %.4f)", worst, worstX, len(xs), goWorst) - if worst > 0.53 { - t.Errorf("log error %.4f ULP at x=%v exceeds 0.53", worst, worstX) - } -} - -func TestLogSpecial(t *testing.T) { - for _, c := range []struct{ x, want float64 }{ - {1, 0}, {0, math.Inf(-1)}, {math.Copysign(0, -1), math.Inf(-1)}, - {math.Inf(1), math.Inf(1)}, {math.E, 1}, {2, math.Ln2}, - } { - if got := log1(c.x); got != c.want { - t.Errorf("log(%v) = %v, want %v", c.x, got, c.want) - } - } - for _, x := range []float64{-1, math.Inf(-1), math.NaN(), -1e-310} { - if got := log1(x); !math.IsNaN(got) { - t.Errorf("log(%v) = %v, want NaN", x, got) - } - } - if got := log1(5e-324); got != -744.4400719213812 { - t.Errorf("log(5e-324) = %v", got) - } -} - -// TestLogLoopMatchesLog1: Log's loop and log1 give the same bits everywhere. -func TestLogLoopMatchesLog1(t *testing.T) { - rng := rand.New(rand.NewSource(4)) - xs := []float64{0, negZero, 1, 1 - 0x1p-4, 1 + 0x1.09p-4, 0x1p-1022, 5e-324, -3, - math.Inf(1), math.Inf(-1), math.NaN(), math.MaxFloat64} - for i := 0; i < 20000; i++ { - xs = append(xs, math.Ldexp(rng.Float64()*2-0.5, rng.Intn(2100)-1074), 1+(rng.Float64()-0.5)*0.3) - } - got := make([]float64, len(xs)) - Log(got, xs) - for i, x := range xs { - w := log1(x) - if math.Float64bits(got[i]) != math.Float64bits(w) && !(math.IsNaN(w) && math.IsNaN(got[i])) { - t.Fatalf("Log(%v) = %v, log1 = %v", x, got[i], w) - } - } - src := randVec(3000, 9) - for i := range src { - src[i] = math.Abs(src[i]) * 100 - } - want, par := make([]float64, len(src)), make([]float64, len(src)) - Log(want, src) - withThresholds(64, 1<<14, func() { LogP(par, src) }) - eqSlice(t, par, want) - LogP(par[:5], src[:5]) - Log10(want, src) - withThresholds(64, 1<<14, func() { Log10P(par, src) }) - eqSlice(t, par, want) -} diff --git a/internal/kernels/logcr.go b/internal/kernels/logcr.go deleted file mode 100644 index 01e00f8..0000000 --- a/internal/kernels/logcr.go +++ /dev/null @@ -1,138 +0,0 @@ -package kernels - -import ( - "math" - "math/bits" -) - -// logCR is CORE-MATH's correctly rounded binary32 logarithm (Alexei Sibidanov -// and Paul Zimmermann, src/binary32/log/logf.c, MIT), one C operation per Go -// operation. x = 2^e·m, m in [1,2); tr[j] ≈ 1/m on 64 intervals so that -// z = m·tr[j] - 1 is small; log x = e·ln2 + tl[j] + log1p(z), with -// tl[j] = -log(tr[j]). A degree-3 polynomial gives a lower bound r of the -// exact value with r + 0x1.f06p-33 an upper bound; when both round to the -// same float32 that float32 is the answer. Otherwise (762 750 of the 2^31 -// positive finite inputs, about 1 in 2800, measured by TestLogCRMatchesLog32) -// a degree-8 polynomial with a split ln2 decides it. -// -// Every product is converted explicitly to float64, which Go defines as -// rounding it, so no architecture fuses it into an FMA: the operations are -// the C source's, unfused, everywhere. TestLogCRMatchesLog32 holds the result -// to log32, already proven correctly rounded, over every float32. -func logCR(x float32) float32 { - ux := math.Float32bits(x) - if ux < 1<<23 || ux >= 0x7f800000 { - if ux == 0 || ux >= 0x7f800000 { - return logCRSpecial(x) - } - n := bits.LeadingZeros32(ux) - 8 // subnormal: normalise the mantissa - ux <<= n - ux -= uint32(n) << 23 - } - if ux == 127<<23 { - return 0 - } - m := ux & (1<<23 - 1) - j := (m + 1<<(23-7)) >> (23 - 6) - e := float64(int32(ux)>>23 - 127) - z := float64(math.Float64frombits((uint64(m)|1023<<23)<<(52-23))*logCRTr[j]) - 1 - z2 := float64(z * z) - r := ((float64(e*0x1.62e42fefa39efp-1) + logCRTl[j]) + float64(z*logCRB[0])) + - float64(z2*(logCRB[1]+float64(z*logCRB[2]))) - ub, lb := float32(r), float32(r+0x1.f06p-33) - if ub == lb { - return ub - } - return logCRSlow(x, z, z2, e, j) -} - -// logCRSlow is logf.c's second step, for the inputs whose bounds straddle a -// float32 rounding boundary. -func logCRSlow(x float32, z, z2, e float64, j uint32) float32 { - c := &logCRC - f := float64(z2 * ((c[0] + float64(z*c[1])) + float64(z2*((c[2]+float64(z*c[3]))+ - float64(z2*(c[4]+float64(z*c[5])+float64(z2*c[6]))))))) - if math.Abs(float64(x-1)) < 0x1p-10 { - return float32(z + f) - } - f -= float64(0x1.0ca86c3898dp-49 * e) - f += z - f += logCRTl[j] - logCRTl[0] - el := float64(e * 0x1.62e42fefa3ap-1) - r := el + f - ub := float32(r) - if math.Float64bits(r)&(1<<28-1) == 0 { - dr := (el - r) + f - r += float64(dr * 64) - ub = float32(r) - } - return ub -} - -// logCRSpecial is logf.c's as_special: ±0 → -Inf, +Inf → +Inf, NaN → NaN, -// negatives → NaN. -func logCRSpecial(x float32) float32 { - switch ux := math.Float32bits(x); { - case ux<<1 == 0: - return float32(math.Inf(-1)) - case ux == 0x7f800000: - return x - case ux<<1 > 0xff000000: - return x + x - } - return float32(math.NaN()) -} - -// logCRSlice writes logCR(src[i]) into dst[i]. -func logCRSlice(dst, src []float32) { - src = src[:len(dst)] - for i, x := range src { - dst[i] = logCR(x) - } -} - -var logCRTr = [65]float64{ - 0x1p+0, 0x1.f81f82p-1, 0x1.f07c1fp-1, 0x1.e9131acp-1, - 0x1.e1e1e1ep-1, 0x1.dae6077p-1, 0x1.d41d41dp-1, 0x1.cd85689p-1, - 0x1.c71c71cp-1, 0x1.c0e0704p-1, 0x1.bacf915p-1, 0x1.b4e81b5p-1, - 0x1.af286bdp-1, 0x1.a98ef6p-1, 0x1.a41a41ap-1, 0x1.9ec8e95p-1, - 0x1.999999ap-1, 0x1.948b0fdp-1, 0x1.8f9c19p-1, 0x1.8acb90fp-1, - 0x1.8618618p-1, 0x1.8181818p-1, 0x1.7d05f41p-1, 0x1.78a4c81p-1, - 0x1.745d174p-1, 0x1.702e05cp-1, 0x1.6c16c17p-1, 0x1.6816817p-1, - 0x1.642c859p-1, 0x1.605816p-1, 0x1.5c9882cp-1, 0x1.58ed231p-1, - 0x1.5555555p-1, 0x1.51d07ebp-1, 0x1.4e5e0a7p-1, 0x1.4afd6ap-1, - 0x1.47ae148p-1, 0x1.446f865p-1, 0x1.4141414p-1, 0x1.3e22cbdp-1, - 0x1.3b13b14p-1, 0x1.3813814p-1, 0x1.3521cfbp-1, 0x1.323e34ap-1, - 0x1.2f684bep-1, 0x1.2c9fb4ep-1, 0x1.29e412ap-1, 0x1.27350b9p-1, - 0x1.2492492p-1, 0x1.21fb781p-1, 0x1.1f7047ep-1, 0x1.1cf06aep-1, - 0x1.1a7b961p-1, 0x1.1811812p-1, 0x1.15b1e5fp-1, 0x1.135c811p-1, - 0x1.1111111p-1, 0x1.0ecf56cp-1, 0x1.0c9715p-1, 0x1.0a6810ap-1, - 0x1.0842108p-1, 0x1.0624dd3p-1, 0x1.041041p-1, 0x1.0204081p-1, 0.5, -} - -var logCRTl = [65]float64{ - -0x1.3b40815cd0628p-45, 0x1.fc0a890fbb514p-7, 0x1.f829b1e780b98p-6, 0x1.77458f532c948p-5, - 0x1.f0a30c2114ef2p-5, 0x1.341d793bbc7f7p-4, 0x1.6f0d28d256172p-4, 0x1.a926d3a6acb89p-4, - 0x1.e2707722ae90cp-4, 0x1.0d77e7a90896cp-3, 0x1.29552f6fff036p-3, 0x1.44d2b6c5b7831p-3, - 0x1.5ff306ee78ee7p-3, 0x1.7ab890410d41cp-3, 0x1.9525a9e3451c7p-3, 0x1.af3c94ed0bb06p-3, - 0x1.c8ff7c59a9535p-3, 0x1.e27076d5aedf9p-3, 0x1.fb9186b5e393ep-3, 0x1.0a324e38b8e6dp-2, - 0x1.1675cacaba398p-2, 0x1.22941fc0f76efp-2, 0x1.2e8e2bc311abap-2, 0x1.3a64c56b14373p-2, - 0x1.4618bc31c5c4cp-2, 0x1.51aad874df5b7p-2, 0x1.5d1bdbea80754p-2, 0x1.686c81d331238p-2, - 0x1.739d7f6dbcd9p-2, 0x1.7eaf83c82ad4dp-2, 0x1.89a3385813fe4p-2, 0x1.947941aa91484p-2, - 0x1.9f323edbf95d5p-2, 0x1.a9cec9a4205d3p-2, 0x1.b44f77c5c8cecp-2, 0x1.beb4d9ea71905p-2, - 0x1.c8ff7c69a97abp-2, 0x1.d32fe7f38e95fp-2, 0x1.dd46a0501c22ap-2, 0x1.e7442617e8511p-2, - 0x1.f128f5eaf0476p-2, 0x1.faf588dd8f0a8p-2, 0x1.02552a5edcfc4p-1, 0x1.0723e5c64de05p-1, - 0x1.0be72e3852947p-1, 0x1.109f39d554b5cp-1, 0x1.154c3d2c4d4aep-1, 0x1.19ee6b38bc834p-1, - 0x1.1e85f5ef03f95p-1, 0x1.23130d7fabe07p-1, 0x1.2795e1219afep-1, 0x1.2c0e9ec9c8d5p-1, - 0x1.307d7337f0f83p-1, 0x1.34e289cb4e098p-1, 0x1.393e0d42e28dep-1, 0x1.3d9026ad555bfp-1, - 0x1.41d8fe8667173p-1, 0x1.4618bc1ec5d87p-1, 0x1.4a4f85d303d8p-1, 0x1.4e7d8127f5a75p-1, - 0x1.52a2d26dbc47p-1, 0x1.56bf9d597f25ep-1, 0x1.5ad404cb59df2p-1, 0x1.5ee02a928153ap-1, - 0x1.62e42fefa38b4p-1, -} - -var logCRB = [3]float64{0x1.00000006342eap+0, -0x1.0001f7fdc3977p-1, 0x1.554a4e5cae9cfp-2} - -var logCRC = [7]float64{ - -0x1p-1, 0x1.55555555571cap-2, -0x1.0000000002d85p-2, 0x1.9999987d0c963p-3, - -0x1.555554059a8bbp-3, 0x1.24aebcf71a38fp-3, -0x1.001c73915d758p-3, -} diff --git a/internal/kernels/logcr_amd64.go b/internal/kernels/logcr_amd64.go deleted file mode 100644 index 340032e..0000000 --- a/internal/kernels/logcr_amd64.go +++ /dev/null @@ -1,53 +0,0 @@ -package kernels - -// logCRConsts is the constant block log32CRAVX2 reads (asmgen/amd64/gen.go, -// log32CRKernel), at the offsets the kernel uses. -type logCRConsts struct { - one, ln2, b0, b1, b2, w [4]float64 // offsets 0, 32, 64, 96, 128, 160 - sub, maxv, mmask, jadd, bias, oneBits [8]uint32 // 192, 224, 256, 288, 320, 352 - expOr [4]uint64 // 384 - tr, tl [65]float64 // 416, 936 -} - -var logCRK = func() *logCRConsts { - k := new(logCRConsts) - for i := range 4 { - k.one[i], k.ln2[i], k.w[i] = 1, 0x1.62e42fefa39efp-1, 0x1.f06p-33 - k.b0[i], k.b1[i], k.b2[i] = logCRB[0], logCRB[1], logCRB[2] - k.expOr[i] = 1023 << 23 - } - for i := range 8 { - // The kernel's positive-normal test: ux + sub (wrapping) <= maxv as int32. - k.sub[i], k.maxv[i], k.mmask[i] = 0x7f800000, 0xfeffffff, 1<<23-1 - k.jadd[i], k.bias[i], k.oneBits[i] = 1<<16, 127, 0x3f800000 - } - k.tr, k.tl = logCRTr, logCRTl - return k -}() - -// useLog32AVX2 selects the AVX2 kernel; a var so tests can force the Go path. -var useLog32AVX2 = hasAVX2() - -// log32Fast is the correctly rounded float32 log: blocks of 8 on the AVX2 -// kernel, which hands back each block it cannot finish (an input that is not -// a positive normal float32, or bounds that round apart, about 1 lane in -// 2800) to logCR, as is the tail. Both run logCR's operations in its order, -// so they agree bit for bit. -func log32Fast(dst, src []float32) { - src = src[:len(dst)] - i := 0 - if useLog32AVX2 { - n := len(dst) &^ 7 - for i < n { - i += log32CRAVX2(&dst[i], &src[i], n-i, logCRK) - if i < n { - logCRSlice(dst[i:i+8], src[i:i+8]) - i += 8 - } - } - } - logCRSlice(dst[i:], src[i:]) -} - -//go:noescape -func log32CRAVX2(dst, src *float32, n int, k *logCRConsts) (done int) diff --git a/internal/kernels/logcr_amd64_test.go b/internal/kernels/logcr_amd64_test.go deleted file mode 100644 index 624f445..0000000 --- a/internal/kernels/logcr_amd64_test.go +++ /dev/null @@ -1,91 +0,0 @@ -package kernels - -import ( - "math" - "math/rand" - "testing" - "unsafe" -) - -// TestLogCRConstsLayout holds logCRConsts to the offsets log32CRKernel uses. -func TestLogCRConstsLayout(t *testing.T) { - var k logCRConsts - for _, c := range []struct { - name string - got, want uintptr - }{ - {"one", unsafe.Offsetof(k.one), 0}, {"ln2", unsafe.Offsetof(k.ln2), 32}, - {"b0", unsafe.Offsetof(k.b0), 64}, {"b1", unsafe.Offsetof(k.b1), 96}, - {"b2", unsafe.Offsetof(k.b2), 128}, {"w", unsafe.Offsetof(k.w), 160}, - {"sub", unsafe.Offsetof(k.sub), 192}, {"maxv", unsafe.Offsetof(k.maxv), 224}, - {"mmask", unsafe.Offsetof(k.mmask), 256}, {"jadd", unsafe.Offsetof(k.jadd), 288}, - {"bias", unsafe.Offsetof(k.bias), 320}, {"oneBits", unsafe.Offsetof(k.oneBits), 352}, - {"expOr", unsafe.Offsetof(k.expOr), 384}, {"tr", unsafe.Offsetof(k.tr), 416}, - {"tl", unsafe.Offsetof(k.tl), 936}, - } { - if c.got != c.want { - t.Errorf("%s at %d, the kernel reads %d", c.name, c.got, c.want) - } - } -} - -// TestLog32AVX2MatchesGo holds log32Fast (AVX2 blocks, Go for the blocks it -// hands back and the tail) to logCR bit for bit (NaN to NaN): random bit -// patterns (a third of them negative, NaN or subnormal, so most blocks bail), -// random positive normals (few bail), 1.0 runs, the special values, and every -// length 0..70. -func TestLog32AVX2MatchesGo(t *testing.T) { - if !useLog32AVX2 { - t.Skip("no AVX2 here") - } - r := rand.New(rand.NewSource(5)) - src := make([]float32, 0, 300000) - for _, v := range []float32{1, 1, 1, 1, 1, 1, 1, 1, 0, float32(math.Copysign(0, -1)), -1, - float32(math.Inf(1)), float32(math.Inf(-1)), float32(math.NaN()), math.SmallestNonzeroFloat32, - 0x1p-126, math.MaxFloat32, 0.99999994, 1.0000001, 2, 0.5} { - src = append(src, v) - } - for len(src) < 100000 { - src = append(src, math.Float32frombits(r.Uint32())) - } - for len(src) < cap(src) { - src = append(src, math.Float32frombits(0x00800000+uint32(r.Int63n(0x7f000000)))) - } - want := make([]float32, len(src)) - logCRSlice(want, src) - got := make([]float32, len(src)) - log32Fast(got, src) - same := func(a, b float32) bool { return math.Float32bits(a) == math.Float32bits(b) || (a != a && b != b) } - for i := range src { - if !same(got[i], want[i]) { - t.Fatalf("log(0x%08x): AVX2 0x%08x, Go 0x%08x", math.Float32bits(src[i]), math.Float32bits(got[i]), math.Float32bits(want[i])) - } - } - for n := 0; n <= 70; n++ { - d := make([]float32, n) - log32Fast(d, src[:n]) - for i := range d { - if !same(d[i], want[i]) { - t.Fatalf("length %d [%d]: 0x%08x, want 0x%08x", n, i, math.Float32bits(d[i]), math.Float32bits(want[i])) - } - } - } -} - -// TestLog32AVX2InPlace runs the kernel with dst == src: a block it hands -// back must still hold its inputs. -func TestLog32AVX2InPlace(t *testing.T) { - r := rand.New(rand.NewSource(6)) - src := make([]float32, 4099) - for i := range src { - src[i] = math.Float32frombits(r.Uint32() &^ (1 << 31)) // NaN, Inf and subnormals among them - } - want := make([]float32, len(src)) - logCRSlice(want, src) - log32Fast(src, src) - for i := range src { - if math.Float32bits(src[i]) != math.Float32bits(want[i]) && !(src[i] != src[i] && want[i] != want[i]) { - t.Fatalf("[%d]: 0x%08x, want 0x%08x", i, math.Float32bits(src[i]), math.Float32bits(want[i])) - } - } -} diff --git a/internal/kernels/logcr_fence_amd64_test.go b/internal/kernels/logcr_fence_amd64_test.go deleted file mode 100644 index 26c1604..0000000 --- a/internal/kernels/logcr_fence_amd64_test.go +++ /dev/null @@ -1,20 +0,0 @@ -//go:build linux || darwin - -package kernels - -import "testing" - -// TestLog32AVX2StaysInside fences src and dst of the AVX2 log kernel -// (fences32 needs mmap/mprotect, hence the build tag of fence_test.go). -func TestLog32AVX2StaysInside(t *testing.T) { - if !useLog32AVX2 { - t.Skip("no AVX2 here") - } - for n := 8; n <= 64; n += 8 { - for _, s := range fences32(t, n) { - for _, d := range fences32(t, n) { - log32CRAVX2(&d[0], &s[0], n, logCRK) - } - } - } -} diff --git a/internal/kernels/logcr_test.go b/internal/kernels/logcr_test.go deleted file mode 100644 index f28152d..0000000 --- a/internal/kernels/logcr_test.go +++ /dev/null @@ -1,150 +0,0 @@ -package kernels - -import ( - "fmt" - "math" - "math/bits" - "runtime" - "sync" - "sync/atomic" - "testing" -) - -// TestLogCRMatchesLog32 holds logCR to log32 bit for bit (NaN to NaN) on -// every float32 input (NDARRAY_EXHAUSTIVE=1; otherwise every 65537th). -// log32 is proven correctly rounded over all 2^32 inputs by -// TestLog32CorrectlyRounded, so equality everywhere proves logCR is too. It -// also counts the inputs logCRSlow decides. -func TestLogCRMatchesLog32(t *testing.T) { - stride := exhaustiveStride() - const batch = 4096 - var next atomic.Uint64 - var bad, slow, checked atomic.Int64 - var mu sync.Mutex - var first []string - var wg sync.WaitGroup - for range runtime.GOMAXPROCS(0) { - wg.Add(1) - go func() { - defer wg.Done() - xs := make([]float32, 0, batch) - want := make([]float32, batch) - for { - start := next.Add(batch*stride) - batch*stride - if start >= 1<<32 { - return - } - xs = xs[:0] - for b := start; b < start+batch*stride && b < 1<<32; b += stride { - xs = append(xs, math.Float32frombits(uint32(b))) - } - log32(want[:len(xs)], xs) - for i, x := range xs { - checked.Add(1) - got := logCR(x) - if logCRTakesSlow(x) { - slow.Add(1) - } - w := want[i] - if math.Float32bits(got) != math.Float32bits(w) && !(got != got && w != w) { - bad.Add(1) - mu.Lock() - if len(first) < 8 { - first = append(first, "logCR("+formatF32(x)+") = "+formatF32(got)+", log32 "+formatF32(w)) - } - mu.Unlock() - } - } - } - }() - } - wg.Wait() - t.Logf("%d inputs checked, %d through the slow path, %d differ from log32", checked.Load(), slow.Load(), bad.Load()) - for _, s := range first { - t.Error(s) - } -} - -// logCRTakesSlow repeats logCR's fast path and reports whether its bounds -// straddle a rounding boundary. -func logCRTakesSlow(x float32) bool { - ux := math.Float32bits(x) - if ux == 0 || ux >= 0x7f800000 || ux == 127<<23 { - return false - } - if ux < 1<<23 { - n := bits.LeadingZeros32(ux) - 8 - ux <<= n - ux -= uint32(n) << 23 - } - m := ux & (1<<23 - 1) - j := (m + 1<<(23-7)) >> (23 - 6) - e := float64(int32(ux)>>23 - 127) - z := float64(math.Float64frombits((uint64(m)|1023<<23)<<(52-23))*logCRTr[j]) - 1 - z2 := float64(z * z) - r := ((float64(e*0x1.62e42fefa39efp-1) + logCRTl[j]) + float64(z*logCRB[0])) + - float64(z2*(logCRB[1]+float64(z*logCRB[2]))) - return float32(r) != float32(r+0x1.f06p-33) -} - -// BenchmarkLog32Kernels compares the float32 log kernels on two inputs in -// [0.25, 24.25] at 16 Ki and 1 Mi: a ramp, and the multiples of 0.25 that -// benchmarks/BenchmarkDTypeExpLog uses (0.25 + (i%97)*0.25). It reports the -// blocks of 8 the AVX2 kernel hands back to logCR per call. -func BenchmarkLog32Kernels(b *testing.B) { - for _, data := range []struct { - name string - at func(i, n int) float32 - }{ - {"ramp", func(i, n int) float32 { return 0.25 + 24*float32(i)/float32(n) }}, - {"steps", func(i, n int) float32 { return float32(0.25 + float64(i%97)*0.25) }}, - } { - for _, n := range []int{16 << 10, 1 << 20} { - src := make([]float32, n) - for i := range src { - src[i] = data.at(i, n) - } - handed := 0 - for i := 0; i+8 <= n; i += 8 { - for _, x := range src[i : i+8] { - if logCRTakesSlow(x) { - handed++ - break - } - } - } - dst := make([]float32, n) - for _, k := range []struct { - name string - f func(dst, src []float32) - }{{"via64", log32}, {"logCR", logCRSlice}, {"fast", log32Fast}} { - b.Run(fmt.Sprintf("%s/n=%d/%s", data.name, n, k.name), func(b *testing.B) { - b.SetBytes(int64(4 * n)) - for b.Loop() { - k.f(dst, src) - } - b.ReportMetric(float64(handed), "handed-blocks/op") - }) - } - } - } -} - -func TestLogCRSpecial(t *testing.T) { - for _, c := range []struct{ x, want float32 }{ - {0, float32(math.Inf(-1))}, - {float32(math.Copysign(0, -1)), float32(math.Inf(-1))}, - {float32(math.Inf(1)), float32(math.Inf(1))}, - {1, 0}, - {math.SmallestNonzeroFloat32, -103.27893}, - } { - if got := logCR(c.x); got != c.want { - t.Errorf("logCR(%v) = %v, want %v", c.x, got, c.want) - } - } - for _, x := range []float32{-1, float32(math.Inf(-1)), float32(math.NaN()), -math.SmallestNonzeroFloat32} { - if got := logCR(x); got == got { - t.Errorf("logCR(%v) = %v, want NaN", x, got) - } - } -} diff --git a/internal/kernels/parallel.go b/internal/kernels/parallel.go index 29de022..85f03f2 100644 --- a/internal/kernels/parallel.go +++ b/internal/kernels/parallel.go @@ -266,12 +266,14 @@ func mapReduceP(a []float64, red func([]float64) float64, combine func([]float64 return combine(partials) } -// SumP returns the sum of all elements of a, computed as the SIMD-summed -// per-chunk partials folded by a final Sum. Floating-point addition is not -// associative, so this grouping can differ from the strictly-sequential Sum by -// a few ULP — the same trade-off numpy's pairwise summation makes; callers that -// need the exact sequential value use Sum. -func SumP(a []float64) float64 { return mapReduceP(a, sumSIMD, Sum) } +// SumP returns the sum of all elements of a in NumPy's order: floats.SumPairwise, +// bit-identical to np.sum of a contiguous 1-D float64 array. It is serial at every +// size: NumPy's pairwise order is a fixed recursion, and summing chunks on several +// cores and folding the partials would regroup it (the result would then depend on +// the core count). An order-preserving parallel split (chunks cut where NumPy's +// recursion cuts) is possible and not done here. Callers that need the exact +// sequential left-to-right value use Sum. +func SumP(a []float64) float64 { return sumSIMD(a) } // reduceP folds per-chunk results of the associative scalar reducer red with // red itself, giving a result identical to the serial kernel. diff --git a/internal/kernels/reduce_amd64.go b/internal/kernels/reduce_amd64.go deleted file mode 100644 index a378763..0000000 --- a/internal/kernels/reduce_amd64.go +++ /dev/null @@ -1,124 +0,0 @@ -package kernels - -// amd64 SSE2 kernels generated by go-asmgen (sum_amd64.s). SSE2 is part of the -// amd64 baseline (GOAMD64=v1), so every kernel here is always callable — there -// is no CPU-feature branch. -// -//go:generate sh -c "cd asmgen/amd64 && GOFLAGS=-mod=mod go run gen.go && mv sum_amd64.s ../../sum_amd64.s" - -// HaveReduceSIMD reports that this build ships SIMD reduction/elementwise -// kernels (so the dispatch wrappers route through them rather than the scalar -// oracle). True on amd64/arm64; false on the four arches without vector double. -const HaveReduceSIMD = true - -// sumSIMD returns the SIMD sum of a (empty -> 0), using four SSE2 accumulators. -// The lane-parallel grouping can differ from sequential Sum by a few ULP. -func sumSIMD(a []float64) float64 { - if len(a) == 0 { - return 0 - } - return sumSSE2(&a[0], len(a)) -} - -// sqrtSIMD writes sqrt(src[i]) into dst[i] with the packed SSE2 SQRTPD kernel. -// SQRTPD computes the same correctly-rounded IEEE square root the scalar oracle -// (math.Sqrt) does, lane-by-lane, so it is bit-identical to sqrtScalar including -// NaN/Inf/signed-zero. len(dst) == len(src) is guaranteed by the caller. -func sqrtSIMD(dst, src []float64) { - if len(dst) == 0 { - return - } - _ = src[len(dst)-1] // the kernel reads len(dst) of src - sqrtSSE2(&dst[0], &src[0], len(dst)) -} - -// maxSIMD returns the NaN-propagating maximum of a (non-empty): if any element -// is NaN the result is NaN, otherwise the largest element — matching the -// math.Max-based scalar oracle (Max) and numpy.max bit-for-bit. The SSE2 kernel -// accumulates a packed running max with MAXPD and a separate "saw NaN" mask via -// CMPUNORDPD, then returns NaN if any lane ever held a NaN. -// -// MAXPD returns its source operand when both inputs are zeros, so the kernel's -// zero is +0 or -0 depending on where the zeros sit. A zero result is the one -// case where that matters, and it needs every element <= 0 with a zero among -// them: the scalar oracle then gives the sign Go's max uses (+0 if any +0). -// Doing the same in the kernel, both operand orders combined at every step, -// cost 53-60% on a Zen 3; this costs nothing on any other result. -func maxSIMD(a []float64) float64 { - if r := maxSSE2(&a[0], len(a)); r != 0 { - return r - } - return Max(a) -} - -// minSIMD returns the NaN-propagating minimum of a (non-empty), the dual of -// maxSIMD via MINPD; bit-identical to the scalar Min oracle and numpy.min. A -// zero result is resolved by the oracle, as in maxSIMD (-0 if any -0). -func minSIMD(a []float64) float64 { - if r := minSSE2(&a[0], len(a)); r != 0 { - return r - } - return Min(a) -} - -// addBin/subBin/mulBin/divBin write a[i] OP b[i] into dst[i] with the packed -// SSE2 kernels (ADDPD/SUBPD/MULPD/DIVPD, 8 doubles/iter + scalar tail). Each -// lane is an exact, independent IEEE-754 op, so the result is bit-identical to -// the scalar Add/Sub/Mul/Div oracle (no reduction grouping). They are the -// serial elementwise inner loop the per-op fast path and runBinaryP use; len -// agreement is guaranteed by the caller. Empty slices are a no-op (no &a[0]). -func addBin(dst, a, b []float64) { - if len(dst) == 0 { - return - } - _, _ = a[len(dst)-1], b[len(dst)-1] // the kernel reads len(dst) of each - addSSE2(&dst[0], &a[0], &b[0], len(dst)) -} - -func subBin(dst, a, b []float64) { - if len(dst) == 0 { - return - } - _, _ = a[len(dst)-1], b[len(dst)-1] // the kernel reads len(dst) of each - subSSE2(&dst[0], &a[0], &b[0], len(dst)) -} - -func mulBin(dst, a, b []float64) { - if len(dst) == 0 { - return - } - _, _ = a[len(dst)-1], b[len(dst)-1] // the kernel reads len(dst) of each - mulSSE2(&dst[0], &a[0], &b[0], len(dst)) -} - -func divBin(dst, a, b []float64) { - if len(dst) == 0 { - return - } - _, _ = a[len(dst)-1], b[len(dst)-1] // the kernel reads len(dst) of each - divSSE2(&dst[0], &a[0], &b[0], len(dst)) -} - -//go:noescape -func sumSSE2(a *float64, n int) float64 - -//go:noescape -func sqrtSSE2(dst, src *float64, n int) - -//go:noescape -func maxSSE2(a *float64, n int) float64 - -//go:noescape -func minSSE2(a *float64, n int) float64 - -//go:noescape -func addSSE2(dst, a, b *float64, n int) - -//go:noescape -func subSSE2(dst, a, b *float64, n int) - -//go:noescape -func mulSSE2(dst, a, b *float64, n int) - -//go:noescape -func divSSE2(dst, a, b *float64, n int) diff --git a/internal/kernels/reduce_arm64.go b/internal/kernels/reduce_arm64.go deleted file mode 100644 index c3654fd..0000000 --- a/internal/kernels/reduce_arm64.go +++ /dev/null @@ -1,118 +0,0 @@ -package kernels - -// arm64 kernels. The sum reduction is a hand-vectorized NEON kernel generated by -// go-asmgen (sum_arm64.s, using VFMLA against a vector of 1.0 — see asmgen). -// -//go:generate sh -c "cd asmgen/arm64 && GOFLAGS=-mod=mod go run gen.go && mv sum_arm64.s ../../sum_arm64.s" -// -// Every kernel here is NEON: sum, dot, sqrt, add/sub/mul/div, max/min and the -// GEMM micro-kernel. The vector sqrt, divide and max/min use mnemonics Go's -// assembler gained in 1.27 (VFSQRT, VFDIV, VFMAX, VFMIN); before that, sqrt was -// a WORD and divide and max/min stayed scalar. - -// HaveReduceSIMD reports that this build ships at least one hand-vectorized SIMD -// kernel (the NEON sum). The per-arch CI execution job validates it against the -// scalar oracle to a tight relative tolerance. -const HaveReduceSIMD = true - -// sumSIMD returns the SIMD sum of a (empty -> 0), using four NEON accumulators. -// The lane-parallel grouping can differ from sequential Sum by a few ULP. -func sumSIMD(a []float64) float64 { - if len(a) == 0 { - return 0 - } - return sumNEON(&a[0], len(a)) -} - -// sqrtSIMD writes sqrt(src[i]) into dst[i] with the packed NEON FSQRT V.2D kernel -// (two doubles per instruction, VFSQRT in Go). FSQRT V.2D is correctly- -// rounded IEEE-754, identical lane-by-lane to the scalar FSQRTD/math.Sqrt, so it -// is bit-identical to sqrtScalar including negatives->NaN, ±Inf and signed zeros. -// It supersedes the prior scalar-FSQRTD loop, halving the per-element sqrt count. -func sqrtSIMD(dst, src []float64) { - if len(dst) == 0 { - return - } - _ = src[len(dst)-1] // the kernel reads len(dst) of src - sqrtNEON(&dst[0], &src[0], len(dst)) -} - -// maxSIMD / minSIMD are the NaN-propagating extrema of a (non-empty): the NEON -// FMAX/FMIN kernel with four D2 accumulators from eight elements up, the -// scalar oracle below that. Bit-identical to the serial Max/Min -// oracle (see extremumKernel in asmgen/arm64/gen.go). -func maxSIMD(a []float64) float64 { - if len(a) < 8 { - return Max(a) - } - return maxNEON(&a[0], len(a)) -} - -func minSIMD(a []float64) float64 { - if len(a) < 8 { - return Min(a) - } - return minNEON(&a[0], len(a)) -} - -// addBin/subBin/mulBin/divBin write a[i] OP b[i] into dst[i] with the NEON D2 -// kernels (VFADD/VFSUB/VFMUL/VFDIV, 8 doubles/iter + scalar tail). Each lane is -// one correctly rounded IEEE-754 operation, so the result is bit-identical to -// the scalar Add/Sub/Mul/Div oracle. They are the serial -// elementwise inner loop the per-op fast path and runBinaryP use. Empty slices -// are a no-op (no &a[0]). -func addBin(dst, a, b []float64) { - if len(dst) == 0 { - return - } - _, _ = a[len(dst)-1], b[len(dst)-1] // the kernel reads len(dst) of each - addNEON(&dst[0], &a[0], &b[0], len(dst)) -} - -func subBin(dst, a, b []float64) { - if len(dst) == 0 { - return - } - _, _ = a[len(dst)-1], b[len(dst)-1] // the kernel reads len(dst) of each - subNEON(&dst[0], &a[0], &b[0], len(dst)) -} - -func mulBin(dst, a, b []float64) { - if len(dst) == 0 { - return - } - _, _ = a[len(dst)-1], b[len(dst)-1] // the kernel reads len(dst) of each - mulNEON(&dst[0], &a[0], &b[0], len(dst)) -} - -func divBin(dst, a, b []float64) { - if len(dst) == 0 { - return - } - _, _ = a[len(dst)-1], b[len(dst)-1] // the kernel reads len(dst) of each - divNEON(&dst[0], &a[0], &b[0], len(dst)) -} - -//go:noescape -func sumNEON(a *float64, n int) float64 - -//go:noescape -func sqrtNEON(dst, src *float64, n int) - -//go:noescape -func addNEON(dst, a, b *float64, n int) - -//go:noescape -func subNEON(dst, a, b *float64, n int) - -//go:noescape -func mulNEON(dst, a, b *float64, n int) - -//go:noescape -func divNEON(dst, a, b *float64, n int) - -//go:noescape -func maxNEON(a *float64, n int) float64 - -//go:noescape -func minNEON(a *float64, n int) float64 diff --git a/internal/kernels/reduce_floats.go b/internal/kernels/reduce_floats.go new file mode 100644 index 0000000..b265220 --- /dev/null +++ b/internal/kernels/reduce_floats.go @@ -0,0 +1,51 @@ +//go:build amd64 || arm64 + +package kernels + +import "github.com/go-simd/floats" + +// On amd64 (AVX2) and arm64 (NEON) the elementwise, square-root and extremum +// kernels are go-simd/floats': K1 of the go-scikit performance campaign moved +// them there, so that a kernel lives once, in the lowest library every user can +// import. Each is bit-identical to the scalar oracle (Add, Sub, Mul, Div, +// sqrtScalar, Max, Min of kernels.go) on every input, signed zeros and NaNs +// included; the tests of this package hold them to it. +// +// amd64 needs AVX2+FMA for them (floats' gate); on a CPU without it floats runs +// its portable loops, which the compiler does not vectorise (the SSE2 kernels +// this package had before ran there). + +// HaveReduceSIMD reports that this build ships SIMD reduction/elementwise +// kernels (so the dispatch wrappers route through them rather than the scalar +// oracle). True on amd64/arm64. +const HaveReduceSIMD = true + +// sqrtSIMD writes sqrt(src[i]) into dst[i] (correctly rounded, like math.Sqrt). +func sqrtSIMD(dst, src []float64) { floats.Sqrt(dst, src[:len(dst)]) } + +// maxSIMD returns the NaN-propagating maximum of a (non-empty), the serial Max +// oracle's value (max(-0, +0) = +0, as NumPy). floats.Max leaves the sign of a +// zero result unspecified (its kernels compare zeros as equal), so a zero result, +// the one case where it matters, is resolved by the oracle, which costs nothing +// on any other result. +func maxSIMD(a []float64) float64 { + if r := floats.Max(a); r != 0 { + return r + } + return Max(a) +} + +// minSIMD returns the NaN-propagating minimum of a (non-empty) (min(-0, +0) = +// -0), resolved like maxSIMD. +func minSIMD(a []float64) float64 { + if r := floats.Min(a); r != 0 { + return r + } + return Min(a) +} + +// addBin/subBin/mulBin/divBin write a[i] OP b[i] into dst[i]; len(dst) elements. +func addBin(dst, a, b []float64) { floats.Add(dst, a[:len(dst)], b[:len(dst)]) } +func subBin(dst, a, b []float64) { floats.Sub(dst, a[:len(dst)], b[:len(dst)]) } +func mulBin(dst, a, b []float64) { floats.Mul(dst, a[:len(dst)], b[:len(dst)]) } +func divBin(dst, a, b []float64) { floats.Div(dst, a[:len(dst)], b[:len(dst)]) } diff --git a/internal/kernels/reduce_generic.go b/internal/kernels/reduce_generic.go index f181c7c..9251cfc 100644 --- a/internal/kernels/reduce_generic.go +++ b/internal/kernels/reduce_generic.go @@ -19,7 +19,6 @@ package kernels // run states which path it validated. const HaveReduceSIMD = false -func sumSIMD(a []float64) float64 { return Sum(a) } func sqrtSIMD(dst, src []float64) { sqrtScalar(dst, src) } func maxSIMD(a []float64) float64 { return maxUnrolled(a) } func minSIMD(a []float64) float64 { return minUnrolled(a) } diff --git a/internal/kernels/reduce_loong64.go b/internal/kernels/reduce_loong64.go index 59f182f..f82cf6f 100644 --- a/internal/kernels/reduce_loong64.go +++ b/internal/kernels/reduce_loong64.go @@ -28,12 +28,8 @@ func lasxFromAuxv() bool { return err == nil && hwcapFrom(auxv)&hwcapLASX != 0 } -func sumSIMD(a []float64) float64 { - if !useLASX || len(a) == 0 { - return Sum(a) - } - return sumLASX(&a[0], len(a)) -} +// The sum is floats.SumPairwise on every target (reduce_sum.go): NumPy's order. +// sumLASX stays in sum_loong64.s until the loong64 kernels move to go-simd/floats. func sqrtSIMD(dst, src []float64) { if !useLASX || len(dst) == 0 { diff --git a/internal/kernels/reduce_ppc64le.go b/internal/kernels/reduce_ppc64le.go index b91e779..95b1c3b 100644 --- a/internal/kernels/reduce_ppc64le.go +++ b/internal/kernels/reduce_ppc64le.go @@ -13,14 +13,8 @@ import "math" // kernels. const HaveReduceSIMD = true -// sumSIMD returns the VSX sum of a (empty -> 0): four accumulators, a valid -// regrouping of the sum, like the other targets' SIMD sums. -func sumSIMD(a []float64) float64 { - if len(a) == 0 { - return 0 - } - return sumVSX(&a[0], len(a)) -} +// The sum is floats.SumPairwise on every target (reduce_sum.go): NumPy's order. +// sumVSX stays in sum_ppc64le.s until the ppc64le kernels move to go-simd/floats. // sqrtSIMD writes sqrt(src[i]) into dst[i] with XVSQRTDP, bit-identical to // math.Sqrt (one correctly rounded operation per lane). diff --git a/internal/kernels/reduce_riscv64.go b/internal/kernels/reduce_riscv64.go index cf39bde..d4a0dd3 100644 --- a/internal/kernels/reduce_riscv64.go +++ b/internal/kernels/reduce_riscv64.go @@ -32,19 +32,10 @@ func rvvFromAuxv() bool { return err == nil && hwcapFrom(auxv)&hwcapV != 0 } -// sumSIMD is the RVV sum from 8 elements up. Below that it is the sequential -// sum, the order NumPy's pairwise sum and the other targets' kernels (whose -// scalar tails take n < 8 whole) use there. The strip-mined kernel would take -// even 3 elements as one vector, and its unordered reduction regrouped -// [MaxInt64, MinInt64, -3] as float64 so that max + (-3) rounded back to max -// and cancelled: the mean was 0, NumPy's -1. Seen on a real SpacemiT X60, -// never under CI's qemu, whose default CPU has no V. -func sumSIMD(a []float64) float64 { - if !useRVV || len(a) < 8 { - return Sum(a) - } - return sumRVV(&a[0], len(a)) -} +// The sum is floats.SumPairwise on every target (reduce_sum.go): NumPy's order. +// sumRVV stays in sum_riscv64.s until the riscv64 kernels move to go-simd/floats +// (it regrouped [MaxInt64, MinInt64, -3] as float64 on a real SpacemiT X60: the +// mean was 0, NumPy's -1; SumPairwise adds n < 8 left to right as NumPy does). func sqrtSIMD(dst, src []float64) { if !useRVV || len(dst) == 0 { diff --git a/internal/kernels/reduce_s390x.go b/internal/kernels/reduce_s390x.go index a01804f..73f4299 100644 --- a/internal/kernels/reduce_s390x.go +++ b/internal/kernels/reduce_s390x.go @@ -30,12 +30,8 @@ func vxFromAuxv() bool { return err == nil && hwcapFrom(auxv)&hwcapVX != 0 } -func sumSIMD(a []float64) float64 { - if !useVX || len(a) == 0 { - return Sum(a) - } - return sumVX(&a[0], len(a)) -} +// The sum is floats.SumPairwise on every target (reduce_sum.go): NumPy's order. +// sumVX stays in sum_s390x.s until the s390x kernels move to go-simd/floats. func sqrtSIMD(dst, src []float64) { if !useVX || len(dst) == 0 { diff --git a/internal/kernels/reduce_sum.go b/internal/kernels/reduce_sum.go new file mode 100644 index 0000000..7b289f6 --- /dev/null +++ b/internal/kernels/reduce_sum.go @@ -0,0 +1,12 @@ +package kernels + +import "github.com/go-simd/floats" + +// sumSIMD returns the sum of a (empty -> 0) in NumPy's order: floats.SumPairwise, +// bit-identical to np.sum of a contiguous 1-D float64 array on every architecture +// (K1 of the go-scikit performance campaign: the kernel lives in go-simd/floats). +// +// This replaces the lane-parallel sums of the per-architecture kernels, whose +// grouping differed from NumPy's (and from one architecture to the next) by a few +// ULP. NumPy's order is the reference, so a changed last bit is the fix. +func sumSIMD(a []float64) float64 { return floats.SumPairwise(a) } diff --git a/internal/kernels/sum_amd64.s b/internal/kernels/sum_amd64.s index 2df024c..2a85df3 100644 --- a/internal/kernels/sum_amd64.s +++ b/internal/kernels/sum_amd64.s @@ -4,395 +4,6 @@ #include "textflag.h" -TEXT ·sumSSE2(SB), NOSPLIT, $0-24 - MOVQ a+0(FP), AX - MOVQ n+8(FP), CX - XORPS X0, X0 - XORPS X1, X1 - XORPS X2, X2 - XORPS X3, X3 - block: - CMPQ CX, $8 - JL fold - MOVUPD (AX), X4 - MOVUPD 16(AX), X5 - MOVUPD 32(AX), X6 - MOVUPD 48(AX), X7 - ADDPD X4, X0 - ADDPD X5, X1 - ADDPD X6, X2 - ADDPD X7, X3 - ADDQ $64, AX - SUBQ $8, CX - JMP block - fold: - ADDPD X1, X0 - ADDPD X3, X2 - ADDPD X2, X0 - MOVAPS X0, X1 - UNPCKHPD X1, X1 - ADDSD X1, X0 - tail: - TESTQ CX, CX - JZ done - ADDSD (AX), X0 - ADDQ $8, AX - DECQ CX - JMP tail - done: - MOVSD X0, ret+16(FP) - RET - -TEXT ·sqrtSSE2(SB), NOSPLIT, $0-24 - MOVQ dst+0(FP), DI - MOVQ src+8(FP), SI - MOVQ n+16(FP), CX - sblock: - CMPQ CX, $8 - JL stail - MOVUPD (SI), X0 - MOVUPD 16(SI), X1 - MOVUPD 32(SI), X2 - MOVUPD 48(SI), X3 - SQRTPD X0, X0 - SQRTPD X1, X1 - SQRTPD X2, X2 - SQRTPD X3, X3 - MOVUPD X0, (DI) - MOVUPD X1, 16(DI) - MOVUPD X2, 32(DI) - MOVUPD X3, 48(DI) - ADDQ $64, SI - ADDQ $64, DI - SUBQ $8, CX - JMP sblock - stail: - TESTQ CX, CX - JZ sdone - MOVSD (SI), X0 - SQRTSD X0, X0 - MOVSD X0, (DI) - ADDQ $8, SI - ADDQ $8, DI - DECQ CX - JMP stail - sdone: - RET - -TEXT ·maxSSE2(SB), NOSPLIT, $0-24 - MOVQ a+0(FP), AX - MOVQ n+8(FP), CX - MOVDDUP (AX), X0 - MOVAPS X0, X1 - MOVAPS X0, X2 - MOVAPS X0, X3 - XORPS X8, X8 - MOVAPS X0, X9 - CMPPD X9, X9, $3 - ORPD X9, X8 - eblock: - CMPQ CX, $8 - JL efold - MOVUPD (AX), X4 - MOVUPD 16(AX), X5 - MOVUPD 32(AX), X6 - MOVUPD 48(AX), X7 - MOVAPS X4, X9 - CMPPD X9, X9, $3 - ORPD X9, X8 - MOVAPS X5, X9 - CMPPD X9, X9, $3 - ORPD X9, X8 - MOVAPS X6, X9 - CMPPD X9, X9, $3 - ORPD X9, X8 - MOVAPS X7, X9 - CMPPD X9, X9, $3 - ORPD X9, X8 - MAXPD X4, X0 - MAXPD X5, X1 - MAXPD X6, X2 - MAXPD X7, X3 - ADDQ $64, AX - SUBQ $8, CX - JMP eblock - efold: - MAXPD X1, X0 - MAXPD X3, X2 - MAXPD X2, X0 - MOVAPS X0, X1 - UNPCKHPD X1, X1 - MAXSD X1, X0 - etail: - TESTQ CX, CX - JZ enan - MOVSD (AX), X4 - MOVAPS X4, X9 - CMPPD X9, X9, $3 - ORPD X9, X8 - MAXSD X4, X0 - ADDQ $8, AX - DECQ CX - JMP etail - enan: - MOVAPS X8, X9 - UNPCKHPD X9, X9 - ORPD X9, X8 - MOVQ X8, DX - TESTQ DX, DX - JZ edone - XORPS X1, X1 - DIVSD X1, X1 - MOVSD X1, X0 - edone: - MOVSD X0, ret+16(FP) - RET - -TEXT ·minSSE2(SB), NOSPLIT, $0-24 - MOVQ a+0(FP), AX - MOVQ n+8(FP), CX - MOVDDUP (AX), X0 - MOVAPS X0, X1 - MOVAPS X0, X2 - MOVAPS X0, X3 - XORPS X8, X8 - MOVAPS X0, X9 - CMPPD X9, X9, $3 - ORPD X9, X8 - eblock: - CMPQ CX, $8 - JL efold - MOVUPD (AX), X4 - MOVUPD 16(AX), X5 - MOVUPD 32(AX), X6 - MOVUPD 48(AX), X7 - MOVAPS X4, X9 - CMPPD X9, X9, $3 - ORPD X9, X8 - MOVAPS X5, X9 - CMPPD X9, X9, $3 - ORPD X9, X8 - MOVAPS X6, X9 - CMPPD X9, X9, $3 - ORPD X9, X8 - MOVAPS X7, X9 - CMPPD X9, X9, $3 - ORPD X9, X8 - MINPD X4, X0 - MINPD X5, X1 - MINPD X6, X2 - MINPD X7, X3 - ADDQ $64, AX - SUBQ $8, CX - JMP eblock - efold: - MINPD X1, X0 - MINPD X3, X2 - MINPD X2, X0 - MOVAPS X0, X1 - UNPCKHPD X1, X1 - MINSD X1, X0 - etail: - TESTQ CX, CX - JZ enan - MOVSD (AX), X4 - MOVAPS X4, X9 - CMPPD X9, X9, $3 - ORPD X9, X8 - MINSD X4, X0 - ADDQ $8, AX - DECQ CX - JMP etail - enan: - MOVAPS X8, X9 - UNPCKHPD X9, X9 - ORPD X9, X8 - MOVQ X8, DX - TESTQ DX, DX - JZ edone - XORPS X1, X1 - DIVSD X1, X1 - MOVSD X1, X0 - edone: - MOVSD X0, ret+16(FP) - RET - -TEXT ·addSSE2(SB), NOSPLIT, $0-32 - MOVQ dst+0(FP), DI - MOVQ a+8(FP), SI - MOVQ b+16(FP), DX - MOVQ n+24(FP), CX - bblock: - CMPQ CX, $8 - JL btail - MOVUPD (SI), X0 - MOVUPD 16(SI), X1 - MOVUPD 32(SI), X2 - MOVUPD 48(SI), X3 - MOVUPD (DX), X4 - MOVUPD 16(DX), X5 - MOVUPD 32(DX), X6 - MOVUPD 48(DX), X7 - ADDPD X4, X0 - ADDPD X5, X1 - ADDPD X6, X2 - ADDPD X7, X3 - MOVUPD X0, (DI) - MOVUPD X1, 16(DI) - MOVUPD X2, 32(DI) - MOVUPD X3, 48(DI) - ADDQ $64, SI - ADDQ $64, DX - ADDQ $64, DI - SUBQ $8, CX - JMP bblock - btail: - TESTQ CX, CX - JZ bdone - MOVSD (SI), X0 - MOVSD (DX), X1 - ADDSD X1, X0 - MOVSD X0, (DI) - ADDQ $8, SI - ADDQ $8, DX - ADDQ $8, DI - DECQ CX - JMP btail - bdone: - RET - -TEXT ·subSSE2(SB), NOSPLIT, $0-32 - MOVQ dst+0(FP), DI - MOVQ a+8(FP), SI - MOVQ b+16(FP), DX - MOVQ n+24(FP), CX - bblock: - CMPQ CX, $8 - JL btail - MOVUPD (SI), X0 - MOVUPD 16(SI), X1 - MOVUPD 32(SI), X2 - MOVUPD 48(SI), X3 - MOVUPD (DX), X4 - MOVUPD 16(DX), X5 - MOVUPD 32(DX), X6 - MOVUPD 48(DX), X7 - SUBPD X4, X0 - SUBPD X5, X1 - SUBPD X6, X2 - SUBPD X7, X3 - MOVUPD X0, (DI) - MOVUPD X1, 16(DI) - MOVUPD X2, 32(DI) - MOVUPD X3, 48(DI) - ADDQ $64, SI - ADDQ $64, DX - ADDQ $64, DI - SUBQ $8, CX - JMP bblock - btail: - TESTQ CX, CX - JZ bdone - MOVSD (SI), X0 - MOVSD (DX), X1 - SUBSD X1, X0 - MOVSD X0, (DI) - ADDQ $8, SI - ADDQ $8, DX - ADDQ $8, DI - DECQ CX - JMP btail - bdone: - RET - -TEXT ·mulSSE2(SB), NOSPLIT, $0-32 - MOVQ dst+0(FP), DI - MOVQ a+8(FP), SI - MOVQ b+16(FP), DX - MOVQ n+24(FP), CX - bblock: - CMPQ CX, $8 - JL btail - MOVUPD (SI), X0 - MOVUPD 16(SI), X1 - MOVUPD 32(SI), X2 - MOVUPD 48(SI), X3 - MOVUPD (DX), X4 - MOVUPD 16(DX), X5 - MOVUPD 32(DX), X6 - MOVUPD 48(DX), X7 - MULPD X4, X0 - MULPD X5, X1 - MULPD X6, X2 - MULPD X7, X3 - MOVUPD X0, (DI) - MOVUPD X1, 16(DI) - MOVUPD X2, 32(DI) - MOVUPD X3, 48(DI) - ADDQ $64, SI - ADDQ $64, DX - ADDQ $64, DI - SUBQ $8, CX - JMP bblock - btail: - TESTQ CX, CX - JZ bdone - MOVSD (SI), X0 - MOVSD (DX), X1 - MULSD X1, X0 - MOVSD X0, (DI) - ADDQ $8, SI - ADDQ $8, DX - ADDQ $8, DI - DECQ CX - JMP btail - bdone: - RET - -TEXT ·divSSE2(SB), NOSPLIT, $0-32 - MOVQ dst+0(FP), DI - MOVQ a+8(FP), SI - MOVQ b+16(FP), DX - MOVQ n+24(FP), CX - bblock: - CMPQ CX, $8 - JL btail - MOVUPD (SI), X0 - MOVUPD 16(SI), X1 - MOVUPD 32(SI), X2 - MOVUPD 48(SI), X3 - MOVUPD (DX), X4 - MOVUPD 16(DX), X5 - MOVUPD 32(DX), X6 - MOVUPD 48(DX), X7 - DIVPD X4, X0 - DIVPD X5, X1 - DIVPD X6, X2 - DIVPD X7, X3 - MOVUPD X0, (DI) - MOVUPD X1, 16(DI) - MOVUPD X2, 32(DI) - MOVUPD X3, 48(DI) - ADDQ $64, SI - ADDQ $64, DX - ADDQ $64, DI - SUBQ $8, CX - JMP bblock - btail: - TESTQ CX, CX - JZ bdone - MOVSD (SI), X0 - MOVSD (DX), X1 - DIVSD X1, X0 - MOVSD X0, (DI) - ADDQ $8, SI - ADDQ $8, DX - ADDQ $8, DI - DECQ CX - JMP btail - bdone: - RET - TEXT ·gemmMicro6x8FMA(SB), NOSPLIT, $0-40 MOVQ kc+0(FP), CX MOVQ pa+8(FP), SI @@ -1692,247 +1303,6 @@ TEXT ·hasFMA(SB), NOSPLIT, $0-1 hasFMA_unsupported: RET -TEXT ·hasAVX2(SB), NOSPLIT, $0-1 - MOVB $0, ret+0(FP) - XORL AX, AX - CPUID - CMPL AX, $7 - JL hasAVX2_unsupported - MOVL $1, AX - CPUID - ANDL $0x18000000, CX // AVX and OSXSAVE - CMPL CX, $0x18000000 - JNE hasAVX2_unsupported - XORL CX, CX - XGETBV - ANDL $6, AX // XMM and YMM state enabled in XCR0 - CMPL AX, $6 - JNE hasAVX2_unsupported - MOVL $7, AX - XORL CX, CX - CPUID - TESTL $0x20, BX // AVX2 - JZ hasAVX2_unsupported - MOVB $1, ret+0(FP) -hasAVX2_unsupported: - RET - -TEXT ·exp32CRAVX2(SB), NOSPLIT, $0-32 - MOVQ dst+0(FP), DI - MOVQ src+8(FP), SI - MOVQ n+16(FP), CX - MOVQ k+24(FP), R8 - LEAQ 384(R8), R9 - eloop: - VMOVUPS (SI), Y0 - VCMPPS $3, Y0, Y0, Y15 - VMAXPS 320(R8), Y0, Y1 - VMINPS 352(R8), Y1, Y1 - VCVTPS2PD X1, Y2 - VEXTRACTF128 $1, Y1, X3 - VCVTPS2PD X3, Y3 - VMOVUPD 0(R8), Y4 - VFMADD231PD 32(R8), Y2, Y4 - VSUBPD 0(R8), Y4, Y5 - VFNMADD231PD 64(R8), Y5, Y2 - VFNMADD231PD 96(R8), Y5, Y2 - VMULPD Y2, Y2, Y6 - VMOVUPD 128(R8), Y7 - VFMADD231PD 160(R8), Y2, Y7 - VMOVUPD 192(R8), Y8 - VFMADD231PD 224(R8), Y2, Y8 - VFMADD231PD Y8, Y6, Y7 - VFMADD231PD Y7, Y6, Y2 - VPAND 256(R8), Y4, Y5 - VPCMPEQQ Y9, Y9, Y9 - VPGATHERQQ Y9, (R9)(Y5*8), Y10 - VPSLLQ $46, Y4, Y4 - VPAND 288(R8), Y4, Y4 - VPOR Y10, Y4, Y4 - VFMADD231PD Y2, Y4, Y4 - VCVTPD2PSY Y4, X11 - VMOVUPD 0(R8), Y4 - VFMADD231PD 32(R8), Y3, Y4 - VSUBPD 0(R8), Y4, Y5 - VFNMADD231PD 64(R8), Y5, Y3 - VFNMADD231PD 96(R8), Y5, Y3 - VMULPD Y3, Y3, Y6 - VMOVUPD 128(R8), Y7 - VFMADD231PD 160(R8), Y3, Y7 - VMOVUPD 192(R8), Y8 - VFMADD231PD 224(R8), Y3, Y8 - VFMADD231PD Y8, Y6, Y7 - VFMADD231PD Y7, Y6, Y3 - VPAND 256(R8), Y4, Y5 - VPCMPEQQ Y9, Y9, Y9 - VPGATHERQQ Y9, (R9)(Y5*8), Y10 - VPSLLQ $46, Y4, Y4 - VPAND 288(R8), Y4, Y4 - VPOR Y10, Y4, Y4 - VFMADD231PD Y3, Y4, Y4 - VCVTPD2PSY Y4, X12 - VINSERTF128 $1, X12, Y11, Y11 - VBLENDVPS Y15, Y0, Y11, Y11 - VMOVUPS Y11, (DI) - ADDQ $32, SI - ADDQ $32, DI - SUBQ $8, CX - JNZ eloop - VZEROUPPER - RET - -TEXT ·log32CRAVX2(SB), NOSPLIT, $0-40 - MOVQ dst+0(FP), DI - MOVQ src+8(FP), SI - MOVQ n+16(FP), CX - MOVQ k+24(FP), R8 - LEAQ 416(R8), R9 - LEAQ 936(R8), R10 - XORQ AX, AX - lloop: - VMOVDQU (SI), Y0 - VPCMPEQD 352(R8), Y0, Y13 - VPAND 256(R8), Y0, Y1 - VPADDD 288(R8), Y1, Y2 - VPSRLD $17, Y2, Y2 - VPSRLD $23, Y0, Y3 - VPSUBD 320(R8), Y3, Y3 - VPMOVZXDQ X1, Y10 - VCVTDQ2PD X3, Y11 - VPOR 384(R8), Y10, Y10 - VPSLLQ $29, Y10, Y10 - VPCMPEQQ Y12, Y12, Y12 - VGATHERDPD Y12, (R9)(X2*8), Y4 - VPCMPEQQ Y12, Y12, Y12 - VGATHERDPD Y12, (R10)(X2*8), Y5 - VMULPD Y4, Y10, Y4 - VSUBPD 0(R8), Y4, Y4 - VMULPD Y4, Y4, Y6 - VMULPD 32(R8), Y11, Y7 - VADDPD Y5, Y7, Y7 - VMULPD 64(R8), Y4, Y8 - VADDPD Y8, Y7, Y7 - VMULPD 128(R8), Y4, Y8 - VADDPD 96(R8), Y8, Y8 - VMULPD Y8, Y6, Y8 - VADDPD Y8, Y7, Y7 - VADDPD 160(R8), Y7, Y9 - VCVTPD2PSY Y7, X14 - VCVTPD2PSY Y9, X15 - VEXTRACTI128 $1, Y1, X1 - VPMOVZXDQ X1, Y10 - VEXTRACTI128 $1, Y3, X3 - VCVTDQ2PD X3, Y11 - VEXTRACTI128 $1, Y2, X2 - VPOR 384(R8), Y10, Y10 - VPSLLQ $29, Y10, Y10 - VPCMPEQQ Y12, Y12, Y12 - VGATHERDPD Y12, (R9)(X2*8), Y4 - VPCMPEQQ Y12, Y12, Y12 - VGATHERDPD Y12, (R10)(X2*8), Y5 - VMULPD Y4, Y10, Y4 - VSUBPD 0(R8), Y4, Y4 - VMULPD Y4, Y4, Y6 - VMULPD 32(R8), Y11, Y7 - VADDPD Y5, Y7, Y7 - VMULPD 64(R8), Y4, Y8 - VADDPD Y8, Y7, Y7 - VMULPD 128(R8), Y4, Y8 - VADDPD 96(R8), Y8, Y8 - VMULPD Y8, Y6, Y8 - VADDPD Y8, Y7, Y7 - VADDPD 160(R8), Y7, Y9 - VCVTPD2PSY Y7, X1 - VCVTPD2PSY Y9, X3 - VINSERTI128 $1, X1, Y14, Y14 - VINSERTI128 $1, X3, Y15, Y15 - VPCMPEQD Y15, Y14, Y15 - VPOR Y13, Y15, Y15 - VPADDD 192(R8), Y0, Y4 - VPCMPGTD 224(R8), Y4, Y4 - VPANDN Y15, Y4, Y15 - VMOVMSKPS Y15, DX - CMPL DX, $0xff - JNE lbail - VPANDN Y14, Y13, Y14 - VMOVDQU Y14, (DI) - ADDQ $32, SI - ADDQ $32, DI - ADDQ $8, AX - CMPQ AX, CX - JLT lloop - lbail: - VZEROUPPER - MOVQ AX, done+32(FP) - RET - -TEXT ·logAVX2(SB), NOSPLIT, $0-48 - MOVQ dst+0(FP), DI - MOVQ src+8(FP), SI - MOVQ n+16(FP), CX - MOVQ k+24(FP), R8 - MOVQ tab+32(FP), R9 - XORQ AX, AX - dloop: - VMOVDQU (SI), Y0 - VPCMPGTQ 32(R8), Y0, Y1 - VMOVDQU 0(R8), Y2 - VPCMPGTQ Y0, Y2, Y2 - VPCMPGTQ 64(R8), Y0, Y3 - VMOVDQU 96(R8), Y4 - VPCMPGTQ Y0, Y4, Y4 - VPAND Y4, Y3, Y3 - VPOR Y3, Y2, Y2 - VPOR Y1, Y2, Y2 - VPTEST Y2, Y2 - JNZ dbail - VPSUBQ 128(R8), Y0, Y1 - VPSRLQ $45, Y1, Y2 - VPAND 160(R8), Y2, Y2 - VPADDQ Y2, Y2, Y2 - VPSRLQ $52, Y1, Y3 - VPXOR 288(R8), Y3, Y3 - VPOR 224(R8), Y3, Y3 - VSUBPD 256(R8), Y3, Y3 - VPAND 192(R8), Y1, Y4 - VPSUBQ Y4, Y0, Y4 - VPCMPEQQ Y5, Y5, Y5 - VPGATHERQQ Y5, (R9)(Y2*8), Y6 - VPCMPEQQ Y5, Y5, Y5 - VPGATHERQQ Y5, 8(R9)(Y2*8), Y7 - VMOVUPD 320(R8), Y8 - VFMADD231PD Y6, Y4, Y8 - VMULPD 352(R8), Y3, Y9 - VADDPD Y7, Y9, Y9 - VADDPD Y8, Y9, Y10 - VSUBPD Y10, Y9, Y11 - VADDPD Y8, Y11, Y11 - VMULPD 384(R8), Y3, Y12 - VADDPD Y12, Y11, Y11 - VMULPD Y8, Y8, Y12 - VMULPD 480(R8), Y8, Y13 - VADDPD 448(R8), Y13, Y13 - VMULPD 544(R8), Y8, Y14 - VADDPD 512(R8), Y14, Y14 - VMULPD Y14, Y12, Y14 - VADDPD Y14, Y13, Y13 - VMULPD Y12, Y8, Y14 - VMULPD Y13, Y14, Y14 - VMULPD 416(R8), Y12, Y13 - VADDPD Y13, Y11, Y11 - VADDPD Y14, Y11, Y11 - VADDPD Y10, Y11, Y11 - VMOVUPD Y11, (DI) - ADDQ $32, SI - ADDQ $32, DI - ADDQ $4, AX - CMPQ AX, CX - JLT dloop - dbail: - VZEROUPPER - MOVQ AX, done+40(FP) - RET - TEXT ·blocks8SSE2(SB), NOSPLIT, $0-24 MOVQ r+0(FP), DI MOVQ a+8(FP), SI diff --git a/internal/kernels/sum_arm64.s b/internal/kernels/sum_arm64.s index 1075603..c856c6a 100644 --- a/internal/kernels/sum_arm64.s +++ b/internal/kernels/sum_arm64.s @@ -4,44 +4,6 @@ #include "textflag.h" -TEXT ·sumNEON(SB), NOSPLIT, $0-24 - MOVD a+0(FP), R0 - MOVD n+8(FP), R1 - MOVD $0x3FF0000000000000, R2 - VDUP R2, V16.D2 - VEOR V0.B16, V0.B16, V0.B16 - VEOR V1.B16, V1.B16, V1.B16 - VEOR V2.B16, V2.B16, V2.B16 - VEOR V3.B16, V3.B16, V3.B16 - block: - CMP $8, R1 - BLT fold - VLD1.P 64(R0), [V4.D2, V5.D2, V6.D2, V7.D2] - VFMLA V4.D2, V16.D2, V0.D2 - VFMLA V5.D2, V16.D2, V1.D2 - VFMLA V6.D2, V16.D2, V2.D2 - VFMLA V7.D2, V16.D2, V3.D2 - SUB $8, R1 - B block - fold: - VFMLA V1.D2, V16.D2, V0.D2 - VFMLA V3.D2, V16.D2, V2.D2 - VFMLA V2.D2, V16.D2, V0.D2 - VMOV V0.D[0], R3 - VMOV V0.D[1], R4 - FMOVD R3, F0 - FMOVD R4, F9 - FADDD F9, F0, F0 - tail: - CBZ R1, done - FMOVD.P 8(R0), F10 - FADDD F10, F0, F0 - SUB $1, R1 - B tail - done: - FMOVD F0, ret+16(FP) - RET - TEXT ·dotNEON(SB), NOSPLIT, $0-32 MOVD a+0(FP), R0 MOVD b+8(FP), R1 @@ -97,143 +59,6 @@ TEXT ·dotNEON(SB), NOSPLIT, $0-32 FMOVD F0, ret+24(FP) RET -TEXT ·sqrtNEON(SB), NOSPLIT, $0-24 - MOVD dst+0(FP), R0 - MOVD src+8(FP), R1 - MOVD n+16(FP), R2 - block: - CMP $8, R2 - BLT tail - VLD1.P 64(R1), [V0.D2, V1.D2, V2.D2, V3.D2] - VFSQRT V0.D2, V0.D2 - VFSQRT V1.D2, V1.D2 - VFSQRT V2.D2, V2.D2 - VFSQRT V3.D2, V3.D2 - VST1.P [V0.D2, V1.D2, V2.D2, V3.D2], 64(R0) - SUB $8, R2 - B block - tail: - CBZ R2, done - FMOVD.P 8(R1), F0 - FSQRTD F0, F0 - FMOVD.P F0, 8(R0) - SUB $1, R2 - B tail - done: - RET - -TEXT ·addNEON(SB), NOSPLIT, $0-32 - MOVD dst+0(FP), R0 - MOVD a+8(FP), R1 - MOVD b+16(FP), R2 - MOVD n+24(FP), R3 - block: - CMP $8, R3 - BLT tail - VLD1.P 64(R1), [V0.D2, V1.D2, V2.D2, V3.D2] - VLD1.P 64(R2), [V4.D2, V5.D2, V6.D2, V7.D2] - VFADD V4.D2, V0.D2, V0.D2 - VFADD V5.D2, V1.D2, V1.D2 - VFADD V6.D2, V2.D2, V2.D2 - VFADD V7.D2, V3.D2, V3.D2 - VST1.P [V0.D2, V1.D2, V2.D2, V3.D2], 64(R0) - SUB $8, R3 - B block - tail: - CBZ R3, done - FMOVD.P 8(R1), F0 - FMOVD.P 8(R2), F1 - FADDD F1, F0, F0 - FMOVD.P F0, 8(R0) - SUB $1, R3 - B tail - done: - RET - -TEXT ·subNEON(SB), NOSPLIT, $0-32 - MOVD dst+0(FP), R0 - MOVD a+8(FP), R1 - MOVD b+16(FP), R2 - MOVD n+24(FP), R3 - block: - CMP $8, R3 - BLT tail - VLD1.P 64(R1), [V0.D2, V1.D2, V2.D2, V3.D2] - VLD1.P 64(R2), [V4.D2, V5.D2, V6.D2, V7.D2] - VFSUB V4.D2, V0.D2, V0.D2 - VFSUB V5.D2, V1.D2, V1.D2 - VFSUB V6.D2, V2.D2, V2.D2 - VFSUB V7.D2, V3.D2, V3.D2 - VST1.P [V0.D2, V1.D2, V2.D2, V3.D2], 64(R0) - SUB $8, R3 - B block - tail: - CBZ R3, done - FMOVD.P 8(R1), F0 - FMOVD.P 8(R2), F1 - FSUBD F1, F0, F0 - FMOVD.P F0, 8(R0) - SUB $1, R3 - B tail - done: - RET - -TEXT ·mulNEON(SB), NOSPLIT, $0-32 - MOVD dst+0(FP), R0 - MOVD a+8(FP), R1 - MOVD b+16(FP), R2 - MOVD n+24(FP), R3 - block: - CMP $8, R3 - BLT tail - VLD1.P 64(R1), [V0.D2, V1.D2, V2.D2, V3.D2] - VLD1.P 64(R2), [V4.D2, V5.D2, V6.D2, V7.D2] - VFMUL V4.D2, V0.D2, V0.D2 - VFMUL V5.D2, V1.D2, V1.D2 - VFMUL V6.D2, V2.D2, V2.D2 - VFMUL V7.D2, V3.D2, V3.D2 - VST1.P [V0.D2, V1.D2, V2.D2, V3.D2], 64(R0) - SUB $8, R3 - B block - tail: - CBZ R3, done - FMOVD.P 8(R1), F0 - FMOVD.P 8(R2), F1 - FMULD F1, F0, F0 - FMOVD.P F0, 8(R0) - SUB $1, R3 - B tail - done: - RET - -TEXT ·divNEON(SB), NOSPLIT, $0-32 - MOVD dst+0(FP), R0 - MOVD a+8(FP), R1 - MOVD b+16(FP), R2 - MOVD n+24(FP), R3 - block: - CMP $8, R3 - BLT tail - VLD1.P 64(R1), [V0.D2, V1.D2, V2.D2, V3.D2] - VLD1.P 64(R2), [V4.D2, V5.D2, V6.D2, V7.D2] - VFDIV V4.D2, V0.D2, V0.D2 - VFDIV V5.D2, V1.D2, V1.D2 - VFDIV V6.D2, V2.D2, V2.D2 - VFDIV V7.D2, V3.D2, V3.D2 - VST1.P [V0.D2, V1.D2, V2.D2, V3.D2], 64(R0) - SUB $8, R3 - B block - tail: - CBZ R3, done - FMOVD.P 8(R1), F0 - FMOVD.P 8(R2), F1 - FDIVD F1, F0, F0 - FMOVD.P F0, 8(R0) - SUB $1, R3 - B tail - done: - RET - TEXT ·add32NEON(SB), NOSPLIT, $0-32 MOVD dst+0(FP), R0 MOVD a+8(FP), R1 @@ -371,74 +196,6 @@ TEXT ·sqrt32NEON(SB), NOSPLIT, $0-24 done: RET -TEXT ·maxNEON(SB), NOSPLIT, $0-24 - MOVD a+0(FP), R0 - MOVD n+8(FP), R1 - VLD1.P 64(R0), [V0.D2, V1.D2, V2.D2, V3.D2] - SUB $8, R1 - block: - CMP $8, R1 - BLT fold - VLD1.P 64(R0), [V4.D2, V5.D2, V6.D2, V7.D2] - VFMAX V4.D2, V0.D2, V0.D2 - VFMAX V5.D2, V1.D2, V1.D2 - VFMAX V6.D2, V2.D2, V2.D2 - VFMAX V7.D2, V3.D2, V3.D2 - SUB $8, R1 - B block - fold: - VFMAX V1.D2, V0.D2, V0.D2 - VFMAX V3.D2, V2.D2, V2.D2 - VFMAX V2.D2, V0.D2, V0.D2 - VMOV V0.D[0], R3 - VMOV V0.D[1], R4 - FMOVD R3, F0 - FMOVD R4, F9 - FMAXD F9, F0, F0 - tail: - CBZ R1, done - FMOVD.P 8(R0), F10 - FMAXD F10, F0, F0 - SUB $1, R1 - B tail - done: - FMOVD F0, ret+16(FP) - RET - -TEXT ·minNEON(SB), NOSPLIT, $0-24 - MOVD a+0(FP), R0 - MOVD n+8(FP), R1 - VLD1.P 64(R0), [V0.D2, V1.D2, V2.D2, V3.D2] - SUB $8, R1 - block: - CMP $8, R1 - BLT fold - VLD1.P 64(R0), [V4.D2, V5.D2, V6.D2, V7.D2] - VFMIN V4.D2, V0.D2, V0.D2 - VFMIN V5.D2, V1.D2, V1.D2 - VFMIN V6.D2, V2.D2, V2.D2 - VFMIN V7.D2, V3.D2, V3.D2 - SUB $8, R1 - B block - fold: - VFMIN V1.D2, V0.D2, V0.D2 - VFMIN V3.D2, V2.D2, V2.D2 - VFMIN V2.D2, V0.D2, V0.D2 - VMOV V0.D[0], R3 - VMOV V0.D[1], R4 - FMOVD R3, F0 - FMOVD R4, F9 - FMIND F9, F0, F0 - tail: - CBZ R1, done - FMOVD.P 8(R0), F10 - FMIND F10, F0, F0 - SUB $1, R1 - B tail - done: - FMOVD F0, ret+16(FP) - RET - TEXT ·max32NEON(SB), NOSPLIT, $0-20 MOVD a+0(FP), R0 MOVD n+8(FP), R1 diff --git a/internal/kernels/sum_numpy_test.go b/internal/kernels/sum_numpy_test.go new file mode 100644 index 0000000..0585b4c --- /dev/null +++ b/internal/kernels/sum_numpy_test.go @@ -0,0 +1,37 @@ +package kernels + +import ( + "math" + "testing" +) + +// TestSumMatchesNumPy holds Sum's value (SumP, sumSIMD) to np.sum bit for bit. +// The golden values are NumPy 2.5.4's (np.sum of a[i] = -1 + 2*frac(0.618033988749 +// 8949*(i+1)), computed with np.arange, not by this package): +// +// uv run --no-project --python 3.12.12 --with numpy==2.5.4 python -I -c \ +// 'import numpy as np; i=np.arange(1000)+1.0; print(float(np.sum(-1+2*((i*0.6180339887498949)%1.0))).hex())' +// +// Since K1 of the go-scikit performance campaign the sum is floats.SumPairwise, +// NumPy's order; the lane-parallel sums it replaced differed in the last bits. +func TestSumMatchesNumPy(t *testing.T) { + for _, c := range []struct { + n int + want float64 + }{ + {1000, 0x1.748cc9f0cbc00p-6}, + {100000, 0x1.cbaf96d9ef5e0p-1}, + } { + a := make([]float64, c.n) + for i := range a { + ph := float64(float64(i+1) * 0.6180339887498949) + a[i] = -1 + float64(2*(ph-math.Floor(ph))) + } + if got := SumP(a); math.Float64bits(got) != math.Float64bits(c.want) { + t.Errorf("SumP n=%d: %v (%#x), NumPy %v (%#x)", c.n, got, math.Float64bits(got), c.want, math.Float64bits(c.want)) + } + if got := sumSIMD(a); math.Float64bits(got) != math.Float64bits(c.want) { + t.Errorf("sumSIMD n=%d: %v, NumPy %v", c.n, got, c.want) + } + } +} diff --git a/internal/kernels/via64.go b/internal/kernels/via64.go deleted file mode 100644 index 615ca01..0000000 --- a/internal/kernels/via64.go +++ /dev/null @@ -1,37 +0,0 @@ -package kernels - -// Correctly rounded float32 Exp and Log. -// -// float32 Exp and Log went through math.Exp and math.Log per element, 6-14x -// slower than NumPy's float32 SIMD on one Zen 3 core. Both now compute in -// float64 and round once, and both are proven correctly rounded over every -// float32 input (via64_exhaustive_test.go, expcr_test.go, logcr_test.go): -// -// - Exp runs expCR (expcr.go), Arm's correctly rounded expf algorithm; -// - Log runs logCR (logcr.go), CORE-MATH's correctly rounded logf; -// -// each on an AVX2 kernel where there is one (expcr_amd64.go, logcr_amd64.go), -// in Go elsewhere. Nothing is allocated. - -const via64Chunk = 256 - -// Exp32 writes the correctly rounded exp(src[i]) into dst[i]. -func Exp32(dst, src []float32) { exp32Fast(dst, src) } - -// Log32 writes the correctly rounded log(src[i]) into dst[i]. -func Log32(dst, src []float32) { log32Fast(dst, src) } - -// Exp32P is Exp32, split across cores above ParThreshold. -func Exp32P(dst, src []float32) { split32(exp32Fast, dst, src) } - -// Log32P is Log32, split across cores above ParThreshold. -func Log32P(dst, src []float32) { split32(log32Fast, dst, src) } - -func split32(f func(dst, src []float32), dst, src []float32) { - n := len(dst) - if n < ParThreshold { - f(dst, src) - return - } - parallelFor(n, numWorkers(n), func(lo, hi int) { f(dst[lo:hi], src[lo:hi]) }) -} diff --git a/internal/kernels/via64_exhaustive_test.go b/internal/kernels/via64_exhaustive_test.go deleted file mode 100644 index 7731014..0000000 --- a/internal/kernels/via64_exhaustive_test.go +++ /dev/null @@ -1,267 +0,0 @@ -package kernels - -import ( - "math" - "math/big" - "os" - "runtime" - "sync" - "sync/atomic" - "testing" -) - -// Exp32 and Log32 are correctly rounded: proven over every float32. -// -// A float32 function has 2^32 inputs, so correct rounding can be checked -// exhaustively instead of sampled. For each input the float64 intermediate y -// (Exp or Log of the widened input, the value Exp32 and Log32 round) has an -// error of at most 0.504 float64 ULP (TestExpULP, TestLogULP). If y is more -// than 2 float64 ULP away from both float32 rounding midpoints around -// float32(y), the exact value rounds to the same float32: certified. The few -// inputs closer than that are decided by the 300-bit references refExp and -// refLog, rounded by big.Float.Float32 (round to nearest, ties to even). -// Results far beyond float32's range are certified by bounds: |y| > 2^129 -// rounds to ±Inf, |y| < 2^-152 to ±0, and so does the exact value. -// -// The full run (2^32 inputs per function) takes seconds on all cores and is -// set off by NDARRAY_EXHAUSTIVE=1; without it every 65537th bit pattern is -// checked, about 65 000 inputs, so CI stays fast. - -// certified32 rounds the float64 intermediate y to float32 and reports -// whether that rounding is provably the correct rounding of the exact value. -func certified32(y float64) (float32, bool) { return certified32Margin(y, 2) } - -// certified32Margin is certified32 for an intermediate whose error may reach -// margin/2 float64 ULP (the margin is kept at twice the measured worst). -func certified32Margin(y float64, margin float64) (float32, bool) { - f := float32(y) - ay := math.Abs(y) - if math.IsNaN(y) || math.IsInf(y, 0) || ay > 0x1p129 || ay < 0x1p-152 { - return f, true - } - eps := margin * (math.Nextafter(ay, math.Inf(1)) - ay) - if math.IsInf(float64(f), 0) { - // y rounded to ±Inf although |y| <= 2^129: the only boundary is - // half an ULP past MaxFloat32, where ties-to-even goes to Inf. - // (TestCertified32 plants exactly this value.) - return f, ay-(math.MaxFloat32+0x1p103) > eps - } - ff := float64(f) - near := func(neighbour float32) bool { - var m float64 - if math.IsInf(float64(neighbour), 0) { // f is ±MaxFloat32: the boundary is half an ULP past it - m = ff + math.Copysign(0x1p103, ff) - } else { - m = (ff + float64(neighbour)) / 2 // exact: adjacent float32s have 24-bit mantissas - } - return math.Abs(y-m) <= eps - } - if near(math.Nextafter32(f, float32(math.Inf(-1)))) || near(math.Nextafter32(f, float32(math.Inf(1)))) { - return f, false - } - return f, true -} - -type exhaustiveCase struct { - name string - k32 func(dst, src []float32) - k64 func(dst, src []float64) - ref func(float64) *big.Float - // special returns the expected result for inputs the reference does not - // take (NaN, ±Inf, and for log 0 and negatives), and ok=false otherwise. - special func(x float32) (want float32, isNaN, ok bool) - // margin is the certifier's distance to a midpoint, in float64 ULP - // (0 means 2, for the 0.504-ULP kernels). - margin float64 -} - -func exp32Special(x float32) (float32, bool, bool) { - switch { - case x != x: - return 0, true, true - case math.IsInf(float64(x), 1): - return float32(math.Inf(1)), false, true - case math.IsInf(float64(x), -1): - return 0, false, true - } - return 0, false, false -} - -func log32Special(x float32) (float32, bool, bool) { - switch { - case x != x || x < 0: - return 0, true, true - case x == 0: - return float32(math.Inf(-1)), false, true - case math.IsInf(float64(x), 1): - return float32(math.Inf(1)), false, true - } - return 0, false, false -} - -func runExhaustive(t *testing.T, c exhaustiveCase, stride uint64) { - const batch = 4096 - workers := runtime.GOMAXPROCS(0) - var next atomic.Uint64 - var bad, hard, checked atomic.Int64 - var mu sync.Mutex - var firstBad []string - var wg sync.WaitGroup - for range workers { - wg.Add(1) - go func() { - defer wg.Done() - xs := make([]float32, 0, batch) - got := make([]float32, batch) - w64 := make([]float64, batch) - y64 := make([]float64, batch) - for { - start := next.Add(batch*stride) - batch*stride - if start >= 1<<32 { - return - } - xs = xs[:0] - for b := start; b < start+batch*stride && b < 1<<32; b += stride { - xs = append(xs, math.Float32frombits(uint32(b))) - } - n := len(xs) - c.k32(got[:n], xs) - for i, x := range xs { - w64[i] = float64(x) - } - c.k64(y64[:n], w64[:n]) - for i, x := range xs { - checked.Add(1) - var want float32 - if w, isNaN, ok := c.special(x); ok { - if isNaN { - if got[i] == got[i] { - bad.Add(1) - mu.Lock() - if len(firstBad) < 8 { - firstBad = append(firstBad, c.name+": NaN wanted") - } - mu.Unlock() - } - continue - } - want = w - } else { - margin := c.margin - if margin == 0 { - margin = 2 - } - f, sure := certified32Margin(y64[i], margin) - if sure { - want = f - } else { - hard.Add(1) - want, _ = c.ref(float64(x)).Float32() - if os.Getenv("NDARRAY_DUMP_HARD") != "" { - mu.Lock() - t.Logf("hard %s 0x%s -> 0x%s", c.name, hex32(math.Float32bits(x)), hex32(math.Float32bits(want))) - mu.Unlock() - } - } - } - if math.Float32bits(got[i]) != math.Float32bits(want) { - bad.Add(1) - mu.Lock() - if len(firstBad) < 8 { - firstBad = append(firstBad, c.name+"("+formatF32(x)+") = "+formatF32(got[i])+", correctly rounded "+formatF32(want)) - } - mu.Unlock() - } - } - } - }() - } - wg.Wait() - t.Logf("%s: %d inputs checked, %d decided by the 300-bit reference, %d not correctly rounded", - c.name, checked.Load(), hard.Load(), bad.Load()) - for _, s := range firstBad { - t.Error(s) - } -} - -func formatF32(x float32) string { - return big.NewFloat(float64(x)).Text('g', 9) + " (0x" + hex32(math.Float32bits(x)) + ")" -} - -func hex32(u uint32) string { - const digits = "0123456789abcdef" - var b [8]byte - for i := 7; i >= 0; i-- { - b[i] = digits[u&15] - u >>= 4 - } - return string(b[:]) -} - -func exhaustiveStride() uint64 { - if os.Getenv("NDARRAY_EXHAUSTIVE") != "" { - return 1 - } - return 65537 -} - -func TestExp32CorrectlyRounded(t *testing.T) { - // Exp32 rounds expCR (1.006 float64 ULP): the certifier's margin is 4. - runExhaustive(t, exhaustiveCase{"exp", Exp32, expCRSlice, refExp, exp32Special, 4}, exhaustiveStride()) -} - -func TestLog32CorrectlyRounded(t *testing.T) { - runExhaustive(t, exhaustiveCase{"log", Log32, Log, refLog, log32Special, 2}, exhaustiveStride()) -} - -// TestLog32Exceptions runs Log32 on every input of the exception table (the -// sampled run rarely meets one): each result must be the table's, and the -// table's must be the 300-bit reference rounded to float32. An input outside -// the table keeps the computed value. -func TestLog32Exceptions(t *testing.T) { - src := make([]float32, len(log32Exceptions)) - for i, e := range log32Exceptions { - src[i] = math.Float32frombits(e[0]) - } - got := make([]float32, len(src)) - Log32(got, src) - for i, e := range log32Exceptions { - want, _ := refLog(float64(src[i])).Float32() - if math.Float32bits(want) != e[1] { - t.Errorf("table entry 0x%s: 0x%s, the reference rounds to 0x%s", hex32(e[0]), hex32(e[1]), hex32(math.Float32bits(want))) - } - if math.Float32bits(got[i]) != e[1] { - t.Errorf("Log32(0x%s) = 0x%s, want the table's 0x%s", hex32(e[0]), hex32(math.Float32bits(got[i])), hex32(e[1])) - } - } - if f := log32Exception(2, 0.5); f != 0.5 { - t.Errorf("an input outside the table changed: %v", f) - } -} - -// TestCertified32 checks the certifier itself on planted values: a float64 -// exactly on a float32 midpoint, and one ULP from it, must be sent to the -// reference; a value in the middle of a float32 interval is certified. -func TestCertified32(t *testing.T) { - one, next := float64(float32(1)), float64(math.Nextafter32(1, 2)) - mid := (one + next) / 2 - if _, sure := certified32(mid); sure { - t.Error("a float32 midpoint was certified") - } - if _, sure := certified32(math.Nextafter(mid, 2)); sure { - t.Error("one float64 ULP above a midpoint was certified") - } - if f, sure := certified32(one + (next-one)/8); !sure || f != 1 { - t.Errorf("1 + ulp/8: %v, %v", f, sure) - } - if f, sure := certified32(0x1p130); !sure || !math.IsInf(float64(f), 1) { - t.Errorf("2^130: %v, %v", f, sure) - } - if f, sure := certified32(0x1p-160); !sure || f != 0 { - t.Errorf("2^-160: %v, %v", f, sure) - } - maxMid := float64(math.MaxFloat32) + 0x1p103 - if _, sure := certified32(maxMid); sure { - t.Error("the overflow boundary was certified") - } -} diff --git a/internal/kernels/via64_oracle_test.go b/internal/kernels/via64_oracle_test.go deleted file mode 100644 index 659a324..0000000 --- a/internal/kernels/via64_oracle_test.go +++ /dev/null @@ -1,68 +0,0 @@ -package kernels - -import "math" - -// log32 is the float32 log of v0.32.0 and v0.33.0, kept as an independent -// oracle for logCR (TestLogCRMatchesLog32): it widens to this package's -// float64 Log (0.504 ULP), rounds once, and takes the 24 inputs whose float64 -// result lies next to a float32 rounding midpoint from a table. It shares no -// code or constant with logCR. -func log32(dst, src []float32) { - var buf [via64Chunk]float64 - src = src[:len(dst)] - for i := 0; i < len(dst); i += via64Chunk { - b := buf[:min(via64Chunk, len(dst)-i)] - s := src[i : i+len(b)] - for j := range b { - b[j] = float64(s[j]) - } - Log(b, b) - d := dst[i : i+len(b)] - for j, v := range b { - f := float32(v) - // Within 8 float64 ULP of a float32 rounding midpoint (the 29 - // bits float32 drops are 0x10000000 +- 8): the double rounding - // may be wrong, so take the table's answer. - if math.Float64bits(v)&(1<<29-1)-(1<<28-8) < 16 { - f = log32Exception(s[j], f) - } - d[j] = f - } - } -} - -// log32Exceptions are the float32 inputs whose float64 log lies within 2 -// float64 ULP of a float32 rounding midpoint, sorted, each with its correctly -// rounded log decided by a 300-bit reference (NDARRAY_EXHAUSTIVE=1 -// NDARRAY_DUMP_HARD=1 on the v0.32.0 TestLog32CorrectlyRounded): 24 inputs, 5 -// of which rounded to the wrong float32 without this table. They stay the -// hardest known inputs for any float32 log (TestLog32Exceptions). -var log32Exceptions = [...][2]uint32{ - {0x0dc8bba4, 0xc289bac4}, {0x111c87f8, 0xc28085df}, {0x1a8446cb, 0xc24d0a45}, - {0x1f116ab8, 0xc233b53e}, {0x28e3fa26, 0xc1fa75e0}, {0x29e6126b, 0xc1ef4c02}, - {0x29fd22f8, 0xc1ee8859}, {0x2c4c24b7, 0xc1d48710}, {0x2e492984, 0xc1be76ee}, - {0x38dcbe38, 0xc1128ba4}, {0x3bf86ef0, 0xc09c399e}, {0x3c413d3a, 0xc08e158f}, - {0x41178feb, 0x400fe5e7}, {0x464d5b2b, 0x4117bcf1}, {0x4665a9a6, 0x41198725}, - {0x4c5d65a5, 0x418f034b}, {0x4d604ebe, 0x419a352c}, {0x4e85f412, 0x41a6b811}, - {0x5ee8984e, 0x422e4a21}, {0x65d890d3, 0x4254d1f9}, {0x66a8c860, 0x42595e46}, - {0x66abbd63, 0x4259700e}, {0x6f31a8ec, 0x42845a89}, {0x79e7ec37, 0x42a1ffb7}, -} - -// log32Exception returns the table's correctly rounded log of x if x is in -// log32Exceptions, and f otherwise. -func log32Exception(x, f float32) float32 { - bits := math.Float32bits(x) - lo, hi := 0, len(log32Exceptions) - for lo < hi { - m := (lo + hi) / 2 - switch e := log32Exceptions[m]; { - case e[0] == bits: - return math.Float32frombits(e[1]) - case e[0] < bits: - lo = m + 1 - default: - hi = m - } - } - return f -} diff --git a/internal/kernels/via64_test.go b/internal/kernels/via64_test.go deleted file mode 100644 index da20764..0000000 --- a/internal/kernels/via64_test.go +++ /dev/null @@ -1,146 +0,0 @@ -package kernels - -import ( - "math" - "math/rand" - "testing" -) - -// ulps32 is the distance between two float32s in units in the last place -// (0 for equal values, including two NaNs). -func ulps32(a, b float32) int64 { - if a == b || (a != a && b != b) { - return 0 - } - ia, ib := int64(math.Float32bits(a)), int64(math.Float32bits(b)) - if ia < 0 { - ia = math.MinInt32 - ia - } - if ib < 0 { - ib = math.MinInt32 - ib - } - if d := ia - ib; d < 0 { - return -d - } else { - return d - } -} - -// TestVia64MatchesMath holds Exp32 and Log32 to float32(math.Exp/Log(x)): -// both round a float64 result once, so they may differ by one float32 ULP -// where math's float64 result (0.88 ULP) and ours (0.504) fall on opposite -// sides of a rounding midpoint, and must agree on at least 99.9% of inputs. -// Every length 0..600 crosses the 256-element chunk boundaries. -func TestVia64MatchesMath(t *testing.T) { - r := rand.New(rand.NewSource(7)) - for _, f := range []struct { - name string - k func(dst, src []float32) - ref func(float64) float64 - gen func() float32 - }{ - {"Exp", Exp32, math.Exp, func() float32 { return float32(r.Float64()*200 - 110) }}, - {"Log", Log32, math.Log, func() float32 { - return float32(math.Exp(r.Float64()*170 - 85)) - }}, - } { - const n = 100000 - src, dst := make([]float32, n), make([]float32, n) - for i := range src { - src[i] = f.gen() - } - f.k(dst, src) - var off int - for i, x := range src { - want := float32(f.ref(float64(x))) - switch d := ulps32(dst[i], want); { - case d > 1: - t.Fatalf("%s(%v) = %v, math gives %v (%d ULP)", f.name, x, dst[i], want, d) - case d == 1: - off++ - } - } - if off > n/1000 { - t.Errorf("%s: %d of %d differ from math by one ULP", f.name, off, n) - } - for m := 0; m <= 600; m++ { - d := make([]float32, m) - f.k(d, src[:m]) - for i := range d { - if d[i] != dst[i] { - t.Fatalf("%s length %d [%d]: %v != %v", f.name, m, i, d[i], dst[i]) - } - } - } - } -} - -// TestVia64Specials checks the edges that come from rounding the float64 -// result: overflow to +Inf, underflow through the subnormals to 0, and the -// IEEE special values. -func TestVia64Specials(t *testing.T) { - inf, nan := float32(math.Inf(1)), float32(math.NaN()) - negZero := float32(math.Copysign(0, -1)) - for _, c := range []struct { - name string - k func(dst, src []float32) - in, out float32 - }{ - {"Exp", Exp32, 0, 1}, - {"Exp", Exp32, negZero, 1}, - {"Exp", Exp32, 89, inf}, // finite in float64, overflows float32 - {"Exp", Exp32, -inf, 0}, - {"Exp", Exp32, inf, inf}, - {"Exp", Exp32, -104, 0}, // below the smallest float32 subnormal - {"Exp", Exp32, -100, float32(math.Exp(-100))}, - {"Log", Log32, 1, 0}, - {"Log", Log32, 0, -inf}, - {"Log", Log32, negZero, -inf}, - {"Log", Log32, inf, inf}, - {"Log", Log32, math.SmallestNonzeroFloat32, float32(math.Log(math.SmallestNonzeroFloat32))}, - } { - d := []float32{0} - c.k(d, []float32{c.in}) - if d[0] != c.out || math.Signbit(float64(d[0])) != math.Signbit(float64(c.out)) { - t.Errorf("%s(%v) = %v, want %v", c.name, c.in, d[0], c.out) - } - } - for _, in := range []float32{nan, -1, -inf} { - d := []float32{0} - Log32(d, []float32{in}) - if d[0] == d[0] { - t.Errorf("Log(%v) = %v, want NaN", in, d[0]) - } - } - d := []float32{0} - Exp32(d, []float32{nan}) - if d[0] == d[0] { - t.Errorf("Exp(NaN) = %v", d[0]) - } -} - -// TestVia64InPlaceAndParallel: dst may be src, and the parallel entries give -// the serial bits. -func TestVia64InPlaceAndParallel(t *testing.T) { - src := make([]float32, 3*ParThreshold+17) - for i := range src { - src[i] = float32(i%97) * 0.25 - } - want := make([]float32, len(src)) - Exp32(want, src) - got := append([]float32(nil), src...) - Exp32(got, got) - par := make([]float32, len(src)) - withMaxProcs(4, func() { Exp32P(par, src) }) - for i := range want { - if got[i] != want[i] || par[i] != want[i] { - t.Fatalf("[%d]: in place %v, parallel %v, serial %v", i, got[i], par[i], want[i]) - } - } - Log32P(par, src) // the serial path below ParThreshold is covered by the tests above - for i, x := range src { - if par[i] != float32(math.Log(float64(x))) && ulps32(par[i], float32(math.Log(float64(x)))) > 1 { - t.Fatalf("Log32P [%d]: %v", i, par[i]) - } - } -} From 5a4412ca8aa77fd5315c97da30bff454fac3e1bd Mon Sep 17 00:00:00 2001 From: tannevaled Date: Sat, 7 Nov 2026 22:32:30 +0100 Subject: [PATCH 2/2] kernels: test the parallel exp/log forms and Log10 (100% coverage gate) Co-Authored-By: Claude Sonnet 5.5 --- internal/kernels/explog_test.go | 52 +++++++++++++++++++++++++++++++++ 1 file changed, 52 insertions(+) create mode 100644 internal/kernels/explog_test.go diff --git a/internal/kernels/explog_test.go b/internal/kernels/explog_test.go new file mode 100644 index 0000000..fd7a987 --- /dev/null +++ b/internal/kernels/explog_test.go @@ -0,0 +1,52 @@ +package kernels + +import ( + "math" + "testing" +) + +// The exp/log kernels are go-simd/floats' (proofs and kernel tests live there); +// what stays here is the parallel split and Log10: every P form gives the serial +// bits, below and above ParThreshold. +func TestExpLogParallelMatchesSerial(t *testing.T) { + src := randVec(5000, 3) + for i := range src { + src[i] = math.Abs(src[i]) * 40 + } + src32 := make([]float32, len(src)) + for i, v := range src { + src32[i] = float32(v) + } + run := func(name string, serial, par func(dst, src []float64)) { + want, got := make([]float64, len(src)), make([]float64, len(src)) + serial(want, src) + withThresholds(64, 1<<14, func() { par(got, src) }) + eqSlice(t, got, want) + par(got[:10], src[:10]) // below the threshold + eqSlice(t, got[:10], want[:10]) + _ = name + } + run("exp", Exp, ExpP) + run("log", Log, LogP) + run("log10", Log10, Log10P) + for _, c := range []struct { + name string + serial, par func(dst, src []float32) + }{{"exp32", Exp32, Exp32P}, {"log32", Log32, Log32P}} { + want, got := make([]float32, len(src32)), make([]float32, len(src32)) + c.serial(want, src32) + withThresholds(64, 1<<14, func() { c.par(got, src32) }) + for i := range want { + if got[i] != want[i] { + t.Fatalf("%sP [%d]: %v, serial %v", c.name, i, got[i], want[i]) + } + } + c.par(got[:10], src32[:10]) + } + // Log10 is log(x)*(1/ln10). + d := make([]float64, 1) + Log10(d, []float64{1000}) + if math.Abs(d[0]-3) > 1e-15 { + t.Errorf("Log10(1000) = %v", d[0]) + } +}