24956675961

Committed 26 Apr 2026 12:30PM UTC coverage: 84.313% (+0.01%) from 84.3%

Build # 24956675961

Build Type

push

github

Committed by

web-flow

Commit Message

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

* feat(math): add constexpr sqrt to sw::math::constexpr_math

Adds constexpr<float|double> sqrt to the constexpr_math facility (#763
Epic, #423/PR #762 substrate). Tier-2 sub-issue: independent of all
other constexpr_math functions, used for normalization, range constants,
and complex magnitude.

Algorithm:
  Decompose x = 2^E * M with M in [1, 2). If E is odd, absorb a factor
  of 2 into M and decrement E so it becomes even. Then
    sqrt(x) = sqrt(M) * 2^(E/2)
  with 2^(E/2) constructed exactly via detail::pow2 and sqrt(M) computed
  via Newton's iteration y_{n+1} = 0.5 * (y_n + M/y_n).

  Newton has quadratic convergence: starting from y_0 = 1.5 (worst-case
  ~2 bits accurate over M in [1, 4)), each iteration roughly doubles
  the bit-precision (2 -> 4 -> 8 -> 16 -> 32 -> 64 -> 128). Six
  iterations for double, four for float -- both well past their
  significand widths.

Special-value handling per IEEE-754:
  sqrt(NaN)   -> NaN
  sqrt(x<0)   -> NaN  (including -inf)
  sqrt(+/-0)  -> +/-0  (sign preserved)
  sqrt(+inf)  -> +inf

Self-contained: only depends on detail::pow2 (already in the facility)
and std::bit_cast.

Verification:
  - 22 static_asserts: perfect squares (bit-exact via Newton convergence
    to fixed point), sqrt(2) and sqrt(3) tolerance, cross-check against
    detail::SQRT2, special values (NaN/+/-inf/+/-0), wide dynamic range
  - Runtime sweep over 23 hand-picked points + subnormal sample +
    round-trip stress (sqrt(x)^2 ~= x)
  - Float sweep over 10 points

Out-of-band accuracy stress (100K random samples in [1e-100, 1e100]):
  max relative error 2.22e-16 = 1.00 x double epsilon (1 ulp)
  -- equivalent to std::sqrt accuracy

cm_log2, cm_exp2, cm_log, cm_pow, cm_sqrt all PASS on gcc and clang.

Resolves #767
Relates to #763 (Epic)

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

* test(math): pin sqrt(-0) sign preservation explicitly (Code... (continued)

Coverage Stats

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

45439 of 53893 relevant lines covered (84.31%)

6358923.24 hits per line

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

89.66

/include/sw/math/constexpr_math/sqrt.hpp

#pragma once
// sqrt.hpp: constexpr square root for IEEE-754 float and 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.
//
// -----------------------------------------------------------------------------
// Algorithm
// -----------------------------------------------------------------------------
// Decompose x = 2^E * M with M in [1, 2). If E is odd, absorb a factor of 2
// into M (giving M in [2, 4)) and decrement E so it becomes even. Then
//
//     sqrt(x) = sqrt(M) * 2^(E/2)
//
// where 2^(E/2) is constructed exactly via detail::pow2 and sqrt(M) is
// computed via Newton's iteration
//
//     y_{n+1} = 0.5 * (y_n + M / y_n)
//
// which has quadratic convergence. From a constant initial guess of 1.5
// (worst-case ~2 bits accurate over M in [1, 4)), six iterations reach
// ~96 bits of precision -- well past double's 53-bit significand. Two extra
// over the strict minimum keeps us safe against the rounding mode of any
// individual multiply/add.
//
// Self-contained: only depends on detail::pow2 (already in the facility) and
// std::bit_cast.

#include <bit>
#include <cstdint>
#include <limits>

#include <math/constexpr_math/detail.hpp>

namespace sw { namespace math { namespace constexpr_math {

// constexpr sqrt(double): returns the square root of a.
// Special values follow IEEE-754 conventions:
//   sqrt(NaN)   -> NaN
//   sqrt(x<0)   -> NaN  (including -inf)
//   sqrt(+/-0)  -> +/-0  (sign preserved)
//   sqrt(+inf)  -> +inf
constexpr double sqrt(double a) {
        if (a != a) return a;                                                // NaN propagation
        if (a < 0.0) return std::numeric_limits<double>::quiet_NaN();
        if (a == 0.0) return a;                                              // +/-0 preserved
        if (a == std::numeric_limits<double>::infinity()) return a;

        // Decompose a = 2^E * M with M in [1, 2). Mirrors log2/exp2's bit handling.
        std::uint64_t bits = std::bit_cast<std::uint64_t>(a);
        int exp_field = static_cast<int>((bits >> 52) & 0x7FFu);

        int subnormal_shift = 0;
        if (exp_field == 0) {
                std::uint64_t mant = bits & ((std::uint64_t{1} << 52) - 1);
                while ((mant & (std::uint64_t{1} << 52)) == 0) {
                        mant <<= 1;
                        ++subnormal_shift;
                }
                bits = mant;
                exp_field = 1;
        }

        int E = exp_field - 1023 - subnormal_shift;
        std::uint64_t m_bits = (bits & ((std::uint64_t{1} << 52) - 1))
                               | (std::uint64_t{1023} << 52);
        double M = std::bit_cast<double>(m_bits);                            // M in [1, 2)

        // If E is odd, absorb a factor of 2 into M so 2^(E/2) is an integer power.
        // Bitwise & 1 correctly identifies odd magnitudes for both signs in two's
        // complement; E -= 1 brings odd-positive down to even and odd-negative
        // further down to even.
        if (E & 1) {
                M *= 2.0;
                E -= 1;
        }
        // Now E is even and M is in [1, 4).

        // Newton's iteration for sqrt(M). Quadratic convergence: starting from
        // 1.5 (worst-case error ~2 bits over [1, 4)), each iteration roughly
        // doubles the bit-precision: 2 -> 4 -> 8 -> 16 -> 32 -> 64 -> 128.
        // Six iterations reach well past double's 53-bit significand.
        double y = 1.5;
        y = 0.5 * (y + M / y);
        y = 0.5 * (y + M / y);
        y = 0.5 * (y + M / y);
        y = 0.5 * (y + M / y);
        y = 0.5 * (y + M / y);
        y = 0.5 * (y + M / y);

        return y * detail::pow2(E / 2);
}

// constexpr sqrt(float): same algorithm at lower precision.
// 4 Newton iterations from 1.5 reach ~32 bits, well past float's 24-bit
// significand.
constexpr float sqrt(float a) {
        if (a != a) return a;
        if (a < 0.0f) return std::numeric_limits<float>::quiet_NaN();
        if (a == 0.0f) return a;
        if (a == std::numeric_limits<float>::infinity()) return a;

        std::uint32_t bits = std::bit_cast<std::uint32_t>(a);
        int exp_field = static_cast<int>((bits >> 23) & 0xFFu);

        int subnormal_shift = 0;
        if (exp_field == 0) {
                std::uint32_t mant = bits & ((std::uint32_t{1} << 23) - 1);
                while ((mant & (std::uint32_t{1} << 23)) == 0) {
                        mant <<= 1;
                        ++subnormal_shift;
                }
                bits = mant;
                exp_field = 1;
        }

        int E = exp_field - 127 - subnormal_shift;
        std::uint32_t m_bits = (bits & ((std::uint32_t{1} << 23) - 1))
                               | (std::uint32_t{127} << 23);
        float M = std::bit_cast<float>(m_bits);

        if (E & 1) {
                M *= 2.0f;
                E -= 1;
        }

        float y = 1.5f;
        y = 0.5f * (y + M / y);
        y = 0.5f * (y + M / y);
        y = 0.5f * (y + M / y);
        y = 0.5f * (y + M / y);

        return y * detail::pow2f(E / 2);
}

}}}  // namespace sw::math::constexpr_math

1	#pragma once
2	// sqrt.hpp: constexpr square root for IEEE-754 float and double
3	//
4	// Copyright (C) 2017 Stillwater Supercomputing, Inc.
5	// SPDX-License-Identifier: MIT
6	//
7	// This file is part of the universal numbers project, which is released under an MIT Open Source license.
8	//
9	// -----------------------------------------------------------------------------
10	// Algorithm
11	// -----------------------------------------------------------------------------
12	// Decompose x = 2^E * M with M in [1, 2). If E is odd, absorb a factor of 2
13	// into M (giving M in [2, 4)) and decrement E so it becomes even. Then
14	//
15	// sqrt(x) = sqrt(M) * 2^(E/2)
16	//
17	// where 2^(E/2) is constructed exactly via detail::pow2 and sqrt(M) is
18	// computed via Newton's iteration
19	//
20	// y_{n+1} = 0.5 * (y_n + M / y_n)
21	//
22	// which has quadratic convergence. From a constant initial guess of 1.5
23	// (worst-case ~2 bits accurate over M in [1, 4)), six iterations reach
24	// ~96 bits of precision -- well past double's 53-bit significand. Two extra
25	// over the strict minimum keeps us safe against the rounding mode of any
26	// individual multiply/add.
27	//
28	// Self-contained: only depends on detail::pow2 (already in the facility) and
29	// std::bit_cast.
30
31	#include <bit>
32	#include <cstdint>
33	#include <limits>
34
35	#include <math/constexpr_math/detail.hpp>
36
37	namespace sw { namespace math { namespace constexpr_math {
38
39	// constexpr sqrt(double): returns the square root of a.
40	// Special values follow IEEE-754 conventions:
41	// sqrt(NaN) -> NaN
42	// sqrt(x<0) -> NaN (including -inf)
43	// sqrt(+/-0) -> +/-0 (sign preserved)
44	// sqrt(+inf) -> +inf
45	constexpr double sqrt(double a) {	34✔
46	if (a != a) return a; // NaN propagation	34✔
47	if (a < 0.0) return std::numeric_limits<double>::quiet_NaN();	34✔
48	if (a == 0.0) return a; // +/-0 preserved	34✔
49	if (a == std::numeric_limits<double>::infinity()) return a;	34✔
50
51	// Decompose a = 2^E * M with M in [1, 2). Mirrors log2/exp2's bit handling.
52	std::uint64_t bits = std::bit_cast<std::uint64_t>(a);	34✔
53	int exp_field = static_cast<int>((bits >> 52) & 0x7FFu);	34✔
54
55	int subnormal_shift = 0;	34✔
56	if (exp_field == 0) {	34✔
57	std::uint64_t mant = bits & ((std::uint64_t{1} << 52) - 1);	1✔
58	while ((mant & (std::uint64_t{1} << 52)) == 0) {	3✔
59	mant <<= 1;	2✔
60	++subnormal_shift;	2✔
61	}
62	bits = mant;	1✔
63	exp_field = 1;	1✔
64	}
65
66	int E = exp_field - 1023 - subnormal_shift;	34✔
67	std::uint64_t m_bits = (bits & ((std::uint64_t{1} << 52) - 1))	34✔
68	\| (std::uint64_t{1023} << 52);	34✔
69	double M = std::bit_cast<double>(m_bits); // M in [1, 2)	34✔
70
71	// If E is odd, absorb a factor of 2 into M so 2^(E/2) is an integer power.
72	// Bitwise & 1 correctly identifies odd magnitudes for both signs in two's
73	// complement; E -= 1 brings odd-positive down to even and odd-negative
74	// further down to even.
75	if (E & 1) {	34✔
76	M *= 2.0;	17✔
77	E -= 1;	17✔
78	}
79	// Now E is even and M is in [1, 4).
80
81	// Newton's iteration for sqrt(M). Quadratic convergence: starting from
82	// 1.5 (worst-case error ~2 bits over [1, 4)), each iteration roughly
83	// doubles the bit-precision: 2 -> 4 -> 8 -> 16 -> 32 -> 64 -> 128.
84	// Six iterations reach well past double's 53-bit significand.
85	double y = 1.5;	34✔
86	y = 0.5 * (y + M / y);	34✔
87	y = 0.5 * (y + M / y);	34✔
88	y = 0.5 * (y + M / y);	34✔
89	y = 0.5 * (y + M / y);	34✔
90	y = 0.5 * (y + M / y);	34✔
91	y = 0.5 * (y + M / y);	34✔
92
93	return y * detail::pow2(E / 2);	34✔
94	}
95
96	// constexpr sqrt(float): same algorithm at lower precision.
97	// 4 Newton iterations from 1.5 reach ~32 bits, well past float's 24-bit
98	// significand.
99	constexpr float sqrt(float a) {	10✔
100	if (a != a) return a;	10✔
101	if (a < 0.0f) return std::numeric_limits<float>::quiet_NaN();	10✔
102	if (a == 0.0f) return a;	10✔
103	if (a == std::numeric_limits<float>::infinity()) return a;	10✔
104
105	std::uint32_t bits = std::bit_cast<std::uint32_t>(a);	10✔
106	int exp_field = static_cast<int>((bits >> 23) & 0xFFu);	10✔
107
108	int subnormal_shift = 0;	10✔
109	if (exp_field == 0) {	10✔
NEW 110	std::uint32_t mant = bits & ((std::uint32_t{1} << 23) - 1);	×
NEW 111	while ((mant & (std::uint32_t{1} << 23)) == 0) {	×
NEW 112	mant <<= 1;	×
NEW 113	++subnormal_shift;	×
114	}
NEW 115	bits = mant;	×
NEW 116	exp_field = 1;	×
117	}
118
119	int E = exp_field - 127 - subnormal_shift;	10✔
120	std::uint32_t m_bits = (bits & ((std::uint32_t{1} << 23) - 1))	10✔
121	\| (std::uint32_t{127} << 23);	10✔
122	float M = std::bit_cast<float>(m_bits);	10✔
123
124	if (E & 1) {	10✔
125	M *= 2.0f;	5✔
126	E -= 1;	5✔
127	}
128
129	float y = 1.5f;	10✔
130	y = 0.5f * (y + M / y);	10✔
131	y = 0.5f * (y + M / y);	10✔
132	y = 0.5f * (y + M / y);	10✔
133	y = 0.5f * (y + M / y);	10✔
134
135	return y * detail::pow2f(E / 2);	10✔
136	}
137
138	}}} // namespace sw::math::constexpr_math

stillwater-sc / universal / 24956675961

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