diff options
Diffstat (limited to 'src/math/__expo2.c')
| -rw-r--r-- | src/math/__expo2.c | 49 |
1 files changed, 7 insertions, 42 deletions
diff --git a/src/math/__expo2.c b/src/math/__expo2.c index ef14e5f5..740ac680 100644 --- a/src/math/__expo2.c +++ b/src/math/__expo2.c @@ -1,51 +1,16 @@ -/* origin: FreeBSD /usr/src/lib/msun/src/k_exp.c */ -/*- - * Copyright (c) 2011 David Schultz <das@FreeBSD.ORG> - * All rights reserved. - * - * Redistribution and use in source and binary forms, with or without - * modification, are permitted provided that the following conditions - * are met: - * 1. Redistributions of source code must retain the above copyright - * notice, this list of conditions and the following disclaimer. - * 2. Redistributions in binary form must reproduce the above copyright - * notice, this list of conditions and the following disclaimer in the - * documentation and/or other materials provided with the distribution. - * - * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND - * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE - * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE - * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE - * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL - * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS - * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) - * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT - * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY - * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF - * SUCH DAMAGE. - */ - #include "libm.h" -/* - * We use exp(x) = exp(x - kln2) * 2**k, - * k is carefully chosen to minimize |exp(kln2) - 2**k| - */ -static const uint32_t k = 1799; -static const double kln2 = 1246.97177782734161156; +/* k is such that k*ln2 has minimal relative error and x - kln2 > log(DBL_MIN) */ +static const int k = 2043; +static const double kln2 = 0x1.62066151add8bp+10; -/* exp(x)/2 when x is huge */ +/* exp(x)/2 for x >= log(DBL_MAX), slightly better than 0.5*exp(x/2)*exp(x/2) */ double __expo2(double x) { double scale; - int n; - /* - * efficient scalbn(y, k-1): - * 2**(k-1) cannot be represented - * so we use that k-1 is even and scale in two steps - */ - n = (k - 1)/2; - INSERT_WORDS(scale, (0x3ff + n) << 20, 0); + /* note that k is odd and scale*scale overflows */ + INSERT_WORDS(scale, (uint32_t)(0x3ff + k/2) << 20, 0); + /* exp(x - k ln2) * 2**(k-1) */ return exp(x - kln2) * scale * scale; } |