/*
|
* e_logf.c - single precision log function
|
*
|
* Copyright (c) 2009-2018, Arm Limited.
|
* SPDX-License-Identifier: MIT
|
*/
|
|
/*
|
* Algorithm was once taken from Cody & Waite, but has been munged
|
* out of all recognition by SGT.
|
*/
|
|
#include <math.h>
|
#include <errno.h>
|
#include "math_private.h"
|
|
float
|
logf(float X)
|
{
|
int N = 0;
|
int aindex;
|
float a, x, s;
|
unsigned ix = fai(X);
|
|
if (__builtin_expect((ix - 0x00800000) >= 0x7f800000 - 0x00800000, 0)) {
|
if ((ix << 1) > 0xff000000) /* NaN */
|
return FLOAT_INFNAN(X);
|
if (ix == 0x7f800000) /* +inf */
|
return X;
|
if (X < 0) { /* anything negative */
|
return MATHERR_LOGF_NEG(X);
|
}
|
if (X == 0) {
|
return MATHERR_LOGF_0(X);
|
}
|
/* That leaves denormals. */
|
N = -23;
|
X *= 0x1p+23F;
|
ix = fai(X);
|
}
|
|
/*
|
* Separate X into three parts:
|
* - 2^N for some integer N
|
* - a number a of the form (1+k/8) for k=0,...,7
|
* - a residual which we compute as s = (x-a)/(x+a), for
|
* x=X/2^N.
|
*
|
* We pick the _nearest_ (N,a) pair, so that (x-a) has magnitude
|
* at most 1/16. Hence, we must round things that are just
|
* _below_ a power of two up to the next power of two, so this
|
* isn't as simple as extracting the raw exponent of the FP
|
* number. Instead we must grab the exponent together with the
|
* top few bits of the mantissa, and round (in integers) there.
|
*/
|
{
|
int rounded = ix + 0x00080000;
|
int Nnew = (rounded >> 23) - 127;
|
aindex = (rounded >> 20) & 7;
|
a = fhex(0x3f800000 + (aindex << 20));
|
N += Nnew;
|
x = fhex(ix - (Nnew << 23));
|
}
|
|
if (!N && !aindex) {
|
/*
|
* Use an alternative strategy for very small |x|, which
|
* avoids the 1ULP of relative error introduced in the
|
* computation of s. If our nearest (N,a) pair is N=0,a=1,
|
* that means we have -1/32 < x-a < 1/16, on which interval
|
* the ordinary series for log(1+z) (setting z-x-a) will
|
* converge adequately fast; so we can simply find an
|
* approximation to log(1+z)/z good on that interval and
|
* scale it by z on the way out.
|
*
|
* Coefficients generated by the command
|
|
./auxiliary/remez.jl --variable=z --suffix=f -- '-1/BigFloat(32)' '+1/BigFloat(16)' 3 0 '(log1p(x)-x)/x^2'
|
|
*/
|
float z = x - 1.0f;
|
float p = z*z * (
|
-4.999999767382730053173434595877399055021398381370452534949864039404089549132551e-01f+z*(3.333416379155995401749506866323446447523793085809161350343357014272193712456391e-01f+z*(-2.501299948811686421962724839011563450757435183422532362736159418564644404218257e-01f+z*(1.903576945606738444146078468935429697455230136403008172485495359631510244557255e-01f)))
|
);
|
|
return z + p;
|
}
|
|
/*
|
* Now we have N, a and x correct, so that |x-a| <= 1/16.
|
* Compute s.
|
*
|
* (Since |x+a| >= 2, this means that |s| will be at most 1/32.)
|
*/
|
s = (x - a) / (x + a);
|
|
/*
|
* The point of computing s = (x-a)/(x+a) was that this makes x
|
* equal to a * (1+s)/(1-s). So we can now compute log(x) by
|
* means of computing log((1+s)/(1-s)) (which has a more
|
* efficiently converging series), and adding log(a) which we
|
* obtain from a lookup table.
|
*
|
* So our full answer to log(X) is now formed by adding together
|
* N*log(2) + log(a) + log((1+s)/(1-s)).
|
*
|
* Now log((1+s)/(1-s)) has the exact Taylor series
|
*
|
* 2s + 2s^3/3 + 2s^5/5 + ...
|
*
|
* and what we do is to compute all but the first term of that
|
* as a polynomial approximation in s^2, then add on the first
|
* term - and all the other bits and pieces above - in
|
* precision-and-a-half so as to keep the total error down.
|
*/
|
{
|
float s2 = s*s;
|
|
/*
|
* We want a polynomial L(s^2) such that
|
*
|
* 2s + s^3*L(s^2) = log((1+s)/(1-s))
|
*
|
* => L(s^2) = (log((1+s)/(1-s)) - 2s) / s^3
|
*
|
* => L(z) = (log((1+sqrt(z))/(1-sqrt(z))) - 2*sqrt(z)) / sqrt(z)^3
|
*
|
* The required range of the polynomial is only [0,1/32^2].
|
*
|
* Our accuracy requirement for the polynomial approximation
|
* is that we don't want to introduce any error more than
|
* about 2^-23 times the _top_ bit of s. But the value of
|
* the polynomial has magnitude about s^3; and since |s| <
|
* 2^-5, this gives us |s^3/s| < 2^-10. In other words,
|
* our approximation only needs to be accurate to 13 bits or
|
* so before its error is lost in the noise when we add it
|
* to everything else.
|
*
|
* Coefficients generated by the command
|
|
./auxiliary/remez.jl --variable=s2 --suffix=f -- '0' '1/BigFloat(32^2)' 1 0 '(abs(x) < 1e-20 ? BigFloat(2)/3 + 2*x/5 + 2*x^2/7 + 2*x^3/9 : (log((1+sqrt(x))/(1-sqrt(x)))-2*sqrt(x))/sqrt(x^3))'
|
|
*/
|
float p = s * s2 * (
|
6.666666325680271091157649745099739739798210281016897722498744752867165138320995e-01f+s2*(4.002792299542401431889592846825025487338520940900492146195427243856292349188402e-01f)
|
);
|
|
static const float log2hi = 0x1.62ep-1F, log2lo = 0x1.0bfbe8p-15F;
|
static const float logahi[8] = { 0x0p+0F, 0x1.e26p-4F, 0x1.c8ep-3F, 0x1.46p-2F, 0x1.9f2p-2F, 0x1.f12p-2F, 0x1.1e8p-1F, 0x1.41cp-1F };
|
static const float logalo[8] = { 0x0p+0F, 0x1.076e2ap-16F, 0x1.f7c79ap-15F, 0x1.8bc21cp-14F, 0x1.23eccp-14F, 0x1.1ebf5ep-15F, 0x1.7d79c2p-15F, 0x1.8fe846p-13F };
|
return (N*log2hi + logahi[aindex]) + (2.0f*s + (N*log2lo + logalo[aindex] + p));
|
}
|
}
|
