24944226215

Committed 26 Apr 2026 12:27AM UTC coverage: 84.328% (+0.01%) from 84.317%

Build # 24944226215

Build Type

push

github

Committed by

web-flow

Commit Message

feat(math): add sw::math::constexpr_math facility with constexpr log2 (#762)

* feat(math): add sw::math::constexpr_math facility with constexpr log2

Introduces a new self-contained library of constexpr math functions for
IEEE-754 float and double, parallel to the existing sw::math::function
facility but specifically for compile-time evaluation. Per the issue
brief: "constexpr is a slightly different use case than templated math
functions ... it would be nice to have this as a separate and
specialized include set ... independently useful outside of Universal's
number systems".

Layout:
  include/sw/math/constexpr_math.hpp       -- aggregator
  include/sw/math/constexpr_math/log2.hpp  -- base-2 logarithm
  internal/constexpr_math/api/log2.cpp     -- regression test (cm_log2)

Namespace: sw::math::constexpr_math

Algorithm (log2):
  1. Decompose x = 2^E * M with M in [1, 2) via std::bit_cast<uint64_t>.
  2. Range-reduce to M in [1/sqrt(2), sqrt(2)] by halving M and
     incrementing E when M > sqrt(2).
  3. With u = (M-1)/(M+1), |u| <= 0.1716, evaluate
       ln(M) = 2 * artanh(u) = 2 * (u + u^3/3 + u^5/5 + ...)
     via Horner. Coefficients are odd reciprocals (1/3, 1/5, ...) --
     derivable from first principles, no external tooling
     (no Sollya, no Remez).
  4. log2(x) = E + ln(M) / ln(2)

Polynomial degree:
  - double: degree 21 (terms up to u^21/21) -> ~5e-17 truncation error
  - float : degree 11 (terms up to u^11/11) -> ~1e-9 truncation error

Special-value handling follows IEEE-754:
  log2(NaN)  -> NaN
  log2(x<0)  -> NaN
  log2(0)    -> -infinity
  log2(+inf) -> +infinity

Subnormals are normalized via a software shift loop (constexpr-safe).

Verification:
  - 22 static_assert checks for exact powers of 2, special values,
    log2(3.14) and log2(0.75) acceptance forms
  - Runtime sweep over 23 hand-picked points + subnormal sample
  - Out-of-band stress sweep (100K random samples in [1e-100, 1e100]):
    max relative error 1.74e-16 = 0.78... (continued)

Coverage Stats

90 of 104 new or added lines in 2 files covered. (86.54%)

4 existing lines in 3 files now uncovered.

45103 of 53485 relevant lines covered (84.33%)

6403422.69 hits per line

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

75.76

/internal/constexpr_math/api/log2.cpp

// log2.cpp: regression for sw::math::constexpr_math::log2(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 <cmath>
#include <iostream>
#include <iomanip>
#include <limits>

#include <math/constexpr_math.hpp>

namespace cm = sw::math::constexpr_math;

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

// Exact powers of 2 -- the integer-exponent path, must be bit-exact.
static_assert(cm::log2(1.0)    == 0.0, "log2(1) == 0");
static_assert(cm::log2(2.0)    == 1.0, "log2(2) == 1");
static_assert(cm::log2(4.0)    == 2.0, "log2(4) == 2");
static_assert(cm::log2(8.0)    == 3.0, "log2(8) == 3");
static_assert(cm::log2(0.5)    == -1.0, "log2(0.5) == -1");
static_assert(cm::log2(0.25)   == -2.0, "log2(0.25) == -2");
static_assert(cm::log2(1024.0) == 10.0, "log2(1024) == 10");

// Float overload, exact powers of 2
static_assert(cm::log2(1.0f) == 0.0f, "log2f(1) == 0");
static_assert(cm::log2(8.0f) == 3.0f, "log2f(8) == 3");
static_assert(cm::log2(0.5f) == -1.0f, "log2f(0.5) == -1");

// Special values
static_assert(cm::log2(0.0)  == -std::numeric_limits<double>::infinity(),
              "log2(0) == -inf");
static_assert(cm::log2(0.0f) == -std::numeric_limits<float>::infinity(),
              "log2f(0) == -inf");
static_assert(cm::log2(std::numeric_limits<double>::infinity())
              == std::numeric_limits<double>::infinity(),
              "log2(+inf) == +inf");

// log2(x) for x < 0 must yield NaN. NaN != NaN, so the test is `result != result`.
static_assert(cm::log2(-1.0) != cm::log2(-1.0), "log2(-1) == NaN");
static_assert(cm::log2(-2.5f) != cm::log2(-2.5f), "log2f(-2.5) == NaN");

// Acceptance form for issue #423: enables `constexpr takum16 t = 3.14f;`
// log2(3.14) = ln(3.14) / ln(2) ~= 1.65076455911...
constexpr double cx_log2_3p14 = cm::log2(3.14);
static_assert(cx_log2_3p14 > 1.65076 && cx_log2_3p14 < 1.65077,
              "log2(3.14) ~= 1.65076455...");

// Acceptance: enables compile-time lns encoding from a float literal.
// For lns<8,4>, log2(0.75) is needed in the value->encoding mapping.
constexpr double cx_log2_0p75 = cm::log2(0.75);
static_assert(cx_log2_0p75 > -0.4151 && cx_log2_0p75 < -0.4149,
              "log2(0.75) ~= -0.41504");

// ----------------------------------------------------------------------------
// Runtime cross-check vs std::log2 over a representative grid
// ----------------------------------------------------------------------------

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

        int errors = 0;
        auto check = [&](const char* name, 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
                                  << "  our=" << std::setprecision(17) << our
                                  << "  ref=" << ref
                                  << "  rel-err=" << err << '\n';
                }
        };

        // Sweep across a wide dynamic range and a fine-grained set inside [1, 2).
        const double points[] = {
                0.5, 0.75, 1.0, 1.25, 1.5, 1.7320508, 1.9999999,
                2.0, 3.0, 3.14, 5.0, 7.0, 10.0, 100.0, 1024.0,
                1e-10, 1e-5, 1e-3, 1e3, 1e10, 1e100, 1e-300, 1e300,
        };
        for (double x : points) {
                double our = cm::log2(x);
                double ref = std::log2(x);
                // Allow up to 4 ulp (double eps ~2.2e-16, so ~1e-15)
                check("double sweep", our, ref, 1e-15);
        }

        // Subnormal sample
        {
                double sub = std::numeric_limits<double>::min() / 4.0;  // subnormal
                double our = cm::log2(sub);
                double ref = std::log2(sub);
                check("double subnormal", our, ref, 1e-14);
        }

        // Float sweep
        const float fpoints[] = {
                0.5f, 0.75f, 1.0f, 1.5f, 2.0f, 3.14f, 1024.0f, 1e10f, 1e-10f,
        };
        for (float x : fpoints) {
                float our = cm::log2(x);
                float ref = std::log2(x);
                double err = (ref == 0.0f) ? std::abs(our) : std::abs((our - ref) / ref);
                // Float epsilon ~1.19e-7; allow 1e-6 relative.
                if (err > 1e-6) {
                        ++errors;
                        std::cout << "FAIL float sweep  x=" << x
                                  << "  our=" << our << "  ref=" << ref
                                  << "  rel-err=" << err << '\n';
                }
        }

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

