Detailed Implementation of Logarithmic Functions in the C Standard Library: log2f, log10f, and logf

Table of Contents

OverviewMathematical Background and Algorithm PrinciplesDetailed Implementation of log2f

• • Algorithm Principles and Code Explanation
• • Numerical Stability Analysis
• • Handling of Subnormal NumbersDetailed Implementation of log10f
• • Algorithm Principles and Code Explanation
• • Numerical Stability Analysis
• • Handling of Subnormal NumbersDetailed Implementation of logf
• • Algorithm Principles and Code Explanation
• • Numerical Stability Analysis
• • Handling of Subnormal NumbersPerformance Optimization StrategiesCompliance with StandardsConclusion

Overview

Logarithmic functions have wide applications in scientific computing, signal processing, data compression, and other fields. The musl libc provides a complete implementation of the logarithmic function family, including <span>log2f</span> (base-2 logarithm), <span>log10f</span> (base-10 logarithm), and <span>logf</span> (natural logarithm). These functions are designed with careful algorithms and numerical optimizations to ensure excellent numerical precision while maintaining performance.

Mathematical Background and Algorithm Principles

Basic Properties of Logarithmic Functions

Logarithmic functions have the following important properties:

• •
• •
• •
• • (Change of Base Formula)

Algorithm Strategy

The implementation of logarithmic functions in musl adopts a unified strategy:

1. 1. Parameter Normalization: Normalize the input to a specific range
2. 2. Table Lookup Approximation: Use precomputed lookup tables to obtain approximate values
3. 3. Polynomial Correction: Use low-order polynomials to correct residuals
4. 4. Result Reconstruction: Combine parts to obtain the final value

Detailed Implementation of log2f

Algorithm Principles and Code Explanation

<span>log2f</span> calculates using a combination of table lookup and polynomial approximation methods.

Data Structures and Constants

#define LOG2F_TABLE_BITS 4
#define LOG2F_POLY_ORDER 4

const struct log2f_data __log2f_data = {
  .tab = {
  { 0x1.661ec79f8f3bep+0, -0x1.efec65b963019p-2 },
  // ... 15 entries
  { 0x1.767dcf5534862p-1, 0x1.ce0a44eb17bccp-2 },
  },
  .poly = {
  -0x1.712b6f70a7e4dp-2, 0x1.ecabf496832ep-2, -0x1.715479ffae3dep-1,
  0x1.715475f35c8b8p0,
  }
};

Lookup Table Design:

• • 16 entries ( )
• • Each entry contains <span>invc</span> (reciprocal value) and <span>logc</span> (logarithmic value)
• • The entries cover the range to ensure that normalized inputs fall within this range

Input Processing and Special Value Checks

ix = asuint(x);
if (WANT_ROUNDING && predict_false(ix == 0x3f800000))
    return 0;
if (predict_false(ix - 0x00800000 >= 0x7f800000 - 0x00800000)) {
    if (ix * 2 == 0)
        return __math_divzerof(1);
    if (ix == 0x7f800000)
        return x;
    if ((ix & 0x80000000) || ix * 2 >= 0xff000000)
        return __math_invalidf(x);
    ix = asuint(x * 0x1p23f);
    ix -= 23 << 23;
}

Handling of Special Cases:

• • :
• • : return (division by zero)
• • : return
• • or NaN: return NaN
• • Subnormal numbers: normalize by scaling

Parameter Normalization and Range Reduction

tmp = ix - OFF;
i = (tmp >> (23 - LOG2F_TABLE_BITS)) % N;
top = tmp & 0xff800000;
iz = ix - top;
k = (int32_t)tmp >> 23;
invc = T[i].invc;
logc = T[i].logc;
z = (double_t)asfloat(iz);

Mathematical Principles:

Decompose as:

where , .

Then:

Further decompose as:

where is the center value from the lookup table, and is a small residual.

Polynomial Approximation

r = z * invc - 1;
y0 = logc + (double_t)k;

r2 = r * r;
y = A[1] * r + A[2];
y = A[0] * r2 + y;
p = A[3] * r + y0;
y = y * r2 + p;

Approximation Formula:

where is approximated using a polynomial.

Numerical Stability Analysis

Error Sources:

1. 1. Entry Error: Precomputed has finite precision
2. 2. Polynomial Approximation Error: Maximum error of the 4th degree polynomial within the target range
3. 3. Rounding Error in Calculations: Rounding of intermediate double precision calculations

Error Control:

