25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis/*
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * CDDL HEADER START
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis *
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * The contents of this file are subject to the terms of the
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * Common Development and Distribution License (the "License").
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * You may not use this file except in compliance with the License.
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis *
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * or http://www.opensolaris.org/os/licensing.
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * See the License for the specific language governing permissions
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * and limitations under the License.
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis *
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * When distributing Covered Code, include this CDDL HEADER in each
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * file and include the License file at usr/src/OPENSOLARIS.LICENSE.
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * If applicable, add the following below this CDDL HEADER, with the
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * fields enclosed by brackets "[]" replaced with your own identifying
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * information: Portions Copyright [yyyy] [name of copyright owner]
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis *
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * CDDL HEADER END
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis/*
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * Copyright 2011 Nexenta Systems, Inc. All rights reserved.
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis/*
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * Copyright 2006 Sun Microsystems, Inc. All rights reserved.
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * Use is subject to license terms.
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis/*
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * expl(x)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * Table driven method
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * Written by K.C. Ng, November 1988.
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * Algorithm :
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * 1. Argument Reduction: given the input x, find r and integer k
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * and j such that
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * x = (32k+j)*ln2 + r, |r| <= (1/64)*ln2 .
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis *
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * 2. expl(x) = 2^k * (2^(j/32) + 2^(j/32)*expm1(r))
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * Note:
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * a. expm1(r) = (2r)/(2-R), R = r - r^2*(t1 + t2*r^2)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * b. 2^(j/32) is represented as
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * _TBL_expl_hi[j]+_TBL_expl_lo[j]
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * where
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * _TBL_expl_hi[j] = 2^(j/32) rounded
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * _TBL_expl_lo[j] = 2^(j/32) - _TBL_expl_hi[j].
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis *
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * Special cases:
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * expl(INF) is INF, expl(NaN) is NaN;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * expl(-INF)= 0;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * for finite argument, only expl(0)=1 is exact.
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis *
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * Accuracy:
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * according to an error analysis, the error is always less than
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * an ulp (unit in the last place).
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis *
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * Misc. info.
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * For 113 bit long double
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * if x > 1.135652340629414394949193107797076342845e+4
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * then expl(x) overflow;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * if x < -1.143346274333629787883724384345262150341e+4
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * then expl(x) underflow
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis *
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * Constants:
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * Only decimal values are given. We assume that the compiler will convert
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * from decimal to binary accurately enough to produce the correct
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * hexadecimal values.
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis
ddc0e0b53c661f6e439e3b7072b3ef353eadb4afRichard Lowe#pragma weak __expl = expl
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis#include "libm.h"
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtisextern const long double _TBL_expl_hi[], _TBL_expl_lo[];
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtisstatic const long double
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtisone = 1.0L,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtistwo = 2.0L,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtisln2_64 = 1.083042469624914545964425189778400898568e-2L,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtisovflthreshold = 1.135652340629414394949193107797076342845e+4L,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtisunflthreshold = -1.143346274333629787883724384345262150341e+4L,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtisinvln2_32 = 4.616624130844682903551758979206054839765e+1L,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtisln2_32hi = 2.166084939249829091928849858592451515688e-2L,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtisln2_32lo = 5.209643502595475652782654157501186731779e-27L;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis/* rational approximation coeffs for [-(ln2)/64,(ln2)/64] */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtisstatic const long double
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtist1 = 1.666666666666666666666666666660876387437e-1L,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtist2 = -2.777777777777777777777707812093173478756e-3L,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtist3 = 6.613756613756613482074280932874221202424e-5L,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtist4 = -1.653439153392139954169609822742235851120e-6L,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtist5 = 4.175314851769539751387852116610973796053e-8L;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtislong double
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtisexpl(long double x) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis int *px = (int *) &x, ix, j, k, m;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis long double t, r;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis ix = px[0]; /* high word of x */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis if (ix >= 0x7fff0000)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis return (x + x); /* NaN of +inf */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis if (((unsigned) ix) >= 0xffff0000)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis return (-one / x); /* NaN or -inf */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis if ((ix & 0x7fffffff) < 0x3fc30000) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis if ((int) x < 1)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis return (one + x); /* |x|<2^-60 */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis }
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis if (ix > 0) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis if (x > ovflthreshold)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis return (scalbnl(x, 20000));
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis k = (int) (invln2_32 * (x + ln2_64));
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis } else {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis if (x < unflthreshold)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis return (scalbnl(-x, -40000));
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis k = (int) (invln2_32 * (x - ln2_64));
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis }
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis j = k&0x1f;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis m = k>>5;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis t = (long double) k;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis x = (x - t * ln2_32hi) - t * ln2_32lo;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis t = x * x;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis r = (x - t * (t1 + t * (t2 + t * (t3 + t * (t4 + t * t5))))) - two;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis x = _TBL_expl_hi[j] - ((_TBL_expl_hi[j] * (x + x)) / r -
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis _TBL_expl_lo[j]);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis return (scalbnl(x, m));
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis}