Line data    Source code

       1             : /*-
       2             :  * Copyright (c) 2011 David Schultz <das@FreeBSD.ORG>
       3             :  * All rights reserved.
       4             :  *
       5             :  * Redistribution and use in source and binary forms, with or without
       6             :  * modification, are permitted provided that the following conditions
       7             :  * are met:
       8             :  * 1. Redistributions of source code must retain the above copyright
       9             :  *    notice, this list of conditions and the following disclaimer.
      10             :  * 2. Redistributions in binary form must reproduce the above copyright
      11             :  *    notice, this list of conditions and the following disclaimer in the
      12             :  *    documentation and/or other materials provided with the distribution.
      13             :  *
      14             :  * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
      15             :  * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
      16             :  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
      17             :  * ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
      18             :  * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
      19             :  * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
      20             :  * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
      21             :  * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
      22             :  * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
      23             :  * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
      24             :  * SUCH DAMAGE.
      25             :  */
      26             : 
      27             : //#include <sys/cdefs.h>
      28             : //__FBSDID("$FreeBSD$");
      29             : 
      30             :  #include "math_private.h"
      31             : 
      32             : static const uint32_t k = 1799;         /* constant for reduction */
      33             : static const double kln2 =  1246.97177782734161156;     /* k * ln2 */
      34             : 
      35             : /*
      36             :  * Compute exp(x), scaled to avoid spurious overflow.  An exponent is
      37             :  * returned separately in 'expt'.
      38             :  *
      39             :  * Input:  ln(DBL_MAX) <= x < ln(2 * DBL_MAX / DBL_MIN_DENORM) ~= 1454.91
      40             :  * Output: 2**1023 <= y < 2**1024
      41             :  */
      42             : static double
      43           0 : __frexp_exp(double x, int *expt)
      44             : {
      45             :         double exp_x;
      46             :         uint32_t hx;
      47             : 
      48             :         /*
      49             :          * We use exp(x) = exp(x - kln2) * 2**k, carefully chosen to
      50             :          * minimize |exp(kln2) - 2**k|.  We also scale the exponent of
      51             :          * exp_x to MAX_EXP so that the result can be multiplied by
      52             :          * a tiny number without losing accuracy due to denormalization.
      53             :          */
      54           0 :         exp_x = exp(x - kln2);
      55           0 :         GET_HIGH_WORD(hx, exp_x);
      56           0 :         *expt = (hx >> 20) - (0x3ff + 1023) + k;
      57           0 :         SET_HIGH_WORD(exp_x, (hx & 0xfffff) | ((0x3ff + 1023) << 20));
      58           0 :         return (exp_x);
      59             : }
      60             : 
      61             : /*
      62             :  * __ldexp_exp(x, expt) and __ldexp_cexp(x, expt) compute exp(x) * 2**expt.
      63             :  * They are intended for large arguments (real part >= ln(DBL_MAX))
      64             :  * where care is needed to avoid overflow.
      65             :  *
      66             :  * The present implementation is narrowly tailored for our hyperbolic and
      67             :  * exponential functions.  We assume expt is small (0 or -1), and the caller
      68             :  * has filtered out very large x, for which overflow would be inevitable.
      69             :  */
      70             : 
      71             : double
      72           0 : __ldexp_exp(double x, int expt)
      73             : {
      74             :         double exp_x, scale;
      75             :         int ex_expt;
      76             : 
      77           0 :         exp_x = __frexp_exp(x, &ex_expt);
      78           0 :         expt += ex_expt;
      79           0 :         INSERT_WORDS(scale, (0x3ff + expt) << 20, 0);
      80           0 :         return (exp_x * scale);
      81             : }