• • Use double precision for intermediate calculations
• • Polynomial coefficients optimized for the target range
• • Entry design ensures that is very small, improving polynomial approximation accuracy

Handling of Subnormal Numbers

Normalization of Subnormal Numbers:

ix = asuint(x * 0x1p23f);
ix -= 23 << 23;

Subnormal numbers are converted to normal numbers by multiplying by while adjusting the exponent for compensation.

Detailed Implementation of log10f

Algorithm Principles and Code Explanation

<span>log10f</span> calculates based on the change of base formula, but optimizes performance using a direct approximation method.

Constant Definitions

static const float
ivln10hi  =  4.3432617188e-01, /* 0x3ede6000 */
ivln10lo  = -3.1689971365e-05, /* 0xb804ead9 */
log10_2hi =  3.0102920532e-01, /* 0x3e9a2080 */
log10_2lo =  7.9034151668e-07, /* 0x355427db */
Lg1 = 0xaaaaaa.0p-24, /* 0.66666662693 */
Lg2 = 0xccce13.0p-25, /* 0.40000972152 */
Lg3 = 0x91e9ee.0p-25, /* 0.28498786688 */
Lg4 = 0xf89e26.0p-26; /* 0.24279078841 */

Constant Design:

• • <span>ivln10hi</span> and <span>ivln10lo</span>: High and low parts decomposition of
• • <span>log10_2hi</span> and <span>log10_2lo</span>: High and low parts decomposition of
• • <span>Lg1</span>–<span>Lg4</span>: Coefficients of the polynomial approximation

Input Processing and Range Reduction

ix = u.i;
k = 0;
if (ix < 0x00800000 || ix>>31) {
    if (ix<<1 == 0)
        return -1/(x*x);
    if (ix>>31)
        return (x-x)/0.0f;
    k -= 25;
    x *= 0x1p25f;
    u.f = x;
    ix = u.i;
}

Special Cases:

• • : return
• • : return NaN
• • Subnormal numbers: normalize by scaling

Range Reduction to

ix += 0x3f800000 - 0x3f3504f3;
k += (int)(ix>>23) - 0x7f;
ix = (ix & 0x007fffff) + 0x3f3504f3;
u.i = ix;
x = u.f;

Normalize to the range while recording the exponent adjustment .

Rational Function Approximation

f = x - 1.0f;
s = f/(2.0f + f);
z = s*s;
w = z*z;
t1= w*(Lg2+w*Lg4);
t2= z*(Lg1+w*Lg3);
R = t2 + t1;
hfsq = 0.5f*f*f;

Mathematical Principles:Using the transformation:

Then:

Actual use of optimized rational function approximation.

High Precision Calculation

hi = f - hfsq;
u.f = hi;
u.i &= 0xfffff000;
hi = u.f;
lo = f - hi - hfsq + s*(hfsq+R);

Using high/low separation techniques to reduce rounding errors.

Final Result Combination

dk = k;
return dk*log10_2lo + (lo+hi)*ivln10lo + lo*ivln10hi + hi*ivln10hi + dk*log10_2hi;

Reconstruction Formula:

where is calculated through .

Numerical Stability Analysis

Precision Assurance:

• • Use high/low separation to reduce rounding errors
• • Rational function approximation error is less than
• • High and low parts decomposition of constants ensures accurate multiplication

Handling of Subnormal Numbers

Normalize subnormal numbers by multiplying by while adjusting the exponent :

k -= 25;
x *= 0x1p25f;

Detailed Implementation of logf

Algorithm Principles and Code Explanation

<span>logf</span> calculates using a similar table lookup + polynomial approximation method as <span>log2f</span>.

Data Structures and Constants

#define LOGF_TABLE_BITS 4
#define LOGF_POLY_ORDER 4

const struct logf_data __logf_data = {
  .tab = {
  { 0x1.661ec79f8f3bep+0, -0x1.57bf7808caadep-2 },
  // ... 15 entries
  { 0x1.767dcf5534862p-1, 0x1.4043057b6ee09p-2 },
  },
  .ln2 = 0x1.62e42fefa39efp-1,
  .poly = {
  -0x1.00ea348b88334p-2, 0x1.5575b0be00b6ap-2, -0x1.ffffef20a4123p-2,
  }
};

Lookup Table Design:

• • 16 entries, covering the range
• • Contains high precision values of
• • 3rd degree polynomial (4th order, but the leading coefficient is 1)

Input Processing and Parameter Normalization