1	// log2.cpp: regression for sw::math::constexpr_math::log2(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 <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	// Exact powers of 2 -- the integer-exponent path, must be bit-exact.
22	static_assert(cm::log2(1.0) == 0.0, "log2(1) == 0");
23	static_assert(cm::log2(2.0) == 1.0, "log2(2) == 1");
24	static_assert(cm::log2(4.0) == 2.0, "log2(4) == 2");
25	static_assert(cm::log2(8.0) == 3.0, "log2(8) == 3");
26	static_assert(cm::log2(0.5) == -1.0, "log2(0.5) == -1");
27	static_assert(cm::log2(0.25) == -2.0, "log2(0.25) == -2");
28	static_assert(cm::log2(1024.0) == 10.0, "log2(1024) == 10");
29
30	// Float overload, exact powers of 2
31	static_assert(cm::log2(1.0f) == 0.0f, "log2f(1) == 0");
32	static_assert(cm::log2(8.0f) == 3.0f, "log2f(8) == 3");
33	static_assert(cm::log2(0.5f) == -1.0f, "log2f(0.5) == -1");
34
35	// Special values
36	static_assert(cm::log2(0.0) == -std::numeric_limits<double>::infinity(),
37	"log2(0) == -inf");
38	static_assert(cm::log2(0.0f) == -std::numeric_limits<float>::infinity(),
39	"log2f(0) == -inf");
40	static_assert(cm::log2(std::numeric_limits<double>::infinity())
41	== std::numeric_limits<double>::infinity(),
42	"log2(+inf) == +inf");
43
44	// log2(x) for x < 0 must yield NaN. NaN != NaN, so the test is `result != result`.
45	static_assert(cm::log2(-1.0) != cm::log2(-1.0), "log2(-1) == NaN");
46	static_assert(cm::log2(-2.5f) != cm::log2(-2.5f), "log2f(-2.5) == NaN");
47
48	// Acceptance form for issue #423: enables `constexpr takum16 t = 3.14f;`
49	// log2(3.14) = ln(3.14) / ln(2) ~= 1.65076455911...
50	constexpr double cx_log2_3p14 = cm::log2(3.14);
51	static_assert(cx_log2_3p14 > 1.65076 && cx_log2_3p14 < 1.65077,
52	"log2(3.14) ~= 1.65076455...");
53
54	// Acceptance: enables compile-time lns encoding from a float literal.
55	// For lns<8,4>, log2(0.75) is needed in the value->encoding mapping.
56	constexpr double cx_log2_0p75 = cm::log2(0.75);
57	static_assert(cx_log2_0p75 > -0.4151 && cx_log2_0p75 < -0.4149,
58	"log2(0.75) ~= -0.41504");
59
60	// ----------------------------------------------------------------------------
61	// Runtime cross-check vs std::log2 over a representative grid
62	// ----------------------------------------------------------------------------
63
64	int main() {	1✔
65	std::cout << "sw::math::constexpr_math::log2 verification\n";	1✔
66
67	int errors = 0;	1✔
68	auto check = [&](const char* name, double our, double ref, double tol) {	24✔
69	double err = (ref == 0.0) ? std::abs(our) : std::abs((our - ref) / ref);	24✔
70	if (err > tol) {	24✔
NEW 71	++errors;	×
72	std::cout << "FAIL " << name
NEW 73	<< " our=" << std::setprecision(17) << our	×
NEW 74	<< " ref=" << ref	×
NEW 75	<< " rel-err=" << err << '\n';	×
76	}
77	};	25✔
78
79	// Sweep across a wide dynamic range and a fine-grained set inside [1, 2).
80	const double points[] = {	1✔
81	0.5, 0.75, 1.0, 1.25, 1.5, 1.7320508, 1.9999999,
82	2.0, 3.0, 3.14, 5.0, 7.0, 10.0, 100.0, 1024.0,
83	1e-10, 1e-5, 1e-3, 1e3, 1e10, 1e100, 1e-300, 1e300,
84	};
85	for (double x : points) {	24✔
86	double our = cm::log2(x);	23✔
87	double ref = std::log2(x);	23✔
88	// Allow up to 4 ulp (double eps ~2.2e-16, so ~1e-15)
89	check("double sweep", our, ref, 1e-15);	23✔
90	}
91
92	// Subnormal sample
93	{
94	double sub = std::numeric_limits<double>::min() / 4.0; // subnormal	1✔
95	double our = cm::log2(sub);	1✔
96	double ref = std::log2(sub);	1✔
97	check("double subnormal", our, ref, 1e-14);	1✔
98	}
99
100	// Float sweep
101	const float fpoints[] = {	1✔
102	0.5f, 0.75f, 1.0f, 1.5f, 2.0f, 3.14f, 1024.0f, 1e10f, 1e-10f,
103	};
104	for (float x : fpoints) {	10✔
105	float our = cm::log2(x);	9✔
106	float ref = std::log2(x);	9✔
107	double err = (ref == 0.0f) ? std::abs(our) : std::abs((our - ref) / ref);	9✔
108	// Float epsilon ~1.19e-7; allow 1e-6 relative.
109	if (err > 1e-6) {	9✔
NEW 110	++errors;	×
NEW 111	std::cout << "FAIL float sweep x=" << x	×
NEW 112	<< " our=" << our << " ref=" << ref	×
NEW 113	<< " rel-err=" << err << '\n';	×
114	}
115	}
116
117	std::cout << "constexpr_math::log2: " << (errors == 0 ? "PASS" : "FAIL")	1✔
118	<< " (" << errors << " error" << (errors == 1 ? "" : "s") << ")\n";	1✔
119	return (errors == 0) ? 0 : 1;	1✔
120	}

stillwater-sc / universal / 24944226215

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