common/Q/expm1l.c

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
ddc0e0b53c661f6e439e3b7072b3ef353eadb4afRichard Lowe#pragma weak __expm1l = expm1l
ddc0e0b53c661f6e439e3b7072b3ef353eadb4afRichard Lowe
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis#if !defined(__sparc)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis#error Unsupported architecture
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis#endif
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis/*
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * expm1l(x)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis *
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * Table driven method
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * Written by K.C. Ng, June 1995.
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * Algorithm :
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis *  1. expm1(x) = x if x<2**-114
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis *  2. if |x| <= 0.0625 = 1/16, use approximation
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis *      expm1(x) = x + x*P/(2-P)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * where
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis *  P = x - z*(P1+z*(P2+z*(P3+z*(P4+z*(P5+z*P6+z*P7))))), z = x*x;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * (this formula is derived from
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis *  2-P+x = R = x*(exp(x)+1)/(exp(x)-1) ~ 2 + x*x/6 - x^4/360 + ...)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis *
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * P1 =   1.66666666666666666666666666666638500528074603030e-0001
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * P2 =  -2.77777777777777777777777759668391122822266551158e-0003
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * P3 =   6.61375661375661375657437408890138814721051293054e-0005
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * P4 =  -1.65343915343915303310185228411892601606669528828e-0006
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * P5 =   4.17535139755122945763580609663414647067443411178e-0008
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * P6 =  -1.05683795988668526689182102605260986731620026832e-0009
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * P7 =   2.67544168821852702827123344217198187229611470514e-0011
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis *
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * Accuracy: |R-x*(exp(x)+1)/(exp(x)-1)|<=2**-119.13
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis *
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis *  3. For 1/16 < |x| < 1.125, choose x(+-i) ~ +-(i+4.5)/64, i=0,..,67
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis *     since
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis *      exp(x) = exp(xi+(x-xi))= exp(xi)*exp((x-xi))
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis *     we have
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis *      expm1(x) = expm1(xi)+(exp(xi))*(expm1(x-xi))
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis *     where
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis *      |s=x-xi| <= 1/128
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis *     and
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis *  expm1(s)=2s/(2-R), R= s-s^2*(T1+s^2*(T2+s^2*(T3+s^2*(T4+s^2*T5))))
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis *
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * T1 =   1.666666666666666666666666666660876387437e-1L,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * T2 =  -2.777777777777777777777707812093173478756e-3L,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * T3 =   6.613756613756613482074280932874221202424e-5L,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * T4 =  -1.653439153392139954169609822742235851120e-6L,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * T5 =   4.175314851769539751387852116610973796053e-8L;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis *
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis *  4. For |x| >= 1.125, return exp(x)-1.
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis *      (see algorithm for exp)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis *
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * Special cases:
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis *  expm1l(INF) is INF, expm1l(NaN) is NaN;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis *  expm1l(-INF)= -1;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis *  for finite argument, only expm1l(0)=0 is exact.
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis *
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * Accuracy:
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis *  according to an error analysis, the error is always less than
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis *  2 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 expm1l(x) overflow;
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
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis#include "libm.h"
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtisextern const long double _TBL_expl_hi[], _TBL_expl_lo[];
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtisextern const long double _TBL_expm1lx[], _TBL_expm1l[];
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtisstatic const long double
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    zero        = +0.0L,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    one     = +1.0L,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    two     = +2.0L,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    ln2_64      = +1.083042469624914545964425189778400898568e-2L,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    ovflthreshold   = +1.135652340629414394949193107797076342845e+4L,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    invln2_32   = +4.616624130844682903551758979206054839765e+1L,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    ln2_32hi    = +2.166084939249829091928849858592451515688e-2L,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    ln2_32lo    = +5.209643502595475652782654157501186731779e-27L,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    huge        = +1.0e4000L,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    tiny        = +1.0e-4000L,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    P1 = +1.66666666666666666666666666666638500528074603030e-0001L,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    P2 = -2.77777777777777777777777759668391122822266551158e-0003L,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    P3 = +6.61375661375661375657437408890138814721051293054e-0005L,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    P4 = -1.65343915343915303310185228411892601606669528828e-0006L,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    P5 = +4.17535139755122945763580609663414647067443411178e-0008L,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    P6 = -1.05683795988668526689182102605260986731620026832e-0009L,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    P7 = +2.67544168821852702827123344217198187229611470514e-0011L,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis/* rational approximation coeffs for [-(ln2)/64,(ln2)/64] */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    T1 = +1.666666666666666666666666666660876387437e-1L,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    T2 = -2.777777777777777777777707812093173478756e-3L,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    T3 = +6.613756613756613482074280932874221202424e-5L,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    T4 = -1.653439153392139954169609822742235851120e-6L,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    T5 = +4.175314851769539751387852116610973796053e-8L;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtislong double
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtisexpm1l(long double x) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    int hx, ix, j, k, m;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    long double t, r, s, w;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    hx = ((int *) &x)[HIXWORD];
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    ix = hx & ~0x80000000;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    if (ix >= 0x7fff0000) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        if (x != x)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis            return (x + x); /* NaN */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        if (x < zero)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis            return (-one);  /* -inf */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        return (x); /* +inf */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    }
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    if (ix < 0x3fff4000) {  /* |x| < 1.25 */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        if (ix < 0x3ffb0000) {  /* |x| < 0.0625 */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis            if (ix < 0x3f8d0000) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis                if ((int) x == 0)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis                    return (x); /* |x|<2^-114 */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis            }
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis            t = x * x;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis            r = (x - t * (P1 + t * (P2 + t * (P3 + t * (P4 + t *
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis                (P5 + t * (P6 + t * P7)))))));
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis            return (x + (x * r) / (two - r));
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        }
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        /* compute i = [64*x] */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        m = 0x4009 - (ix >> 16);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        j = ((ix & 0x0000ffff) | 0x10000) >> m; /* j=4,...,67 */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        if (hx < 0)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis            j += 82;            /* negative */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        s = x - _TBL_expm1lx[j];
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        t = s * s;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        r = s - t * (T1 + t * (T2 + t * (T3 + t * (T4 + t * T5))));
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        r = (s + s) / (two - r);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        w = _TBL_expm1l[j];
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        return (w + (w + one) * r);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    }
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    if (hx > 0) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        if (x > ovflthreshold)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis            return (huge * huge);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        k = (int) (invln2_32 * (x + ln2_64));
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    } else {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        if (x < -80.0)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis            return (tiny - x / x);
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) - one);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis}