24956271367

Committed 26 Apr 2026 12:04PM UTC coverage: 84.299% (-0.001%) from 84.3%

Build # 24956271367

Build Type

Pull #772

github

Committed by

web-flow

Commit Message

Merge 33d58a242 into a18a12f84

Pull Request Pull Request #772: feat(math): add constexpr sqrt to sw::math::constexpr_math

Coverage Stats

82 of 97 new or added lines in 2 files covered. (84.54%)

4 existing lines in 2 files now uncovered.

45431 of 53893 relevant lines covered (84.3%)

6353662.69 hits per line

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

76.92

/internal/constexpr_math/api/sqrt.cpp

// sqrt.cpp: regression for sw::math::constexpr_math::sqrt(float|double)
//
// 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 <bit>
#include <cmath>
#include <cstdint>
#include <iostream>
#include <iomanip>
#include <limits>

#include <math/constexpr_math.hpp>

namespace cm = sw::math::constexpr_math;

// Constexpr-friendly negative-zero predicate. `value == -0.0` is true for
// both +0 and -0 (signed zeros compare equal in IEEE-754), so == cannot
// validate sign preservation. std::signbit isn't constexpr until C++23 and
// Universal targets C++20, so we extract the sign bit via std::bit_cast.
constexpr bool is_negative_zero(double v) {
        return std::bit_cast<std::uint64_t>(v) == (std::uint64_t{1} << 63);
}
constexpr bool is_negative_zero(float v) {
        return std::bit_cast<std::uint32_t>(v) == (std::uint32_t{1} << 31);
}

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

// Perfect squares: must be bit-exact (Newton converges to the exact answer).
static_assert(cm::sqrt(0.0)    == 0.0,   "sqrt(0) == 0");
static_assert(cm::sqrt(1.0)    == 1.0,   "sqrt(1) == 1");
static_assert(cm::sqrt(4.0)    == 2.0,   "sqrt(4) == 2");
static_assert(cm::sqrt(9.0)    == 3.0,   "sqrt(9) == 3");
static_assert(cm::sqrt(16.0)   == 4.0,   "sqrt(16) == 4");
static_assert(cm::sqrt(25.0)   == 5.0,   "sqrt(25) == 5");
static_assert(cm::sqrt(100.0)  == 10.0,  "sqrt(100) == 10");
static_assert(cm::sqrt(10000.0) == 100.0, "sqrt(10000) == 100");
static_assert(cm::sqrt(0.25)   == 0.5,   "sqrt(0.25) == 0.5");

// Float overload, perfect squares
static_assert(cm::sqrt(1.0f)   == 1.0f,  "sqrtf(1) == 1");
static_assert(cm::sqrt(4.0f)   == 2.0f,  "sqrtf(4) == 2");
static_assert(cm::sqrt(81.0f)  == 9.0f,  "sqrtf(81) == 9");

// Non-perfect squares within a few ulp
constexpr double cx_sqrt2 = cm::sqrt(2.0);
static_assert(cx_sqrt2 > 1.4142135623 && cx_sqrt2 < 1.4142135624,
              "sqrt(2) ~= 1.4142135623730951");

constexpr double cx_sqrt3 = cm::sqrt(3.0);
static_assert(cx_sqrt3 > 1.7320508075 && cx_sqrt3 < 1.7320508076,
              "sqrt(3) ~= 1.7320508075688772");

// Cross-check with detail::SQRT2 constant
static_assert(cm::sqrt(2.0) > cm::detail::SQRT2 - 1e-15
           && cm::sqrt(2.0) < cm::detail::SQRT2 + 1e-15,
              "sqrt(2) matches detail::SQRT2 within machine epsilon");

// Special values
static_assert(cm::sqrt(std::numeric_limits<double>::infinity())
              == std::numeric_limits<double>::infinity(),
              "sqrt(+inf) == +inf");