ix = asuint(x);
if (WANT_ROUNDING && predict_false(ix == 0x3f800000))
    return 0;
if (predict_false(ix - 0x00800000 >= 0x7f800000 - 0x00800000)) {
    // Special case handling, same as log2f
}

Special case handling is identical to that of <span>log2f</span>.

Range Reduction

tmp = ix - OFF;
i = (tmp >> (23 - LOGF_TABLE_BITS)) % N;
k = (int32_t)tmp >> 23;
iz = ix - (tmp & 0xff800000);
invc = T[i].invc;
logc = T[i].logc;
z = (double_t)asfloat(iz);

Decompose as:

where .

Polynomial Approximation

r = z * invc - 1;
y0 = logc + (double_t)k * Ln2;

r2 = r * r;
y = A[1] * r + A[2];
y = A[0] * r2 + y;
y = y * r2 + (y0 + r);

Approximation Formula:

where is approximated using a 3rd degree polynomial.

Numerical Stability Analysis

Error Control:

• • Entries and polynomial coefficients optimized for the target range
• • Use double precision for intermediate calculations
• • Polynomial design ensures that the function value and first derivative match at

Precision Assurance: Error is less than 2 ULP across the entire domain.

Handling of Subnormal Numbers

Normalization method for subnormal numbers is the same as that of <span>log2f</span>:

ix = asuint(x * 0x1p23f);
ix -= 23 << 23;

Performance Optimization Strategies

Algorithm-Level Optimizations

1. 1. Unified Infrastructure: The three logarithmic functions share a similar table lookup and polynomial approximation framework
2. 2. Range Reduction: Convert arbitrary inputs to a small range to improve polynomial approximation accuracy
3. 3. Table Optimization: 16-entry tables balance precision and cache efficiency

Implementation-Level Optimizations

1. 1. Bit Manipulation: Directly manipulate floating-point bit patterns to quickly extract exponents and mantissas
2. 2. Double Precision Intermediate Calculations: Use double precision for critical calculations to reduce rounding errors
3. 3. Pipelined Polynomial Calculations: Optimize the order of polynomial evaluations to reduce instruction dependencies

Fast Path for Special Values

• • Provide quick returns for cases
• • Early checks and handling for special values (0, infinity, NaN)
• • Use <span>predict_false</span> to hint the compiler to optimize rare paths

Compliance with Standards

Compliance with IEEE 754 Standard

1. 1. Handling of Special Values:

• •
• •
• •
• • (when )
• •

2. Precision Requirements: Error is less than 2 ULP, meeting standard requirements

3. Rounding Modes: Supports all IEEE 754 rounding modes

Preservation of Mathematical Properties

• • Monotonically increasing function
• • Satisfies the basic algebraic properties of logarithmic functions
• • Maintains continuity within the domain

Conclusion

The implementation of the logarithmic function family in musl libc demonstrates a high level of engineering practice in numerical computation:

Subtlety of Algorithm Design

1. 1. Unified Architecture: The three logarithmic functions share a core algorithm framework, reducing code duplication
2. 2. Precise Range Reduction: Quickly normalize inputs to the optimal range through bit manipulation
3. 3. Optimized Approximations: Use the most suitable approximation method based on the characteristics of different functions

Guarantee of Numerical Stability

Through multi-layer control mechanisms:

• • Numerical optimization of entries and polynomial coefficients
• • Precision maintenance of double precision intermediate calculations
• • Error control through high/low separation
• • Special handling of boundary cases

Engineering Excellence

1. 1. Balance of Performance and Precision: Optimize performance while ensuring precision
2. 2. Portability: Independent of specific hardware features
3. 3. Maintainability: Clear structural design and ample comments

The implementations of these three logarithmic functions are a model for numerical computation library design, with design philosophies including:

• • Layered Processing: Use different precision algorithms based on input ranges
• • Utilization of Mathematical Identities: Optimize calculations based on the algebraic properties of logarithmic functions
• • Error Analysis Driven: Guide algorithm selection based on rigorous error analysis

Compared to directly using the change of base formula, these specially optimized implementations:

• • Provide higher precision within their respective domains
• • Offer correct behavior for boundary cases
• • Maintain good performance while ensuring precision

This perfect balance between algorithmic innovation, numerical stability, and engineering practice makes the logarithmic function family in musl libc a reliable foundation for high-quality mathematical computation, providing a solid numerical basis for applications in scientific computing, engineering analysis, and machine learning.