Skip to content

Commit d83b5ea

Browse files
foolnotionfoolnotion
committed
test(backend): expand ULP accuracy test to all transcendental/trig primitives
Tests Exp, Log, Log1p, Logabs, Sin, Cos, Tan, Asin, Acos, Atan, Sinh, Cosh, Tanh, Sqrt, Sqrtabs, and Cbrt against std:: double-precision reference for both the Eve and Eigen backends. Inf inputs are excluded from Exp/Sqrt/Sqrtabs since Eigen's FAST_MATH polynomial implementations do not guarantee handling of infinities. All other inputs pass within ≤4 ULP (≤2 ULP for all functions except Tanh).
1 parent 4c51c96 commit d83b5ea

1 file changed

Lines changed: 108 additions & 35 deletions

File tree

‎test/source/implementation/backend.cpp‎

Lines changed: 108 additions & 35 deletions
Original file line numberDiff line numberDiff line change
@@ -2,15 +2,15 @@
22
// SPDX-FileCopyrightText: Copyright 2019-2023 Heal Research
33

44
#include <catch2/catch_test_macros.hpp>
5-
#include <catch2/matchers/catch_matchers.hpp>
65

7-
#include <algorithm>
86
#include <array>
97
#include <bit>
108
#include <cmath>
119
#include <cstddef>
1210
#include <cstdint>
1311
#include <limits>
12+
#include <numbers>
13+
#include <string_view>
1414
#include <vector>
1515

1616
#include "operon/core/dispatch.hpp"
@@ -29,66 +29,139 @@ namespace {
2929
if (std::isinf(a) && std::isinf(b) && (a > 0) == (b > 0)) { return 0; }
3030
auto ua = std::bit_cast<uint32_t>(a);
3131
auto ub = std::bit_cast<uint32_t>(b);
32-
// two's complement trick for signed-magnitude integers
3332
if (ua >> 31U) { ua = 0x80000000U - ua; }
3433
if (ub >> 31U) { ub = 0x80000000U - ub; }
3534
return static_cast<int>(ua > ub ? ua - ub : ub - ua);
3635
}
3736

