Build has been canceled!

24955899585

Committed 26 Apr 2026 11:50AM UTC coverage: 84.3% (-0.04%) from 84.337%

Build # 24955899585

Build Type

push

github

Committed by

web-flow

Commit Message

feat(math): add constexpr pow to sw::math::constexpr_math (#770)

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

Adds constexpr pow with two overload families to the constexpr_math
facility (#763 Epic, #423/PR #762, #764/PR #769 substrates). Tier-1
sub-issue: unblocks lns saturation/range constants (consumed by #722).

Two overloads per type (4 total: float and double):

1. Integer-exponent fast path: pow(T, int)
   - Fast exponentiation by squaring; pure integer logic.
   - O(log |n|) multiplications, no transcendentals.
   - Negative exponent via 1 / pow(x, -n) with safe long-long promotion
     to avoid INT_MIN overflow on negation.

2. General overload: pow(T, T)
   - Dispatches to the integer fast path when the floating-point
     exponent is exactly an integer (typical for compile-time constants).
   - For x > 0, finite y: pow(x, y) = exp2(y * log2(x)).
   - For x < 0, integer y: pow(|x|, y) with sign from y's parity.
   - For x < 0, non-integer y: NaN.
   - Full IEEE-754 special-case table per C99 7.12.7.4 (~20 cases):
     pow(x, 0) -> 1 (incl. NaN), pow(1, y) -> 1 (incl. NaN), pow(NaN, y),
     pow(x, NaN), pow(+/-0, y) by sign and y-parity, pow(+/-inf, y),
     pow(x, +/-inf) by |x| relative to 1, pow(-1, +/-inf), etc.

New helpers in detail.hpp:
  - is_integer(T)       -- exact-integer test, safe for huge doubles
                           (above 2^53 every value is necessarily integer)
  - is_odd_integer(T)   -- for negative-base sign rule
  - pow_by_squaring<T>  -- the squaring loop, used by both overloads

Test coverage in cm_pow:
  - 38 static_asserts: integer fast path positive/negative, negative
    base parity, all C99 7.12.7.4 special cases (NaN/inf/+-0 mixes)
  - Acceptance forms: pow(2, 10) == 1024, pow(2, 0.5) ~= sqrt(2),
    pow(8, 1/3) ~= 2
  - Runtime sweep over 21 (base, exponent) pairs vs std::pow
  - Float sweep over 7 (base, exponent) pairs

Out-of-band accuracy stress (100K random samples, positive base in
[... (continued)

Coverage Stats

109 of 155 new or added lines in 3 files covered. (70.32%)

6 existing lines in 2 files now uncovered.

45350 of 53796 relevant lines covered (84.3%)

6366762.75 hits per line

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

62.79

/include/sw/math/constexpr_math/pow.hpp

#pragma once
// pow.hpp: constexpr power function 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.
//
// -----------------------------------------------------------------------------
// Two overload families
// -----------------------------------------------------------------------------
//
// 1. Integer-exponent fast path: pow(T, int)
//      - Fast exponentiation by squaring; pure integer logic.
//      - O(log |n|) multiplications, no transcendentals.
//      - Selected when the user knows the exponent is an integer.
//
// 2. General overload: pow(T, T)
//      - First tries the integer fast path when the floating-point exponent
//        is exactly an integer (typical for compile-time constants).
//      - Otherwise computes pow(x, y) = exp2(y * log2(x)) for x > 0.
//      - Negative base with integer exponent: pow(|x|, y) with sign from y.
//      - Negative base with non-integer exponent: NaN.
//      - Full IEEE-754 special-case table per C99 7.12.7.4.
//
// -----------------------------------------------------------------------------
// IEEE-754 special-value handling (C99 7.12.7.4)
// -----------------------------------------------------------------------------
//   pow(x,    0)           -> 1                            (any x including NaN, inf)
//   pow(1,    y)           -> 1                            (any y including NaN, inf)
//   pow(NaN,  y)           -> NaN                          (y != 0)
//   pow(x,    NaN)         -> NaN                          (x != 1)
//   pow(+/-0, y<0, odd int)-> +/-inf  (signed)
//   pow(+/-0, y<0, else)   -> +inf
//   pow(+/-0, y>0, odd int)-> +/-0    (signed)
//   pow(+/-0, y>0, else)   -> +0
//   pow(-1,   +/-inf)      -> 1
//   pow(|x|<1,+inf)        -> +0
//   pow(|x|<1,-inf)        -> +inf
//   pow(|x|>1,+inf)        -> +inf
//   pow(|x|>1,-inf)        -> +0
//   pow(+inf, y<0)         -> +0
//   pow(+inf, y>0)         -> +inf
//   pow(-inf, y<0, odd int)-> -0
//   pow(-inf, y<0, else)   -> +0
//   pow(-inf, y>0, odd int)-> -inf
//   pow(-inf, y>0, else)   -> +inf
//   pow(x<0,  finite non-integer y) -> NaN

#include <limits>

#include <math/constexpr_math/detail.hpp>
#include <math/constexpr_math/exp2.hpp>
#include <math/constexpr_math/log2.hpp>

namespace sw { namespace math { namespace constexpr_math {

// -----------------------------------------------------------------------------
// Integer-exponent fast path
// -----------------------------------------------------------------------------

constexpr double pow(double x, int n) {
        if (n == 0) return 1.0;
        if (n < 0) {
                // Promote to long long to safely negate INT_MIN.
                unsigned long long exp = static_cast<unsigned long long>(-static_cast<long long>(n));
                return 1.0 / detail::pow_by_squaring(x, exp);
        }
        return detail::pow_by_squaring(x, static_cast<unsigned long long>(n));
}

constexpr float pow(float x, int n) {
        if (n == 0) return 1.0f;
        if (n < 0) {
                unsigned long long exp = static_cast<unsigned long long>(-static_cast<long long>(n));
                return 1.0f / detail::pow_by_squaring(x, exp);
        }
        return detail::pow_by_squaring(x, static_cast<unsigned long long>(n));
}

// -----------------------------------------------------------------------------
// General overload
// -----------------------------------------------------------------------------

constexpr double pow(double x, double y) {
        // Cases that win regardless of NaN: pow(x, 0) and pow(1, y).
        if (y == 0.0) return 1.0;
        if (x == 1.0) return 1.0;

        // NaN propagation (after the two NaN-overriding cases above).
        if (x != x || y != y) return std::numeric_limits<double>::quiet_NaN();

        const double pinf = std::numeric_limits<double>::infinity();

        // Special case: pow(-1, +/-inf) == 1
        if (x == -1.0 && (y == pinf || y == -pinf)) return 1.0;

        // y is +/- infinity
        if (y == pinf) {
                double ax = (x < 0.0) ? -x : x;
                if (ax < 1.0) return 0.0;
                // |x| > 1 (==1 was handled by pow(1,y) early-out and x==-1 above)
                return pinf;
        }
        if (y == -pinf) {
                double ax = (x < 0.0) ? -x : x;
                if (ax < 1.0) return pinf;
                return 0.0;
        }

        // x is +/- infinity (y is finite, non-zero)
        if (x == pinf) {
                return (y < 0.0) ? 0.0 : pinf;
        }
        if (x == -pinf) {
                bool y_odd = detail::is_odd_integer(y);
                if (y < 0.0) return y_odd ? -0.0 : 0.0;
                return y_odd ? -pinf : pinf;
        }

        // x is +/- 0 (y is finite, non-zero, not infinity)
        if (x == 0.0) {
                bool y_odd = detail::is_odd_integer(y);
                // Distinguish +0 from -0 to honor the signed-zero base rule.
                bool x_is_neg = std::bit_cast<std::uint64_t>(x) >> 63;
                if (y < 0.0) {
                        return (y_odd && x_is_neg) ? -pinf : pinf;
                }
                // y > 0: pow(+/-0, y>0) is +/-0 if y is odd integer, else +0
                return (y_odd && x_is_neg) ? -0.0 : 0.0;
        }

        // x < 0: integer exponent uses |x|^y with sign; non-integer is NaN.
        // Use the squaring fast path for exact integer-power semantics when y fits
        // safely in long long; fall back to the transcendental path only when y is
        // outside that range (where y is necessarily an even integer above 2^53,
        // so the sign flip via is_odd_integer correctly returns false).
        if (x < 0.0) {
                if (!detail::is_integer(y)) return std::numeric_limits<double>::quiet_NaN();
                bool y_odd = detail::is_odd_integer(y);
                if (y > -detail::LL_BOUND_DOUBLE && y < detail::LL_BOUND_DOUBLE) {
                        long long n = static_cast<long long>(y);
                        unsigned long long m = (n >= 0) ? static_cast<unsigned long long>(n)
                                                        : static_cast<unsigned long long>(-n);
                        double mag = detail::pow_by_squaring(-x, m);
                        if (n < 0) mag = 1.0 / mag;
                        return y_odd ? -mag : mag;
                }
                double mag = exp2(y * log2(-x));
                return y_odd ? -mag : mag;
        }

        // x > 0, finite y: integer fast path when applicable, else exp2(y*log2(x)).
        // Bounds are exclusive on both sides so the cast and the subsequent unary
        // minus are always defined: 2^63 is representable as double but outside
        // long long; -2^63 cast then negated overflows.
        if (detail::is_integer(y)) {
                if (y > -detail::LL_BOUND_DOUBLE && y < detail::LL_BOUND_DOUBLE) {
                        long long n = static_cast<long long>(y);
                        if (n >= 0) return detail::pow_by_squaring(x, static_cast<unsigned long long>(n));
                        unsigned long long m = static_cast<unsigned long long>(-n);
                        return 1.0 / detail::pow_by_squaring(x, m);
                }
        }
        return exp2(y * log2(x));
}

constexpr float pow(float x, float y) {
        if (y == 0.0f) return 1.0f;
        if (x == 1.0f) return 1.0f;
        if (x != x || y != y) return std::numeric_limits<float>::quiet_NaN();

        const float pinf = std::numeric_limits<float>::infinity();

        if (x == -1.0f && (y == pinf || y == -pinf)) return 1.0f;

        if (y == pinf) {
                float ax = (x < 0.0f) ? -x : x;
                return (ax < 1.0f) ? 0.0f : pinf;
        }
        if (y == -pinf) {
                float ax = (x < 0.0f) ? -x : x;
                return (ax < 1.0f) ? pinf : 0.0f;
        }

        if (x == pinf) {
                return (y < 0.0f) ? 0.0f : pinf;
        }
        if (x == -pinf) {
                bool y_odd = detail::is_odd_integer(y);
                if (y < 0.0f) return y_odd ? -0.0f : 0.0f;
                return y_odd ? -pinf : pinf;
        }

        if (x == 0.0f) {
                bool y_odd = detail::is_odd_integer(y);
                bool x_is_neg = std::bit_cast<std::uint32_t>(x) >> 31;
                if (y < 0.0f) {
                        return (y_odd && x_is_neg) ? -pinf : pinf;
                }
                return (y_odd && x_is_neg) ? -0.0f : 0.0f;
        }

        if (x < 0.0f) {
                if (!detail::is_integer(y)) return std::numeric_limits<float>::quiet_NaN();
                bool y_odd = detail::is_odd_integer(y);
                if (y > -detail::LL_BOUND_FLOAT && y < detail::LL_BOUND_FLOAT) {
                        long long n = static_cast<long long>(y);
                        unsigned long long m = (n >= 0) ? static_cast<unsigned long long>(n)
                                                        : static_cast<unsigned long long>(-n);
                        float mag = detail::pow_by_squaring(-x, m);
                        if (n < 0) mag = 1.0f / mag;
                        return y_odd ? -mag : mag;
                }
                float mag = exp2(y * log2(-x));
                return y_odd ? -mag : mag;
        }

        if (detail::is_integer(y)) {
                if (y > -detail::LL_BOUND_FLOAT && y < detail::LL_BOUND_FLOAT) {
                        long long n = static_cast<long long>(y);
                        if (n >= 0) return detail::pow_by_squaring(x, static_cast<unsigned long long>(n));
                        unsigned long long m = static_cast<unsigned long long>(-n);
                        return 1.0f / detail::pow_by_squaring(x, m);
                }
        }
        return exp2(y * log2(x));
}

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

1	#pragma once
2	// pow.hpp: constexpr power function 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	// Two overload families
11	// -----------------------------------------------------------------------------
12	//
13	// 1. Integer-exponent fast path: pow(T, int)
14	// - Fast exponentiation by squaring; pure integer logic.
15	// - O(log \|n\|) multiplications, no transcendentals.
16	// - Selected when the user knows the exponent is an integer.
17	//
18	// 2. General overload: pow(T, T)
19	// - First tries the integer fast path when the floating-point exponent
20	// is exactly an integer (typical for compile-time constants).
21	// - Otherwise computes pow(x, y) = exp2(y * log2(x)) for x > 0.
22	// - Negative base with integer exponent: pow(\|x\|, y) with sign from y.
23	// - Negative base with non-integer exponent: NaN.
24	// - Full IEEE-754 special-case table per C99 7.12.7.4.
25	//
26	// -----------------------------------------------------------------------------
27	// IEEE-754 special-value handling (C99 7.12.7.4)
28	// -----------------------------------------------------------------------------
29	// pow(x, 0) -> 1 (any x including NaN, inf)
30	// pow(1, y) -> 1 (any y including NaN, inf)
31	// pow(NaN, y) -> NaN (y != 0)
32	// pow(x, NaN) -> NaN (x != 1)
33	// pow(+/-0, y<0, odd int)-> +/-inf (signed)
34	// pow(+/-0, y<0, else) -> +inf
35	// pow(+/-0, y>0, odd int)-> +/-0 (signed)
36	// pow(+/-0, y>0, else) -> +0
37	// pow(-1, +/-inf) -> 1
38	// pow(\|x\|<1,+inf) -> +0
39	// pow(\|x\|<1,-inf) -> +inf
40	// pow(\|x\|>1,+inf) -> +inf
41	// pow(\|x\|>1,-inf) -> +0
42	// pow(+inf, y<0) -> +0
43	// pow(+inf, y>0) -> +inf
44	// pow(-inf, y<0, odd int)-> -0
45	// pow(-inf, y<0, else) -> +0
46	// pow(-inf, y>0, odd int)-> -inf
47	// pow(-inf, y>0, else) -> +inf
48	// pow(x<0, finite non-integer y) -> NaN
49
50	#include <limits>
51
52	#include <math/constexpr_math/detail.hpp>
53	#include <math/constexpr_math/exp2.hpp>
54	#include <math/constexpr_math/log2.hpp>
55
56	namespace sw { namespace math { namespace constexpr_math {
57
58	// -----------------------------------------------------------------------------
59	// Integer-exponent fast path
60	// -----------------------------------------------------------------------------
61
62	constexpr double pow(double x, int n) {
63	if (n == 0) return 1.0;
64	if (n < 0) {
65	// Promote to long long to safely negate INT_MIN.
66	unsigned long long exp = static_cast<unsigned long long>(-static_cast<long long>(n));
67	return 1.0 / detail::pow_by_squaring(x, exp);
68	}
69	return detail::pow_by_squaring(x, static_cast<unsigned long long>(n));
70	}
71
72	constexpr float pow(float x, int n) {
73	if (n == 0) return 1.0f;
74	if (n < 0) {
75	unsigned long long exp = static_cast<unsigned long long>(-static_cast<long long>(n));
76	return 1.0f / detail::pow_by_squaring(x, exp);
77	}
78	return detail::pow_by_squaring(x, static_cast<unsigned long long>(n));
79	}
80
81	// -----------------------------------------------------------------------------
82	// General overload
83	// -----------------------------------------------------------------------------
84
85	constexpr double pow(double x, double y) {	21✔
86	// Cases that win regardless of NaN: pow(x, 0) and pow(1, y).
87	if (y == 0.0) return 1.0;	21✔
88	if (x == 1.0) return 1.0;	21✔
89
90	// NaN propagation (after the two NaN-overriding cases above).
91	if (x != x \|\| y != y) return std::numeric_limits<double>::quiet_NaN();	21✔
92
93	const double pinf = std::numeric_limits<double>::infinity();	21✔
94
95	// Special case: pow(-1, +/-inf) == 1
96	if (x == -1.0 && (y == pinf \|\| y == -pinf)) return 1.0;	21✔
97
98	// y is +/- infinity
99	if (y == pinf) {	21✔
NEW 100	double ax = (x < 0.0) ? -x : x;	×
NEW 101	if (ax < 1.0) return 0.0;	×
102	// \|x\| > 1 (==1 was handled by pow(1,y) early-out and x==-1 above)
NEW 103	return pinf;	×
104	}
105	if (y == -pinf) {	21✔
NEW 106	double ax = (x < 0.0) ? -x : x;	×
NEW 107	if (ax < 1.0) return pinf;	×
NEW 108	return 0.0;	×
109	}
110
111	// x is +/- infinity (y is finite, non-zero)
112	if (x == pinf) {	21✔
NEW 113	return (y < 0.0) ? 0.0 : pinf;	×
114	}
115	if (x == -pinf) {	21✔
NEW 116	bool y_odd = detail::is_odd_integer(y);	×
NEW 117	if (y < 0.0) return y_odd ? -0.0 : 0.0;	×
NEW 118	return y_odd ? -pinf : pinf;	×
119	}
120
121	// x is +/- 0 (y is finite, non-zero, not infinity)
122	if (x == 0.0) {	21✔
NEW 123	bool y_odd = detail::is_odd_integer(y);	×
124	// Distinguish +0 from -0 to honor the signed-zero base rule.
NEW 125	bool x_is_neg = std::bit_cast<std::uint64_t>(x) >> 63;	×
NEW 126	if (y < 0.0) {	×
NEW 127	return (y_odd && x_is_neg) ? -pinf : pinf;	×
128	}
129	// y > 0: pow(+/-0, y>0) is +/-0 if y is odd integer, else +0
NEW 130	return (y_odd && x_is_neg) ? -0.0 : 0.0;	×
131	}
132
133	// x < 0: integer exponent uses \|x\|^y with sign; non-integer is NaN.
134	// Use the squaring fast path for exact integer-power semantics when y fits
135	// safely in long long; fall back to the transcendental path only when y is
136	// outside that range (where y is necessarily an even integer above 2^53,
137	// so the sign flip via is_odd_integer correctly returns false).
138	if (x < 0.0) {	21✔
139	if (!detail::is_integer(y)) return std::numeric_limits<double>::quiet_NaN();	4✔
140	bool y_odd = detail::is_odd_integer(y);	4✔
141	if (y > -detail::LL_BOUND_DOUBLE && y < detail::LL_BOUND_DOUBLE) {	4✔
142	long long n = static_cast<long long>(y);	4✔
143	unsigned long long m = (n >= 0) ? static_cast<unsigned long long>(n)	4✔
144	: static_cast<unsigned long long>(-n);
145	double mag = detail::pow_by_squaring(-x, m);	4✔
146	if (n < 0) mag = 1.0 / mag;	4✔
147	return y_odd ? -mag : mag;	4✔
148	}
NEW 149	double mag = exp2(y * log2(-x));	×
NEW 150	return y_odd ? -mag : mag;	×
151	}
152
153	// x > 0, finite y: integer fast path when applicable, else exp2(y*log2(x)).
154	// Bounds are exclusive on both sides so the cast and the subsequent unary
155	// minus are always defined: 2^63 is representable as double but outside
156	// long long; -2^63 cast then negated overflows.
157	if (detail::is_integer(y)) {	17✔
158	if (y > -detail::LL_BOUND_DOUBLE && y < detail::LL_BOUND_DOUBLE) {	11✔
159	long long n = static_cast<long long>(y);	11✔
160	if (n >= 0) return detail::pow_by_squaring(x, static_cast<unsigned long long>(n));	11✔
161	unsigned long long m = static_cast<unsigned long long>(-n);	2✔
162	return 1.0 / detail::pow_by_squaring(x, m);	2✔
163	}
164	}
165	return exp2(y * log2(x));	6✔
166	}
167
168	constexpr float pow(float x, float y) {	7✔
169	if (y == 0.0f) return 1.0f;	7✔
170	if (x == 1.0f) return 1.0f;	7✔
171	if (x != x \|\| y != y) return std::numeric_limits<float>::quiet_NaN();	7✔
172
173	const float pinf = std::numeric_limits<float>::infinity();	7✔
174
175	if (x == -1.0f && (y == pinf \|\| y == -pinf)) return 1.0f;	7✔
176
177	if (y == pinf) {	7✔
NEW 178	float ax = (x < 0.0f) ? -x : x;	×
NEW 179	return (ax < 1.0f) ? 0.0f : pinf;	×
180	}
181	if (y == -pinf) {	7✔
NEW 182	float ax = (x < 0.0f) ? -x : x;	×
NEW 183	return (ax < 1.0f) ? pinf : 0.0f;	×
184	}
185
186	if (x == pinf) {	7✔
NEW 187	return (y < 0.0f) ? 0.0f : pinf;	×
188	}
189	if (x == -pinf) {	7✔
NEW 190	bool y_odd = detail::is_odd_integer(y);	×
NEW 191	if (y < 0.0f) return y_odd ? -0.0f : 0.0f;	×
NEW 192	return y_odd ? -pinf : pinf;	×
193	}
194
195	if (x == 0.0f) {	7✔
NEW 196	bool y_odd = detail::is_odd_integer(y);	×
NEW 197	bool x_is_neg = std::bit_cast<std::uint32_t>(x) >> 31;	×
NEW 198	if (y < 0.0f) {	×
NEW 199	return (y_odd && x_is_neg) ? -pinf : pinf;	×
200	}
NEW 201	return (y_odd && x_is_neg) ? -0.0f : 0.0f;	×
202	}
203
204	if (x < 0.0f) {	7✔
205	if (!detail::is_integer(y)) return std::numeric_limits<float>::quiet_NaN();	2✔
206	bool y_odd = detail::is_odd_integer(y);	2✔
207	if (y > -detail::LL_BOUND_FLOAT && y < detail::LL_BOUND_FLOAT) {	2✔
208	long long n = static_cast<long long>(y);	2✔
209	unsigned long long m = (n >= 0) ? static_cast<unsigned long long>(n)	2✔
210	: static_cast<unsigned long long>(-n);
211	float mag = detail::pow_by_squaring(-x, m);	2✔
212	if (n < 0) mag = 1.0f / mag;	2✔
213	return y_odd ? -mag : mag;	2✔
214	}
NEW 215	float mag = exp2(y * log2(-x));	×
NEW 216	return y_odd ? -mag : mag;	×
217	}
218
219	if (detail::is_integer(y)) {	5✔
220	if (y > -detail::LL_BOUND_FLOAT && y < detail::LL_BOUND_FLOAT) {	3✔
221	long long n = static_cast<long long>(y);	3✔
222	if (n >= 0) return detail::pow_by_squaring(x, static_cast<unsigned long long>(n));	3✔
223	unsigned long long m = static_cast<unsigned long long>(-n);	1✔
224	return 1.0f / detail::pow_by_squaring(x, m);	1✔
225	}
226	}
227	return exp2(y * log2(x));	2✔
228	}
229
230	}}} // namespace sw::math::constexpr_math

stillwater-sc / universal / 24955899585

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