ldexp_exp.jl

"""
    frexp_exp(x)

Compute ``exp(x)``, scaled to avoid spurious overflow.  An exponent is
returned separately in ``expt``.

# Examples
```jldoctest
julia> frexp_exp(1354.9)
(1.4678022081231224e308, 931)
```
"""
function frexp_exp(x::Float64)
    # Input:  ln(DBL_MAX) <= x < ln(2 * DBL_MAX / DBL_MIN_DENORM) ~= 1454.91
    # Output: 2**1023 <= y < 2**1024
    # We use exp(x) = exp(x - kln2) * 2**k, carefully chosen to
    # minimize |exp(kln2) - 2**k|.  We also scale the exponent of
    # exp_x to MAX_EXP so that the result can be multiplied by
    # a tiny number without losing accuracy due to denormalization.

    exp_x = exp(x - 1246.97177782734161156)
    hx, lx = words(exp_x)
    expt = (hx >> 20) - (0x3ff + 1023) + 0x00000707

    exp_x = from_words((hx & 0xfffff) | ((0x3ff + 1023) << 20), lx)
    return exp_x, expt
end

"""
    ldexp_exp(x, n)

Compute ``exp(x) * 2^n``. They are intended for large arguments (real part >= ln(DBL_MAX))
where care is needed to avoid overflow.

# Examples
```jldoctest
julia> frexp_exp(1354.9)
(1.4678022081231224e308, 931)
```
"""
function ldexp_exp(x::Float64, expt)
    # The present implementation is narrowly tailored for our hyperbolic and
    # exponential functions.  We assume expt is small (0 or -1), and the caller
    # has filtered out very large x, for which overflow would be inevitable.
    exp_x, ex_expt = frexp_exp(x)
    expt += ex_expt
    scale = from_words((0x3ff + expt) << 20, 0)
    return exp_x * scale
end