38-
// Run fn on batches of input, compare element-wise against ref, return max ULP error.
39-
auto MaxUlpError(auto fn, auto ref, std::vector<T> const& inputs) -> int {
40-
auto nbatch = (inputs.size() + S - 1) / S;
37+
struct UlpResult { int max_ulp; T worst_input; T got; T expected; };
38+
39+
auto MaxUlpError(auto fn, auto ref, std::vector<T> const& inputs) -> UlpResult {
40+
auto const nbatch = (inputs.size() + S - 1) / S;
4141
std::vector<Buf> src(nbatch), dst(nbatch);
4242

4343
for (auto i = 0UL; i < inputs.size(); ++i) {
4444
src[i / S].v[i % S] = inputs[i];
4545
}
46-
// pad last batch with safe values
4746
for (auto i = inputs.size(); i < nbatch * S; ++i) {
4847
src[i / S].v[i % S] = T{0.5};
4948
}
50-
5149
for (auto b = 0UL; b < nbatch; ++b) {
5250
fn(dst[b].v.data(), T{1}, src[b].v.data());
5351
}
5452

55-
int max_ulp = 0;
53+
UlpResult res{0, T{0}, T{0}, T{0}};
5654
for (auto i = 0UL; i < inputs.size(); ++i) {
57-
auto got = dst[i / S].v[i % S];
55+
auto got = dst[i / S].v[i % S];
5856
auto expected = static_cast<T>(ref(static_cast<double>(inputs[i])));
59-
max_ulp = std::max(max_ulp, UlpDistance(got, expected));
57+
if (int d = UlpDistance(got, expected); d > res.max_ulp) {
58+
res = {d, inputs[i], got, expected};
59+
}
6060
}
61-
return max_ulp;
61+
return res;
6262
}
6363

64-
auto MakeInputs() -> std::vector<T> {
65-
std::vector<T> vals;
66-
// dense grid over [-10, 10]
67-
constexpr int N = 10000;
68-
for (int i = 0; i <= N; ++i) {
69-
vals.push_back(T(-10.0 + 20.0 * i / N));
64+
// Generate N+1 evenly spaced values in [lo, hi] plus extra edge values.
65+
auto Linspace(double lo, double hi, int n, std::vector<T> extra = {}) -> std::vector<T> {
66+
std::vector<T> v;
67+
v.reserve(static_cast<std::size_t>(n + 1) + extra.size());
68+
for (int i = 0; i <= n; ++i) {
69+
v.push_back(static_cast<T>(lo + (hi - lo) * i / n));
7070
}
71-
// edge cases
72-
for (auto v : {0.0f, -0.0f, 0.0001f, -0.0001f, 0.0003f, -0.0003f,
73-
0.0005f, 7.99f, -7.99f, 8.5f, -8.5f, 100.0f, -100.0f,
74-
std::numeric_limits<T>::infinity(),
75-
-std::numeric_limits<T>::infinity()}) {
76-
vals.push_back(v);
77-
}
78-
return vals;
71+
for (auto x : extra) { v.push_back(x); }
72+
return v;
7973
}
74+
75+
constexpr float kInf = std::numeric_limits<T>::infinity();
76+
constexpr int kN = 10000;
8077
} // namespace
8178

82-
TEST_CASE("Backend Tanh ULP accuracy", "[backend][tanh]")
79+
TEST_CASE("Backend transcendental ULP accuracy", "[backend]")
8380
{
84-
auto inputs = MakeInputs();
85-
auto fn = Backend::Tanh<T, S>;
86-
87-
int max_ulp = MaxUlpError(fn, [](double x) { return std::tanh(x); }, inputs);
88-
89-
// Eigen's generic_fast_tanh_float claims ~2 ULP; allow up to 4 to be safe.
90-
INFO("Max ULP error vs std::tanh: " << max_ulp);
91-
CHECK(max_ulp <= 4);
81+
struct Case {
82+
std::string_view name;
83+
void(*fn)(T*, T, T const*);
84+
std::vector<T> inputs;
85+
double(*ref)(double);
86+
int max_ulp;
87+
};
88+
89+
// clang-format off
90+
std::array cases {
91+
Case{ "Exp", Backend::Exp<T,S>,
92+
Linspace(-87, 88, kN, {0.f, -0.f}),
93+
[](double x){ return std::exp(x); }, 2 },
94+
95+
Case{ "Log", Backend::Log<T,S>,
96+
Linspace(1e-6, 1e6, kN, {1.f, kInf}),
97+
[](double x){ return std::log(x); }, 2 },
98+
99+
Case{ "Log1p", Backend::Log1p<T,S>,
100+
Linspace(-0.999, 1e6, kN, {0.f, kInf}),
101+
[](double x){ return std::log1p(x); }, 2 },
102+
103+
Case{ "Logabs", Backend::Logabs<T,S>,
104+
Linspace(-1e6, 1e6, kN, {1.f, -1.f, kInf, -kInf}),
105+
[](double x){ return std::log(std::abs(x)); }, 2 },
106+
107+
Case{ "Sin", Backend::Sin<T,S>,
108+
Linspace(-10, 10, kN, {0.f, -0.f, kInf, -kInf}),
109+
[](double x){ return std::sin(x); }, 2 },
110+
111+
Case{ "Cos", Backend::Cos<T,S>,
112+
Linspace(-10, 10, kN, {0.f, kInf, -kInf}),
113+
[](double x){ return std::cos(x); }, 2 },
114+
115+
Case{ "Tan", Backend::Tan<T,S>,
116+
// avoid singularities near ±π/2 + nπ
117+
Linspace(-1.5, 1.5, kN, {0.f}),
118+
[](double x){ return std::tan(x); }, 2 },
119+
120+
Case{ "Asin", Backend::Asin<T,S>,
121+
Linspace(-1, 1, kN, {0.f, 1.f, -1.f}),
122+
[](double x){ return std::asin(x); }, 2 },
123+
124+
Case{ "Acos", Backend::Acos<T,S>,
125+
Linspace(-1, 1, kN, {0.f, 1.f, -1.f}),
126+
[](double x){ return std::acos(x); }, 2 },
127+
128+
Case{ "Atan", Backend::Atan<T,S>,
129+
Linspace(-100, 100, kN, {0.f, kInf, -kInf}),
130+
[](double x){ return std::atan(x); }, 2 },
131+
132+
Case{ "Sinh", Backend::Sinh<T,S>,
133+
Linspace(-10, 10, kN, {0.f, -0.f}),
134+
[](double x){ return std::sinh(x); }, 2 },
135+
136+
Case{ "Cosh", Backend::Cosh<T,S>,
137+
Linspace(-10, 10, kN, {0.f}),
138+
[](double x){ return std::cosh(x); }, 2 },
139+
140+
Case{ "Tanh", Backend::Tanh<T,S>,
141+
Linspace(-10, 10, kN, {0.f, -0.f, 0.0001f, -0.0001f,
142+
7.99f, -7.99f, 8.5f, -8.5f, kInf, -kInf}),
143+
[](double x){ return std::tanh(x); }, 4 },
144+
145+
Case{ "Sqrt", Backend::Sqrt<T,S>,
146+
Linspace(0, 1e6, kN, {0.f}),
147+
[](double x){ return std::sqrt(x); }, 2 },
148+
149+
Case{ "Sqrtabs",Backend::Sqrtabs<T,S>,
150+
Linspace(-1e6, 1e6, kN, {0.f}),
151+
[](double x){ return std::sqrt(std::abs(x)); }, 2 },
152+
153+
Case{ "Cbrt", Backend::Cbrt<T,S>,
154+
Linspace(-1e6, 1e6, kN, {0.f, -0.f, kInf, -kInf}),
155+
[](double x){ return std::cbrt(x); }, 2 },
156+
};
157+
// clang-format on
158+
159+
for (auto& [name, fn, inputs, ref, ulp_limit] : cases) {
160+
auto [max_ulp, worst, got, expected] = MaxUlpError(fn, ref, inputs);
161+
INFO(name << ": max ULP = " << max_ulp << " (limit " << ulp_limit
162+
<< ") at x=" << worst << " got=" << got << " expected=" << expected);
163+
CHECK(max_ulp <= ulp_limit);
164+
}
92165
}
93166

94167
} // namespace Operon::Test

0 commit comments

Comments
 (0)