Introduction

In high-performance computing, leveraging SIMD (Single Instruction Multiple Data) instructions can lead to significant performance improvements. Among these, AVX2 (Advanced Vector Extensions 2) is a widely used instruction set for Intel and AMD processors. A common mathematical operation like the base-2 logarithm can be made highly efficient on these architectures when implemented using the __m256d data type, which represents a 256-bit vector containing four double-precision floating-point values.

Overview of AVX2 Instructions

AVX2 extends the earlier AVX instruction set, improving integer and floating-point performance. It allows simultaneous computations over wide 256-bit registers, meaning you can perform operations on four doubles or eight floats at once. For operations like log2, AVX2 provides vectorized arithmetic that can dramatically reduce computation time compared to scalar execution, often achieving close to 4x throughput improvement.

Understanding __m256d

The __m256d type is the intrinsic data type for operations on double-precision floating-point vectors:

• It holds exactly four 64-bit double values in a single 256-bit register.
• Operations involving this type process all four elements in parallel with a single instruction.
• Common operations include _mm256_add_pd, _mm256_mul_pd, _mm256_and_pd, and _mm256_castpd_si256 for bitwise manipulation.

Algorithm for Vectorized log2

The mathematical foundation relies on the IEEE 754 representation of floating-point numbers. Any positive double $x$ can be written as:

$$x = m \times 2^e$$

where $m$ is the mantissa (in the range $[1, 2)$) and $e$ is the exponent. Taking $\log_2$:

$$\log_2(x) = e + \log_2(m)$$

Since $m \in [1, 2)$, we have $\log_2(m) \in [0, 1)$, which is a small range ideal for polynomial approximation.

Step 1: Extract Exponent and Mantissa

For IEEE 754 doubles, the exponent is stored in bits 52 through 62, biased by 1023. Using AVX2 integer operations, we extract the exponent and normalize the mantissa:

1// Extract exponent (as integer, then convert to double)
2__m256i xi = _mm256_castpd_si256(x);
3__m256i exp_bits = _mm256_srli_epi64(xi, 52);
4__m256d exponent = _mm256_sub_pd(
5    _mm256_cvtepi64_pd(exp_bits),  // convert to double
6    _mm256_set1_pd(1023.0)          // remove bias
7);
8
9// Normalize mantissa to [1, 2)
10__m256i mantissa_mask = _mm256_set1_epi64x(0x000FFFFFFFFFFFFFLL);
11__m256i mantissa_bits = _mm256_or_si256(
12    _mm256_and_si256(xi, mantissa_mask),
13    _mm256_set1_epi64x(0x3FF0000000000000LL)  // exponent = 0 (biased 1023)
14);
15__m256d m = _mm256_castsi256_pd(mantissa_bits);

Step 2: Polynomial Approximation of $\log_2(m)$

With $m \in [1, 2)$, we approximate $\log_2(m)$ using a minimax polynomial. A degree-5 polynomial provides roughly 20 bits of precision, while degree-7 reaches close to full double precision. For many applications, a degree-5 or degree-6 polynomial strikes a good balance:

$$\log_2(m) \approx c_0 + c_1(m-1) + c_2(m-1)^2 + c_3(m-1)^3 + \ldots$$

Using Horner's method for efficient evaluation:

1__m256d t = _mm256_sub_pd(m, _mm256_set1_pd(1.0));
2
3// Horner's method: ((c3 * t + c2) * t + c1) * t + c0
4__m256d result = _mm256_set1_pd(c5);
5result = _mm256_fmadd_pd(result, t, _mm256_set1_pd(c4));
6result = _mm256_fmadd_pd(result, t, _mm256_set1_pd(c3));
7result = _mm256_fmadd_pd(result, t, _mm256_set1_pd(c2));
8result = _mm256_fmadd_pd(result, t, _mm256_set1_pd(c1));
9result = _mm256_fmadd_pd(result, t, _mm256_set1_pd(c0));

The _mm256_fmadd_pd instruction (fused multiply-add) computes $a \times b + c$ in a single instruction with only one rounding step, improving both speed and precision.

Step 3: Combine Results

The final result combines the integer exponent and the polynomial approximation:

__m256d log2_result = _mm256_add_pd(exponent, result);

Precision vs Performance Trade-offs

For graphics and audio applications, degree 3 to 5 is often sufficient. Scientific computing typically requires degree 7 or higher.

Handling Edge Cases

A production implementation must handle several special cases:

• Zero: $\log_2(0) = -\infty$. Detect using comparison and blend with _mm256_set1_pd(-INFINITY).
• Negative numbers: Result is NaN. Detect with a sign-bit check.
• Subnormal numbers: Very small numbers with a zero exponent field. Multiply by $2^{52}$ first, compute $\log_2$, then subtract 52 from the result.
• Infinity: $\log_2(\infty) = \infty$. Pass through unchanged.
• NaN: Input NaN should propagate to output NaN.

Performance Considerations

• Throughput: A well-optimized AVX2 log2 processes four doubles per call with a latency of roughly 15 to 25 clock cycles, compared to about 80 to 120 cycles for a scalar log2 from the math library.
• Alignment: Ensure data is 32-byte aligned for _mm256_load_pd instead of _mm256_loadu_pd to avoid penalties on some architectures.
• Interleaving: When computing log2 over large arrays, process multiple vectors per loop iteration to hide instruction latency and keep execution units busy.

Summary

Implementing log2 for __m256d in AVX2 follows a three-step process: extract the exponent from the IEEE 754 representation, approximate $\log_2(m)$ for the normalized mantissa using a minimax polynomial evaluated via Horner's method with fused multiply-add, and combine the two parts. The polynomial degree controls the precision-performance trade-off. Proper handling of edge cases (zero, negative, subnormal, infinity, NaN) is essential for production use. This approach delivers roughly 4x throughput improvement over scalar implementations.