static_assert(cm::sqrt(-1.0) != cm::sqrt(-1.0), "sqrt(-1) == NaN");
static_assert(cm::sqrt(-2.5f) != cm::sqrt(-2.5f), "sqrtf(-2.5) == NaN");
static_assert(cm::sqrt(-std::numeric_limits<double>::infinity())
              != cm::sqrt(-std::numeric_limits<double>::infinity()),
              "sqrt(-inf) == NaN");
static_assert(cm::sqrt(std::numeric_limits<double>::quiet_NaN())
              != cm::sqrt(std::numeric_limits<double>::quiet_NaN()),
              "sqrt(NaN) == NaN");

// Signed-zero preservation: sqrt(-0) == -0 per IEEE-754, NOT NaN.
// Use is_negative_zero so the sign bit is asserted explicitly. (==-0.0 alone
// is true for both +0 and -0 and so cannot detect a sign regression.)
static_assert(is_negative_zero(cm::sqrt(-0.0)),  "sqrt(-0) preserves sign bit");
static_assert(is_negative_zero(cm::sqrt(-0.0f)), "sqrtf(-0) preserves sign bit");

// Wide dynamic range
constexpr double cx_sqrt_big = cm::sqrt(1e100);
static_assert(cx_sqrt_big > 9.9999e49 && cx_sqrt_big < 1.0001e50,
              "sqrt(1e100) ~= 1e50");
