24957906871

Committed 26 Apr 2026 01:33PM UTC coverage: 84.282% (+0.008%) from 84.274%

Build # 24957906871

Build Type

push

github

Committed by

web-flow

Commit Message

feat(math): add constexpr exp (natural exponential) to sw::math::constexpr_math (#775)

* feat(math): add constexpr exp (natural exponential) to sw::math::constexpr_math

Adds constexpr<float|double> exp to the constexpr_math facility (#763
Epic). Tier-2 sub-issue: thin wrapper exp(x) = exp2(x * LOG2E) on top of
the merged exp2 machinery (#764 / PR #769).

Implementation: special-value cases (NaN, +/-inf) handled before the
multiply because clang's stricter constexpr evaluator rejects "arithmetic
produces a NaN" -- mirrors the same pattern used in log.hpp.

Verification:
  - 13 static_asserts: exp(0) == 1, exp(1) ~= e, exp(-1) ~= 1/e,
    exp(ln(2)) ~= 2 (cross-check via detail::LN2), log(exp(pi)) ~= pi
    (round-trip with log), all special values (NaN, +/-inf) for both
    float and double, overflow saturation (exp(710) -> +inf),
    underflow (exp(-746) -> 0)
  - Runtime sweep over 17 hand-picked points spanning [-700, 700]
  - Round-trip log(exp(x)) for 10 points
  - Float sweep over 8 points

Out-of-band accuracy stress (100K random samples in [-700, 700]):
  - max relative error 4.92e-14 = 221.57 x double epsilon
  - bounded by the x * LOG2E multiply: at large |x| (e.g. +/-700),
    the magnitude (~1010) amplifies sub-ulp rounding to ~3-5e-14 rel
  - comparable to typical std::exp transcendental accuracy
  - vastly better than lns/takum precision needs

cm_log2, cm_exp2, cm_log, cm_pow, cm_sqrt, cm_exp -- all 6 facility
tests PASS on gcc and clang.

This is the **last** function in Epic #763. The constexpr_math facility
is now complete: log2, exp2, log, exp, pow, sqrt all implemented.

Resolves #768
Closes Epic #763 (functionally complete; admin closure separate)

Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>

* test(math): add off-by-one guards for exp saturation boundaries (CodeRabbit)

Per CodeRabbit on PR #775: the existing static_asserts only check the
hard saturation points at +/-710 / -746, leaving the "just inside"... (continued)

Coverage Stats

35 of 40 new or added lines in 2 files covered. (87.5%)

3 existing lines in 1 file now uncovered.

45467 of 53946 relevant lines covered (84.28%)

6349277.67 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

83.33

/internal/constexpr_math/api/exp.cpp

// exp.cpp: regression for sw::math::constexpr_math::exp (natural exponential)
//
// Copyright (C) 2017 Stillwater Supercomputing, Inc.
// SPDX-License-Identifier: MIT
//
// This file is part of the universal numbers project, which is released under an MIT Open Source license.

#include <cmath>
#include <iostream>
#include <iomanip>
#include <limits>

#include <math/constexpr_math.hpp>

namespace cm = sw::math::constexpr_math;

// ============================================================================
// Compile-time correctness via static_assert
// ============================================================================

// exp(0) == 1 exactly (the exp2(0) == 1 base case propagates)
static_assert(cm::exp(0.0)  == 1.0,  "exp(0) == 1");
static_assert(cm::exp(0.0f) == 1.0f, "expf(0) == 1");

// exp(1) ~= e = 2.71828182845904523536
constexpr double cx_e = cm::exp(1.0);
static_assert(cx_e > 2.71828182 && cx_e < 2.71828183,
              "exp(1) ~= 2.71828182845904523536");

// exp(-1) ~= 1/e = 0.36787944117144232
constexpr double cx_inv_e = cm::exp(-1.0);
static_assert(cx_inv_e > 0.36787944 && cx_inv_e < 0.36787945,
              "exp(-1) ~= 0.367879441171442");

// exp(ln(2)) ~= 2  (cross-check via the LN2 constant in detail::)
constexpr double cx_2 = cm::exp(cm::detail::LN2);
static_assert(cx_2 > 1.9999999 && cx_2 < 2.0000001,
              "exp(ln(2)) ~= 2");

// log(exp(x)) ~= x  (round-trip with log -- the typical encode/decode loop)
constexpr double cx_rt_pi = cm::log(cm::exp(3.14159265358979323846));
static_assert(cx_rt_pi > 3.1415926 && cx_rt_pi < 3.1415927,
              "log(exp(pi)) ~= pi");

// Special values
static_assert(cm::exp(std::numeric_limits<double>::infinity())
              == std::numeric_limits<double>::infinity(),
              "exp(+inf) == +inf");
static_assert(cm::exp(-std::numeric_limits<double>::infinity()) == 0.0,
              "exp(-inf) == 0");
static_assert(cm::exp(std::numeric_limits<double>::quiet_NaN())
              != cm::exp(std::numeric_limits<double>::quiet_NaN()),
              "exp(NaN) == NaN");

// Float specials
static_assert(cm::exp(std::numeric_limits<float>::infinity())
              == std::numeric_limits<float>::infinity(),
              "expf(+inf) == +inf");
static_assert(cm::exp(-std::numeric_limits<float>::infinity()) == 0.0f,
              "expf(-inf) == 0");
static_assert(cm::exp(std::numeric_limits<float>::quiet_NaN())
              != cm::exp(std::numeric_limits<float>::quiet_NaN()),
              "expf(NaN) == NaN");

// Overflow / underflow saturation (at the bounds of double's exponent range)
//   exp(710) -> +inf  (since 710 * log2(e) > 1024)
//   exp(-746) -> 0    (since -746 * log2(e) < -1075)
static_assert(cm::exp(710.0)  == std::numeric_limits<double>::infinity(),
              "exp(710) saturates to +inf");
static_assert(cm::exp(-746.0) == 0.0, "exp(-746) underflows to 0");

// Off-by-one regression guards: pin values just inside the true saturation
// boundaries to catch a future tightening of the saturation thresholds.
// True top boundary: x ~= 709.78 (= 1024 / LOG2E). At x = 709,
//   x * LOG2E ~= 1022.83, comfortably inside exp2's finite range.
// True bottom boundary: x ~= -744.44 (= log(smallest subnormal)). At x = -744,
//   the result is a positive subnormal just above denorm_min.
static_assert(cm::exp(709.0) < std::numeric_limits<double>::infinity()
           && cm::exp(709.0) > 0.0,
              "exp(709) is finite and positive (off-by-one guard for top saturation)");
static_assert(cm::exp(-744.0) > 0.0,
              "exp(-744) is positive subnormal (off-by-one guard for bottom underflow)");

// ============================================================================
// Runtime cross-check vs std::exp
// ============================================================================

int main() {
        std::cout << "sw::math::constexpr_math::exp verification\n";

        int errors = 0;
        // Generic lambda (C++20 auto-parameters) so float and double sweeps share
        // the same failure-printing path with one tolerance type per call.
        auto check = [&](const char* name, auto x, auto our, auto ref, auto tol) {
                auto err = (ref == decltype(ref){0}) ? std::abs(our)
                                                     : std::abs((our - ref) / ref);
                if (err > tol) {
                        ++errors;
                        std::cout << "FAIL " << name
                                  << "  x=" << std::setprecision(17) << x
                                  << "  our=" << our
                                  << "  ref=" << ref
                                  << "  rel-err=" << err << '\n';
                }
        };

        // Sweep across a representative range.
        const double points[] = {
                -700.0, -100.0, -10.0, -3.0, -1.5, -0.5, -0.0001,
                 0.0,    0.0001, 0.5,   1.5,  3.0,   10.0, 100.0, 700.0,
                 0.6931471805599453,  // ln(2) -> should give 2
                 2.302585092994046,   // ln(10) -> should give 10
        };
        for (double x : points) {
                double our = cm::exp(x);
                double ref = std::exp(x);
                // At extreme |x| (e.g. +/-700), the x * LOG2E multiply is sub-ulp
                // accurate but the magnitude (~1010) amplifies the rounding error to
                // ~3e-14. 1e-13 leaves comfortable margin without masking algorithmic
                // regressions.
                check("double sweep", x, our, ref, 1e-13);
        }

        // Round-trip: log(exp(x)) ~= x
        const double rt_points[] = {
                -50.0, -10.0, -1.0, -0.1, 0.0, 0.1, 1.0, 10.0, 50.0, 100.0,
        };
        for (double x : rt_points) {
                double rt = cm::log(cm::exp(x));
                check("round-trip log(exp(x))", x, rt, x, 1e-13);
        }

        // Float sweep -- shares the generic check lambda above.
        const float fpoints[] = {
                -80.0f, -10.0f, -1.0f, 0.0f, 1.0f, 2.71828183f, 10.0f, 80.0f,
        };
        for (float x : fpoints) {
                float our = cm::exp(x);
                float ref = std::exp(x);
                check("float sweep", x, our, ref, 1e-6f);
        }

        std::cout << "constexpr_math::exp: " << (errors == 0 ? "PASS" : "FAIL")
                  << " (" << errors << " error" << (errors == 1 ? "" : "s") << ")\n";
        return (errors == 0) ? 0 : 1;
}

1	// exp.cpp: regression for sw::math::constexpr_math::exp (natural exponential)
2	//
3	// Copyright (C) 2017 Stillwater Supercomputing, Inc.
4	// SPDX-License-Identifier: MIT
5	//
6	// This file is part of the universal numbers project, which is released under an MIT Open Source license.
7
8	#include <cmath>
9	#include <iostream>
10	#include <iomanip>
11	#include <limits>
12
13	#include <math/constexpr_math.hpp>
14
15	namespace cm = sw::math::constexpr_math;
16
17	// ============================================================================
18	// Compile-time correctness via static_assert
19	// ============================================================================
20
21	// exp(0) == 1 exactly (the exp2(0) == 1 base case propagates)
22	static_assert(cm::exp(0.0) == 1.0, "exp(0) == 1");
23	static_assert(cm::exp(0.0f) == 1.0f, "expf(0) == 1");
24
25	// exp(1) ~= e = 2.71828182845904523536
26	constexpr double cx_e = cm::exp(1.0);
27	static_assert(cx_e > 2.71828182 && cx_e < 2.71828183,
28	"exp(1) ~= 2.71828182845904523536");
29
30	// exp(-1) ~= 1/e = 0.36787944117144232
31	constexpr double cx_inv_e = cm::exp(-1.0);
32	static_assert(cx_inv_e > 0.36787944 && cx_inv_e < 0.36787945,
33	"exp(-1) ~= 0.367879441171442");
34
35	// exp(ln(2)) ~= 2 (cross-check via the LN2 constant in detail::)
36	constexpr double cx_2 = cm::exp(cm::detail::LN2);
37	static_assert(cx_2 > 1.9999999 && cx_2 < 2.0000001,
38	"exp(ln(2)) ~= 2");
39
40	// log(exp(x)) ~= x (round-trip with log -- the typical encode/decode loop)
41	constexpr double cx_rt_pi = cm::log(cm::exp(3.14159265358979323846));
42	static_assert(cx_rt_pi > 3.1415926 && cx_rt_pi < 3.1415927,
43	"log(exp(pi)) ~= pi");
44
45	// Special values
46	static_assert(cm::exp(std::numeric_limits<double>::infinity())
47	== std::numeric_limits<double>::infinity(),
48	"exp(+inf) == +inf");
49	static_assert(cm::exp(-std::numeric_limits<double>::infinity()) == 0.0,
50	"exp(-inf) == 0");
51	static_assert(cm::exp(std::numeric_limits<double>::quiet_NaN())
52	!= cm::exp(std::numeric_limits<double>::quiet_NaN()),
53	"exp(NaN) == NaN");
54
55	// Float specials
56	static_assert(cm::exp(std::numeric_limits<float>::infinity())
57	== std::numeric_limits<float>::infinity(),
58	"expf(+inf) == +inf");
59	static_assert(cm::exp(-std::numeric_limits<float>::infinity()) == 0.0f,
60	"expf(-inf) == 0");
61	static_assert(cm::exp(std::numeric_limits<float>::quiet_NaN())
62	!= cm::exp(std::numeric_limits<float>::quiet_NaN()),
63	"expf(NaN) == NaN");
64
65	// Overflow / underflow saturation (at the bounds of double's exponent range)
66	// exp(710) -> +inf (since 710 * log2(e) > 1024)
67	// exp(-746) -> 0 (since -746 * log2(e) < -1075)
68	static_assert(cm::exp(710.0) == std::numeric_limits<double>::infinity(),
69	"exp(710) saturates to +inf");
70	static_assert(cm::exp(-746.0) == 0.0, "exp(-746) underflows to 0");
71
72	// Off-by-one regression guards: pin values just inside the true saturation
73	// boundaries to catch a future tightening of the saturation thresholds.
74	// True top boundary: x ~= 709.78 (= 1024 / LOG2E). At x = 709,
75	// x * LOG2E ~= 1022.83, comfortably inside exp2's finite range.
76	// True bottom boundary: x ~= -744.44 (= log(smallest subnormal)). At x = -744,
77	// the result is a positive subnormal just above denorm_min.
78	static_assert(cm::exp(709.0) < std::numeric_limits<double>::infinity()
79	&& cm::exp(709.0) > 0.0,
80	"exp(709) is finite and positive (off-by-one guard for top saturation)");
81	static_assert(cm::exp(-744.0) > 0.0,
82	"exp(-744) is positive subnormal (off-by-one guard for bottom underflow)");
83
84	// ============================================================================
85	// Runtime cross-check vs std::exp
86	// ============================================================================
87
88	int main() {	1✔
89	std::cout << "sw::math::constexpr_math::exp verification\n";	1✔
90
91	int errors = 0;	1✔
92	// Generic lambda (C++20 auto-parameters) so float and double sweeps share
93	// the same failure-printing path with one tolerance type per call.
94	auto check = [&](const char* name, auto x, auto our, auto ref, auto tol) {	35✔
95	auto err = (ref == decltype(ref){0}) ? std::abs(our)	35✔
96	: std::abs((our - ref) / ref);	34✔
97	if (err > tol) {	35✔
NEW 98	++errors;	×
99	std::cout << "FAIL " << name
NEW 100	<< " x=" << std::setprecision(17) << x	×
NEW 101	<< " our=" << our	×
NEW 102	<< " ref=" << ref	×
NEW 103	<< " rel-err=" << err << '\n';	×
104	}
105	};	36✔
106
107	// Sweep across a representative range.
108	const double points[] = {	1✔
109	-700.0, -100.0, -10.0, -3.0, -1.5, -0.5, -0.0001,
110	0.0, 0.0001, 0.5, 1.5, 3.0, 10.0, 100.0, 700.0,
111	0.6931471805599453, // ln(2) -> should give 2
112	2.302585092994046, // ln(10) -> should give 10
113	};
114	for (double x : points) {	18✔
115	double our = cm::exp(x);	17✔
116	double ref = std::exp(x);	17✔
117	// At extreme \|x\| (e.g. +/-700), the x * LOG2E multiply is sub-ulp
118	// accurate but the magnitude (~1010) amplifies the rounding error to
119	// ~3e-14. 1e-13 leaves comfortable margin without masking algorithmic
120	// regressions.
121	check("double sweep", x, our, ref, 1e-13);	17✔
122	}
123
124	// Round-trip: log(exp(x)) ~= x
125	const double rt_points[] = {	1✔
126	-50.0, -10.0, -1.0, -0.1, 0.0, 0.1, 1.0, 10.0, 50.0, 100.0,
127	};
128	for (double x : rt_points) {	11✔
129	double rt = cm::log(cm::exp(x));	10✔
130	check("round-trip log(exp(x))", x, rt, x, 1e-13);	10✔
131	}
132
133	// Float sweep -- shares the generic check lambda above.
134	const float fpoints[] = {	1✔
135	-80.0f, -10.0f, -1.0f, 0.0f, 1.0f, 2.71828183f, 10.0f, 80.0f,
136	};
137	for (float x : fpoints) {	9✔
138	float our = cm::exp(x);	8✔
139	float ref = std::exp(x);	8✔
140	check("float sweep", x, our, ref, 1e-6f);	8✔
141	}
142
143	std::cout << "constexpr_math::exp: " << (errors == 0 ? "PASS" : "FAIL")	1✔
144	<< " (" << errors << " error" << (errors == 1 ? "" : "s") << ")\n";	1✔
145	return (errors == 0) ? 0 : 1;	1✔
146	}

stillwater-sc / universal / 24957906871

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous