lib/msun/src/k_exp.c

   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 <complex.h>
  31
  32 #include "math.h"
  33 #include "math_private.h"
  34
  35 static const uint32_t k = 1799;         /* constant for reduction */
  36 static const double kln2 =  1246.97177782734161156;     /* k * ln2 */
  37
  38 /*
  39  * Compute exp(x), scaled to avoid spurious overflow.  An exponent is
  40  * returned separately in 'expt'.
  41  *
  42  * Input:  ln(DBL_MAX) <= x < ln(2 * DBL_MAX / DBL_MIN_DENORM) ~= 1454.91
  43  * Output: 2**1023 <= y < 2**1024
  44  */
  45 static double
  46 __frexp_exp(double x, int *expt)
  47 {
  48         double exp_x;
  49         uint32_t hx;
  50
  51         /*
  52          * We use exp(x) = exp(x - kln2) * 2**k, carefully chosen to
  53          * minimize |exp(kln2) - 2**k|.  We also scale the exponent of
  54          * exp_x to MAX_EXP so that the result can be multiplied by
  55          * a tiny number without losing accuracy due to denormalization.
  56          */
  57         exp_x = exp(x - kln2);
  58         GET_HIGH_WORD(hx, exp_x);
  59         *expt = (hx >> 20) - (0x3ff + 1023) + k;
  60         SET_HIGH_WORD(exp_x, (hx & 0xfffff) | ((0x3ff + 1023) << 20));
  61         return (exp_x);
  62 }
  63
  64 /*
  65  * __ldexp_exp(x, expt) and __ldexp_cexp(x, expt) compute exp(x) * 2**expt.
  66  * They are intended for large arguments (real part >= ln(DBL_MAX))
  67  * where care is needed to avoid overflow.
  68  *
  69  * The present implementation is narrowly tailored for our hyperbolic and
  70  * exponential functions.  We assume expt is small (0 or -1), and the caller
  71  * has filtered out very large x, for which overflow would be inevitable.
  72  */
  73
  74 double
  75 __ldexp_exp(double x, int expt)
  76 {
  77         double exp_x, scale;
  78         int ex_expt;
  79
  80         exp_x = __frexp_exp(x, &ex_expt);
  81         expt += ex_expt;
  82         INSERT_WORDS(scale, (0x3ff + expt) << 20, 0);
  83         return (exp_x * scale);
  84 }
  85
  86 double complex
  87 __ldexp_cexp(double complex z, int expt)
  88 {
  89         double x, y, exp_x, scale1, scale2;
  90         int ex_expt, half_expt;
  91
  92         x = creal(z);
  93         y = cimag(z);
  94         exp_x = __frexp_exp(x, &ex_expt);
  95         expt += ex_expt;
  96
  97         /*
  98          * Arrange so that scale1 * scale2 == 2**expt.  We use this to
  99          * compensate for scalbn being horrendously slow.
 100          */
 101         half_expt = expt / 2;
 102         INSERT_WORDS(scale1, (0x3ff + half_expt) << 20, 0);
 103         half_expt = expt - half_expt;
 104         INSERT_WORDS(scale2, (0x3ff + half_expt) << 20, 0);
 105
 106         return (cpack(cos(y) * exp_x * scale1 * scale2,
 107             sin(y) * exp_x * scale1 * scale2));
 108 }