constexpr double cx_sqrt_small = cm::sqrt(1e-100);
static_assert(cx_sqrt_small > 9.9999e-51 && cx_sqrt_small < 1.0001e-50,
              "sqrt(1e-100) ~= 1e-50");

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

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

        int errors = 0;
        auto check = [&](const char* name, double x, double our, double ref, double tol) {
                double err = (ref == 0.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 dynamic range.
        const double points[] = {
                1e-300, 1e-100, 1e-50, 1e-10, 1e-5, 1e-3, 0.1, 0.5, 0.999, 1.0, 1.001,
                1.5, 2.0, 3.14, 7.0, 10.0, 100.0, 1024.0, 1e3, 1e10, 1e50, 1e100, 1e300,
        };
        for (double x : points) {
                double our = cm::sqrt(x);
                double ref = std::sqrt(x);
                check("double sweep", x, our, ref, 1e-15);
        }

        // Subnormal sample (after normalization in sqrt's bit handling).
        {
                double sub = std::numeric_limits<double>::min() / 4.0;
                double our = cm::sqrt(sub);
                double ref = std::sqrt(sub);
                check("double subnormal", sub, our, ref, 1e-14);
        }

        // Round-trip stress: sqrt(x)^2 ~= x.
        const double rt_points[] = {
                1e-50, 1e-10, 1e-5, 0.5, 1.5, 2.0, 3.14, 100.0, 1e10, 1e50,
        };
        for (double x : rt_points) {
                double s = cm::sqrt(x);
                double rt = s * s;
                check("round-trip sqrt(x)^2", x, rt, x, 1e-14);
        }

        // Float sweep
        const float fpoints[] = {
                1e-30f, 1e-10f, 0.001f, 0.5f, 1.0f, 2.0f, 3.14f, 100.0f, 1e10f, 1e30f,
        };
        for (float x : fpoints) {
                float our = cm::sqrt(x);
                float ref = std::sqrt(x);
                double err = (ref == 0.0f) ? std::abs(our) : std::abs((our - ref) / ref);
                if (err > 1e-6) {
                        ++errors;
                        std::cout << "FAIL float sweep  x=" << x
                                  << "  our=" << our << "  ref=" << ref
                                  << "  rel-err=" << err << '\n';
                }
        }

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

1	// sqrt.cpp: regression for sw::math::constexpr_math::sqrt(float\|double)
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 <bit>
9	#include <cmath>
10	#include <cstdint>
11	#include <iostream>
12	#include <iomanip>
13	#include <limits>
14
15	#include <math/constexpr_math.hpp>
16
17	namespace cm = sw::math::constexpr_math;
18
19	// Constexpr-friendly negative-zero predicate. `value == -0.0` is true for
20	// both +0 and -0 (signed zeros compare equal in IEEE-754), so == cannot
21	// validate sign preservation. std::signbit isn't constexpr until C++23 and
22	// Universal targets C++20, so we extract the sign bit via std::bit_cast.
23	constexpr bool is_negative_zero(double v) {
24	return std::bit_cast<std::uint64_t>(v) == (std::uint64_t{1} << 63);
25	}
26	constexpr bool is_negative_zero(float v) {
27	return std::bit_cast<std::uint32_t>(v) == (std::uint32_t{1} << 31);
28	}
29
30	// ============================================================================
31	// Compile-time correctness via static_assert
32	// ============================================================================
33
34	// Perfect squares: must be bit-exact (Newton converges to the exact answer).
35	static_assert(cm::sqrt(0.0) == 0.0, "sqrt(0) == 0");
36	static_assert(cm::sqrt(1.0) == 1.0, "sqrt(1) == 1");
37	static_assert(cm::sqrt(4.0) == 2.0, "sqrt(4) == 2");
38	static_assert(cm::sqrt(9.0) == 3.0, "sqrt(9) == 3");
39	static_assert(cm::sqrt(16.0) == 4.0, "sqrt(16) == 4");
40	static_assert(cm::sqrt(25.0) == 5.0, "sqrt(25) == 5");
41	static_assert(cm::sqrt(100.0) == 10.0, "sqrt(100) == 10");
42	static_assert(cm::sqrt(10000.0) == 100.0, "sqrt(10000) == 100");
43	static_assert(cm::sqrt(0.25) == 0.5, "sqrt(0.25) == 0.5");
44
45	// Float overload, perfect squares
46	static_assert(cm::sqrt(1.0f) == 1.0f, "sqrtf(1) == 1");
47	static_assert(cm::sqrt(4.0f) == 2.0f, "sqrtf(4) == 2");
48	static_assert(cm::sqrt(81.0f) == 9.0f, "sqrtf(81) == 9");
49
50	// Non-perfect squares within a few ulp
51	constexpr double cx_sqrt2 = cm::sqrt(2.0);
52	static_assert(cx_sqrt2 > 1.4142135623 && cx_sqrt2 < 1.4142135624,
53	"sqrt(2) ~= 1.4142135623730951");
54
55	constexpr double cx_sqrt3 = cm::sqrt(3.0);
56	static_assert(cx_sqrt3 > 1.7320508075 && cx_sqrt3 < 1.7320508076,
57	"sqrt(3) ~= 1.7320508075688772");
58
59	// Cross-check with detail::SQRT2 constant
60	static_assert(cm::sqrt(2.0) > cm::detail::SQRT2 - 1e-15
61	&& cm::sqrt(2.0) < cm::detail::SQRT2 + 1e-15,
62	"sqrt(2) matches detail::SQRT2 within machine epsilon");
63
64	// Special values
65	static_assert(cm::sqrt(std::numeric_limits<double>::infinity())
66	== std::numeric_limits<double>::infinity(),
67	"sqrt(+inf) == +inf");
68	static_assert(cm::sqrt(-1.0) != cm::sqrt(-1.0), "sqrt(-1) == NaN");
69	static_assert(cm::sqrt(-2.5f) != cm::sqrt(-2.5f), "sqrtf(-2.5) == NaN");
70	static_assert(cm::sqrt(-std::numeric_limits<double>::infinity())
71	!= cm::sqrt(-std::numeric_limits<double>::infinity()),
72	"sqrt(-inf) == NaN");
73	static_assert(cm::sqrt(std::numeric_limits<double>::quiet_NaN())
74	!= cm::sqrt(std::numeric_limits<double>::quiet_NaN()),
75	"sqrt(NaN) == NaN");
76
77	// Signed-zero preservation: sqrt(-0) == -0 per IEEE-754, NOT NaN.
78	// Use is_negative_zero so the sign bit is asserted explicitly. (==-0.0 alone
79	// is true for both +0 and -0 and so cannot detect a sign regression.)
80	static_assert(is_negative_zero(cm::sqrt(-0.0)), "sqrt(-0) preserves sign bit");
81	static_assert(is_negative_zero(cm::sqrt(-0.0f)), "sqrtf(-0) preserves sign bit");
82
83	// Wide dynamic range
84	constexpr double cx_sqrt_big = cm::sqrt(1e100);
85	static_assert(cx_sqrt_big > 9.9999e49 && cx_sqrt_big < 1.0001e50,
86	"sqrt(1e100) ~= 1e50");
87	constexpr double cx_sqrt_small = cm::sqrt(1e-100);
88	static_assert(cx_sqrt_small > 9.9999e-51 && cx_sqrt_small < 1.0001e-50,
89	"sqrt(1e-100) ~= 1e-50");
90
91	// ============================================================================
92	// Runtime cross-check vs std::sqrt
93	// ============================================================================
94
95	int main() {	1✔
96	std::cout << "sw::math::constexpr_math::sqrt verification\n";	1✔
97
98	int errors = 0;	1✔
99	auto check = [&](const char* name, double x, double our, double ref, double tol) {	34✔
100	double err = (ref == 0.0) ? std::abs(our) : std::abs((our - ref) / ref);	34✔
101	if (err > tol) {	34✔
NEW 102	++errors;	×
103	std::cout << "FAIL " << name
NEW 104	<< " x=" << std::setprecision(17) << x	×
NEW 105	<< " our=" << our	×
NEW 106	<< " ref=" << ref	×
NEW 107	<< " rel-err=" << err << '\n';	×
108	}
109	};	35✔
110
111	// Sweep across a representative dynamic range.
112	const double points[] = {	1✔
113	1e-300, 1e-100, 1e-50, 1e-10, 1e-5, 1e-3, 0.1, 0.5, 0.999, 1.0, 1.001,
114	1.5, 2.0, 3.14, 7.0, 10.0, 100.0, 1024.0, 1e3, 1e10, 1e50, 1e100, 1e300,
115	};
116	for (double x : points) {	24✔
117	double our = cm::sqrt(x);	23✔
118	double ref = std::sqrt(x);	23✔
119	check("double sweep", x, our, ref, 1e-15);	23✔
120	}
121
122	// Subnormal sample (after normalization in sqrt's bit handling).
123	{
124	double sub = std::numeric_limits<double>::min() / 4.0;	1✔
125	double our = cm::sqrt(sub);	1✔
126	double ref = std::sqrt(sub);	1✔
127	check("double subnormal", sub, our, ref, 1e-14);	1✔
128	}
129
130	// Round-trip stress: sqrt(x)^2 ~= x.
131	const double rt_points[] = {	1✔
132	1e-50, 1e-10, 1e-5, 0.5, 1.5, 2.0, 3.14, 100.0, 1e10, 1e50,
133	};
134	for (double x : rt_points) {	11✔
135	double s = cm::sqrt(x);	10✔
136	double rt = s * s;	10✔
137	check("round-trip sqrt(x)^2", x, rt, x, 1e-14);	10✔
138	}
139
140	// Float sweep
141	const float fpoints[] = {	1✔
142	1e-30f, 1e-10f, 0.001f, 0.5f, 1.0f, 2.0f, 3.14f, 100.0f, 1e10f, 1e30f,
143	};
144	for (float x : fpoints) {	11✔
145	float our = cm::sqrt(x);	10✔
146	float ref = std::sqrt(x);	10✔
147	double err = (ref == 0.0f) ? std::abs(our) : std::abs((our - ref) / ref);	10✔
148	if (err > 1e-6) {	10✔
NEW 149	++errors;	×
NEW 150	std::cout << "FAIL float sweep x=" << x	×
NEW 151	<< " our=" << our << " ref=" << ref	×
NEW 152	<< " rel-err=" << err << '\n';	×
153	}
154	}
155
156	std::cout << "constexpr_math::sqrt: " << (errors == 0 ? "PASS" : "FAIL")	1✔
157	<< " (" << errors << " error" << (errors == 1 ? "" : "s") << ")\n";	1✔
158	return (errors == 0) ? 0 : 1;	1✔
159	}

stillwater-sc / universal / 24956271367

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