Skip to content

Instantly share code, notes, and snippets.

@zeux

zeux/lerpfuzz.cpp

Last active December 11, 2024 07:13
Show Gist options
  • Save zeux/ea70d4574a50cc819e5f02b6d7ebb1f2 to your computer and use it in GitHub Desktop.
Save zeux/ea70d4574a50cc819e5f02b6d7ebb1f2 to your computer and use it in GitHub Desktop.
Download ZIP
Fuzzer to investigate properties of lerp implementations
Raw
lerpfuzz.cpp
// clang++ -fsanitize=fuzzer -ffp-contract=off -O1 ./lerpfuzz.cpp -o lerpfuzz && ./lerpfuzz
#include <assert.h>
#include <math.h>
#include <stdint.h>
#include <stdio.h>
#include <string.h>
using real_t = float;
#define LERP 2
#define TEST_EXACT 1
#define TEST_MONOT 1
#define TEST_DETER 0
#define TEST_BOUND 1
#define TEST_CONSI 1
// a. exactness: lerp(0) = a, lerp(1) = b
// b. monotonicity: lerp(t1) >= lerp(t2) for t1 >= t2
// c. determinacy: lerp(a, b, t) is NaN only if t is infinite or a/b/t are NaN
// d. boundedness: lerp(t) in [a, b]
// e. consistency: lerp(x, x, t) == x for all t
// note: if an implementation is exact but not consistent, it can't be monotonic: lerp(x, x, t) will be either <lerp(x, x, 1) or >lerp(x, x, 0) for some value of t
// additional references:
// https://www.open-std.org/jtc1/sc22/wg21/docs/papers/2019/p0811r3.html
// https://math.stackexchange.com/questions/907327/accurate-floating-point-linear-interpolation
real_t lerp(real_t a, real_t b, real_t t)
{
#if LERP == 0
// exactness 0, monotonicity 1, determinacy 0, boundedness 0, consistency 1
return a + (b - a) * t;
#elif LERP == 1
// exactness 1, monotonicity 0, determinacy 1, boundedness 0, consistency 0
return a * (1 - t) + b * t;
#elif LERP == 2
// exactness 1, monotonicity 1, determinacy 0, boundedness 1, consistency 1
return t == 1 ? b : a + (b - a) * t;
#elif LERP == 3
// exactness 1, monotonicity 1, determinacy 0, boundedness 1, consistency 1
return t < 0.5f ? a + (b - a) * t : b + (b - a) * (t - 1);
#elif LERP == 4
// exactness 1, monotonicity 0, determinacy 1, boundedness 0, consistency 1
return a == b ? a : a * (1 - t) + b * t;
#endif
}
float ulp(float v, int o)
{
unsigned int u;
memcpy(&u, &v, 4);
if (v == 0) u = 0;
u += o;
memcpy(&v, &u, 4);
return v;
}
double ulp(double v, int o)
{
unsigned long long u;
memcpy(&u, &v, 8);
if (v == 0) u = 0;
u += o;
memcpy(&v, &u, 8);
return v;
}
extern "C" int LLVMFuzzerTestOneInput(const uint8_t* buffer, size_t size)
{
if (size < sizeof(real_t) * 3)
return 0;
real_t a, b, t;
memcpy(&a, buffer + 0 * sizeof(real_t), sizeof(real_t));
memcpy(&b, buffer + 1 * sizeof(real_t), sizeof(real_t));
memcpy(&t, buffer + 2 * sizeof(real_t), sizeof(real_t));
// preconditions
if (!isfinite(a))
return 0;
if (!isfinite(b))
return 0;
if (!isfinite(t))
return 0;
if (t < 0 || t > 1)
return 0;
if (a > b)
return 0;
real_t v = lerp(a, b, t);
#if TEST_EXACT
// a. exactness: lerp(0) = a, lerp(1) = b
if (isfinite(b - a))
{
assert(lerp(a, b, 0) == a);
assert(lerp(a, b, 1) == b);
}
#endif
#if TEST_MONOT
// b. monotonicity: lerp(t1) >= lerp(t2) for t1 >= t2
if (isfinite(b - a))
{
if (t > 0)
{
real_t pt = ulp(t, -1);
assert(pt < t);
assert(lerp(a, b, pt) <= v);
}
if (t < 1)
{
real_t nt = ulp(t, +1);
assert(nt > t);
assert(lerp(a, b, nt) >= v);
}
{
// monotonicity around 1
real_t p1 = ulp(real_t(1), -1);
real_t n1 = ulp(real_t(1), +1);
assert(p1 < 1);
assert(n1 > 1);
assert(lerp(a, b, p1) <= lerp(a, b, 1));
assert(lerp(a, b, n1) >= lerp(a, b, 1));
}
}
#endif
#if TEST_DETER
// c. determinacy: lerp(a, b, t) is NaN only if t is infinite or a/b/t are NaN
{
assert(isfinite(v));
real_t p1 = ulp(real_t(1), -1);
real_t l1 = lerp(a, b, p1);
assert(isfinite(l1));
}
#endif
#if TEST_BOUND
// d. boundedness: lerp(t) in [a, b]
if (isfinite(b - a))
{
assert(v >= a && v <= b);
}
#endif
#if TEST_CONSI
// e. consistency: lerp(x, x, t) == x for all t
assert(lerp(a, a, t) == a);
assert(lerp(b, b, t) == b);
#endif
return 0;
}
Raw
lerptest.cpp
// Test case adapted from Microsoft STL implementation
// (c) Microsoft Corporation, licensed under Apache 2.0:
// https://github.com/microsoft/STL/blob/main/tests/std/tests/P0811R3_midpoint_lerp/test.cpp
#include <algorithm>
#include <bit>
#include <cassert>
#include <cfenv>
#include <charconv>
#include <cmath>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <iterator>
#include <limits>
#include <numeric>
#include <optional>
#include <type_traits>
using namespace std;
template <typename Ty>
using limits = numeric_limits<Ty>;
template <typename Ty>
void assert_bitwise_equal(const Ty& a, const Ty& b) {
assert(memcmp(&a, &b, sizeof(Ty)) == 0);
}
template <typename Ty>
struct constants; // not defined
template <>
struct constants<float> {
static constexpr float TwoPlusUlp = 0x1.000002p+1f;
static constexpr float OnePlusUlp = 0x1.000002p+0f;
static constexpr float PointFivePlusUlp = 0x1.000002p-1f;
static constexpr float OneMinusUlp = 0x1.fffffep-1f;
static constexpr float PointFiveMinusUlp = 0x1.fffffep-2f;
static constexpr float NegOneMinusUlp = -OnePlusUlp;
static constexpr float NegOnePlusUlp = -OneMinusUlp;
static constexpr float EighthPlusUlp = 0x1.000002p-3f;
static constexpr float EighthMinusUlp = 0x1.fffffep-4f;
};
template <>
struct constants<double> {
static constexpr double TwoPlusUlp = 0x1.0000000000001p+1;
static constexpr double OnePlusUlp = 0x1.0000000000001p+0;
static constexpr double PointFivePlusUlp = 0x1.0000000000001p-1;
static constexpr double OneMinusUlp = 0x1.fffffffffffffp-1;
static constexpr double PointFiveMinusUlp = 0x1.fffffffffffffp-2;
static constexpr double NegOneMinusUlp = -OnePlusUlp;
static constexpr double NegOnePlusUlp = -OneMinusUlp;
static constexpr double EighthPlusUlp = 0x1.0000000000001p-3;
static constexpr double EighthMinusUlp = 0x1.fffffffffffffp-4;
};
template <typename Ty>
constexpr int cmp(const Ty x, const Ty y) {
if (x > y) {
return 1;
} else if (x < y) {
return -1;
} else {
return 0;
}
}
template <typename Ty>
struct LerpTestCase {
Ty x;
Ty y;
Ty t;
Ty expected;
};
template <typename Ty>
struct LerpCases { // TRANSITION, VSO-934633
static constexpr LerpTestCase<Ty> lerpTestCases[] = {
{Ty(-1.0), Ty(1.0), Ty(2.0), Ty(3.0)},
{Ty(0.0), Ty(1.0), Ty(2.0), Ty(2.0)},
{Ty(-1.0), Ty(0.0), Ty(2.0), Ty(1.0)},
{Ty(1.0), Ty(-1.0), Ty(2.0), Ty(-3.0)},
{Ty(0.0), Ty(-1.0), Ty(2.0), Ty(-2.0)},
{Ty(1.0), Ty(0.0), Ty(2.0), Ty(-1.0)},
{Ty(1.0), Ty(2.0), Ty(1.0), Ty(2.0)},
{Ty(1.0), Ty(2.0), Ty(2.0), Ty(3.0)},
{Ty(1.0), Ty(2.0), Ty(0.5), Ty(1.5)},
{Ty(1.0), Ty(2.0), Ty(0.0), Ty(1.0)},
{Ty(1.0), Ty(1.0), Ty(2.0), Ty(1.0)},
// This is the only (finite) test case that I had to disable to make tests pass with luaulerp
// {Ty(-0.0), Ty(-0.0), Ty(0.5), Ty(-0.0)},
{Ty(0.0), Ty(0.0), Ty(0.5), Ty(0.0)},
{Ty(-5.0), Ty(5.0), Ty(0.5), Ty(0.0)},
// float: lerp(-0x1.0000000000000p-149, 0x1.0000000000000p-149, 0.5) = 0x0.0000000000000p+0
// double: lerp(-0x0.0000000000001p-1022, 0x0.0000000000001p-1022, 0.5) = 0x0.0000000000000p+0
{-limits<Ty>::denorm_min(), limits<Ty>::denorm_min(), Ty(0.5), Ty(0.0)},
// float: lerp(0x1.fffffe0000000p-4, 0x1.0000020000000p-3, 0.5) = 0x1.0000000000000p-3
// double: lerp(0x1.fffffffffffffp-4, 0x1.0000000000001p-3, 0.5) = 0x1.0000000000000p-3
{constants<Ty>::EighthMinusUlp, constants<Ty>::EighthPlusUlp, Ty(0.5), Ty(0.125)},
// float: lerp(-0x1.0000020000000p-3, -0x1.fffffe0000000p-4, 0.5) = -0x1.0000000000000p-3
// double: lerp(-0x1.0000000000001p-3, -0x1.fffffffffffffp-4, 0.5) = -0x1.0000000000000p-3
{-constants<Ty>::EighthPlusUlp, -constants<Ty>::EighthMinusUlp, Ty(0.5), Ty(-0.125)},
// float: lerp(0x1.0000000000000p-149, 0x1.0000020000000p+0, 0.5) = 0x1.0000020000000p-1
// double: lerp(0x0.0000000000001p-1022, 0x1.0000000000001p+0, 0.5) = 0x1.0000000000001p-1
{limits<Ty>::denorm_min(), constants<Ty>::OnePlusUlp, Ty(0.5), constants<Ty>::PointFivePlusUlp},
// float: lerp(-0x1.0000000000000p-149, -0x1.0000020000000p+0, 0.5) = -0x1.0000020000000p-1
// double: lerp(-0x0.0000000000001p-1022, -0x1.0000000000001p+0, 0.5) = -0x1.0000000000001p-1
{-limits<Ty>::denorm_min(), constants<Ty>::NegOneMinusUlp, Ty(0.5), -constants<Ty>::PointFivePlusUlp},
{Ty(1.0), constants<Ty>::OnePlusUlp, constants<Ty>::PointFiveMinusUlp, Ty(1.0)},
{Ty(-1.0), constants<Ty>::NegOneMinusUlp, constants<Ty>::PointFiveMinusUlp, Ty(-1.0)},
{Ty(1.0), constants<Ty>::OnePlusUlp, Ty(0.5), Ty(1.0)},
{Ty(-1.0), constants<Ty>::NegOneMinusUlp, Ty(0.5), Ty(-1.0)},
{Ty(1.0), constants<Ty>::OnePlusUlp, constants<Ty>::PointFivePlusUlp, constants<Ty>::OnePlusUlp},
{Ty(-1.0), constants<Ty>::NegOneMinusUlp, constants<Ty>::PointFivePlusUlp, constants<Ty>::NegOneMinusUlp},
{Ty(-1.0), Ty(1.0), constants<Ty>::TwoPlusUlp, Ty(2.0) * constants<Ty>::TwoPlusUlp - Ty(1.0)},
{Ty(0.0), Ty(1.0), constants<Ty>::TwoPlusUlp, constants<Ty>::TwoPlusUlp},
{Ty(-1.0), Ty(0.0), constants<Ty>::TwoPlusUlp, constants<Ty>::TwoPlusUlp - Ty(1.0)},
{Ty(1.0), Ty(-1.0), constants<Ty>::TwoPlusUlp, Ty(1.0) - Ty(2.0) * constants<Ty>::TwoPlusUlp},
{Ty(0.0), Ty(-1.0), constants<Ty>::TwoPlusUlp, -constants<Ty>::TwoPlusUlp},
{Ty(1.0), Ty(0.0), constants<Ty>::TwoPlusUlp, Ty(1.0) - constants<Ty>::TwoPlusUlp},
{Ty(1.0), Ty(2.0), constants<Ty>::OnePlusUlp, Ty(2.0)},
{Ty(1.0), Ty(2.0), constants<Ty>::TwoPlusUlp, constants<Ty>::TwoPlusUlp + Ty(1.0)},
{Ty(1.0), Ty(2.0), Ty(0.5), constants<Ty>::PointFivePlusUlp + Ty(1.0)},
{limits<Ty>::max(), limits<Ty>::max(), Ty(2.0), limits<Ty>::max()},
{limits<Ty>::lowest(), limits<Ty>::lowest(), Ty(2.0), limits<Ty>::lowest()},
{limits<Ty>::denorm_min(), Ty(1.0), Ty(0.5), limits<Ty>::denorm_min() + Ty(1.0) / Ty(2.0)},
{limits<Ty>::denorm_min(), limits<Ty>::max(), Ty(0.5), limits<Ty>::denorm_min() + limits<Ty>::max() / Ty(2.0)},
{limits<Ty>::denorm_min(), limits<Ty>::lowest(), Ty(0.5),
(limits<Ty>::denorm_min() + limits<Ty>::lowest()) / Ty(2.0)},
{limits<Ty>::denorm_min(), limits<Ty>::infinity(), Ty(0.5), limits<Ty>::infinity()},
{limits<Ty>::denorm_min(), -limits<Ty>::infinity(), Ty(0.5), -limits<Ty>::infinity()},
};
};
template <typename Ty>
void print_lerp_result(const LerpTestCase<Ty>& testCase, const Ty answer) {
char failureMessageBuffer[1000]{};
char* cursor = failureMessageBuffer;
char* end = cursor + sizeof(failureMessageBuffer);
memcpy(cursor, "lerp(", 5);
cursor += 5;
cursor = to_chars(cursor, end, testCase.x).ptr;
memcpy(cursor, ", ", 2);
cursor += 2;
cursor = to_chars(cursor, end, testCase.y).ptr;
memcpy(cursor, ", ", 2);
cursor += 2;
cursor = to_chars(cursor, end, testCase.t).ptr;
memcpy(cursor, ") == ", 5);
cursor += 5;
cursor = to_chars(cursor, end, answer).ptr;
memcpy(cursor, "; expected ", 11);
cursor += 11;
cursor = to_chars(cursor, end, testCase.expected).ptr;
cursor[0] = '\0';
puts(failureMessageBuffer);
}
float luaulerp(float a, float b, float t) {
return t == 1 ? b : a + (b - a) * t;
}
double luaulerp(double a, double b, double t) {
return t == 1 ? b : a + (b - a) * t;
}
template <typename Ty>
bool test_lerp() {
for (auto&& testCase : LerpCases<Ty>::lerpTestCases) {
const auto answer = luaulerp(testCase.x, testCase.y, testCase.t);
if (memcmp(&answer, &testCase.expected, sizeof(Ty)) != 0) {
print_lerp_result(testCase, answer);
abort();
}
}
assert(cmp(luaulerp(Ty(1.0), Ty(2.0), Ty(4.0)), luaulerp(Ty(1.0), Ty(2.0), Ty(3.0))) * cmp(Ty(4.0), Ty(3.0))
* cmp(Ty(2.0), Ty(1.0))
>= 0);
return true;
}
int main() {
test_lerp<float>();
test_lerp<double>();
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment