diff options
Diffstat (limited to 'src/math')
| -rw-r--r-- | src/math/__expo2.c | 49 | ||||
| -rw-r--r-- | src/math/__expo2f.c | 49 |
2 files changed, 14 insertions, 84 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; } diff --git a/src/math/__expo2f.c b/src/math/__expo2f.c index 192838f7..5163e418 100644 --- a/src/math/__expo2f.c +++ b/src/math/__expo2f.c @@ -1,51 +1,16 @@ -/* origin: FreeBSD /usr/src/lib/msun/src/k_expf.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 = 235; -static const float kln2 = 162.88958740f; +/* k is such that k*ln2 has minimal relative error and x - kln2 > log(FLT_MIN) */ +static const int k = 235; +static const float kln2 = 0x1.45c778p+7f; -/* expf(x)/2 when x is huge */ +/* expf(x)/2 for x >= log(FLT_MAX), slightly better than 0.5f*expf(x/2)*expf(x/2) */ float __expo2f(float x) { float scale; - int n; - /* - * efficient scalbnf(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; - SET_FLOAT_WORD(scale, (0x7f + n) << 23); + /* note that k is odd and scale*scale overflows */ + SET_FLOAT_WORD(scale, (uint32_t)(0x7f + k/2) << 23); + /* exp(x - k ln2) * 2**(k-1) */ return expf(x - kln2) * scale * scale; } |