Implementation:ClickHouse ClickHouse Musl Powl
| Knowledge Sources | |
|---|---|
| Domains | Mathematics, Portability |
| Last Updated | 2026-02-08 00:00 GMT |
Overview
Long double precision power function using lookup tables and pseudo extended precision arithmetic.
Description
This implementation computes `powl(x, y)` for long double arguments using the formula `x^y = exp(y * log(x))`. It employs a 32-entry lookup table of `2^(-i/32)` values and extended precision techniques to achieve several extra bits of accuracy beyond native long double precision. The code handles integer exponents efficiently through binary exponentiation, uses rational approximations for logarithm and exponential calculations, and includes comprehensive special case handling for edge conditions. Originally from OpenBSD and Stephen L. Moshier.
Usage
Use this function when computing powers with extended precision requirements, particularly in numerical analysis, scientific computing, or financial calculations where double precision is insufficient and platform-specific long double behavior needs to be consistent.
Code Reference
Source Location
- Repository: ClickHouse
- File: base/glibc-compatibility/musl/powl.c
- Lines: 1-525
Signature
long double powl(long double x, long double y);
// Internal helpers
static long double reducl(long double x);
static long double powil(long double x, int n);
Import
#include <math.h>
I/O Contract
| Input Pattern | Output | Precision Notes |
|---|---|---|
| powl(x, y) for normal values | x^y | Relative error ~2.3e-21 * y |
| powl(±0, negative) | ±∞ | Sign depends on exponent parity |
| powl(±0, positive) | ±0 | Sign depends on exponent parity |
| powl(-x, integer) | Computed via `powil` | y| < 32768 |
| powl(x, y) where x ≥ 2^16384 | Overflow → ∞ | Largest representable exponent |
| powl(x, y) where result < 2^-16384 | Underflow → 0 | Smallest representable value |
| powl(NaN, y) or powl(x, NaN) | NaN | Except special cases like 1^NaN |
| powl(1, y) for any y | 1 | Even for infinite y |
| powl(-1, ±∞) | 1 | Oscillating limit |
Usage Examples
// High-precision power computation
long double x = 1.5L;
long double y = 100.0L;
long double result = powl(x, y);
printf("1.5^100 = %.20Lf\n", result);
// Integer exponents use fast path
long double base = 2.0L;
long double power_of_two = powl(base, 63.0L); // Uses powil for efficiency
// Extended precision for small values
long double small = 0.9999L;
long double many_times = powl(small, 10000.0L);
// Maintains precision through extended arithmetic
// Root extraction with high precision
long double cube_root = powl(27.0L, 1.0L/3.0L);
long double tenth_root = powl(1024.0L, 0.1L);
// Special cases
long double one = powl(1.0L, INFINITY); // 1.0
long double neg_one_inf = powl(-1.0L, INFINITY); // 1.0
// Negative base with integer exponent
long double neg_result = powl(-2.0L, 5.0L); // -32.0
long double pos_result = powl(-2.0L, 6.0L); // 64.0
// Overflow detection
long double huge_exp = 20000.0L;
long double overflow_result = powl(2.0L, huge_exp);
if (isinf(overflow_result)) {
printf("Result overflowed\n");
}
// Compound interest with extended precision
long double compound_interest_precise(long double principal,
long double rate,
long double periods) {
return principal * powl(1.0L + rate, periods);
}
// Growth calculation over long timescales
long double exponential_growth(long double initial,
long double rate,
long double time) {
return initial * powl(M_El, rate * time); // e^(rate*time)
}
// Check for reasonable exponent range
int is_safe_exponent(long double y) {
return fabsl(y) < 10000.0L; // Avoid extreme overflow/underflow
}
// Computing very large factorials in log space then exponentiating
long double approx_factorial(int n) {
if (n < 0) return 0.0L/0.0L; // NaN
if (n == 0 || n == 1) return 1.0L;
// Stirling approximation: n! ≈ sqrt(2πn) * (n/e)^n
long double log_fact = 0.5L * logl(2.0L * M_PIl * n)
+ n * logl(n) - n;
return powl(M_El, log_fact);
}