compiler/stdc/math/ld80/s_log1pl.c

   1 /*      $OpenBSD: s_log1pl.c,v 1.3 2013/11/12 20:35:19 martynas Exp $   */
   2
   3 /*
   4  * Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net>
   5  *
   6  * Permission to use, copy, modify, and distribute this software for any
   7  * purpose with or without fee is hereby granted, provided that the above
   8  * copyright notice and this permission notice appear in all copies.
   9  *
  10  * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
  11  * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
  12  * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
  13  * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
  14  * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
  15  * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
  16  * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
  17  */
  18
  19 /*                                                      log1pl.c
  20  *
  21  *      Relative error logarithm
  22  *      Natural logarithm of 1+x, long double precision
  23  *
  24  *
  25  *
  26  * SYNOPSIS:
  27  *
  28  * long double x, y, log1pl();
  29  *
  30  * y = log1pl( x );
  31  *
  32  *
  33  *
  34  * DESCRIPTION:
  35  *
  36  * Returns the base e (2.718...) logarithm of 1+x.
  37  *
  38  * The argument 1+x is separated into its exponent and fractional
  39  * parts.  If the exponent is between -1 and +1, the logarithm
  40  * of the fraction is approximated by
  41  *
  42  *     log(1+x) = x - 0.5 x^2 + x^3 P(x)/Q(x).
  43  *
  44  * Otherwise, setting  z = 2(x-1)/x+1),
  45  *
  46  *     log(x) = z + z^3 P(z)/Q(z).
  47  *
  48  *
  49  *
  50  * ACCURACY:
  51  *
  52  *                      Relative error:
  53  * arithmetic   domain     # trials      peak         rms
  54  *    IEEE     -1.0, 9.0    100000      8.2e-20    2.5e-20
  55  *
  56  * ERROR MESSAGES:
  57  *
  58  * log singularity:  x-1 = 0; returns -INFINITY
  59  * log domain:       x-1 < 0; returns NAN
  60  */
  61
  62 #include "math.h"
  63
  64 #include "math_private.h"
  65
  66 /* Coefficients for log(1+x) = x - x^2 / 2 + x^3 P(x)/Q(x)
  67  * 1/sqrt(2) <= x < sqrt(2)
  68  * Theoretical peak relative error = 2.32e-20
  69  */
  70
  71 static const long double P[] = {
  72  4.5270000862445199635215E-5L,
  73  4.9854102823193375972212E-1L,
  74  6.5787325942061044846969E0L,
  75  2.9911919328553073277375E1L,
  76  6.0949667980987787057556E1L,
  77  5.7112963590585538103336E1L,
  78  2.0039553499201281259648E1L,
  79 };
  80 static const long double Q[] = {
  81 /* 1.0000000000000000000000E0,*/
  82  1.5062909083469192043167E1L,
  83  8.3047565967967209469434E1L,
  84  2.2176239823732856465394E2L,
  85  3.0909872225312059774938E2L,
  86  2.1642788614495947685003E2L,
  87  6.0118660497603843919306E1L,
  88 };
  89
  90 /* Coefficients for log(x) = z + z^3 P(z^2)/Q(z^2),
  91  * where z = 2(x-1)/(x+1)
  92  * 1/sqrt(2) <= x < sqrt(2)
  93  * Theoretical peak relative error = 6.16e-22
  94  */
  95
  96 static const long double R[4] = {
  97  1.9757429581415468984296E-3L,
  98 -7.1990767473014147232598E-1L,
  99  1.0777257190312272158094E1L,
 100 -3.5717684488096787370998E1L,
 101 };
 102 static const long double S[4] = {
 103 /* 1.00000000000000000000E0L,*/
 104 -2.6201045551331104417768E1L,
 105  1.9361891836232102174846E2L,
 106 -4.2861221385716144629696E2L,
 107 };
 108 static const long double C1 = 6.9314575195312500000000E-1L;
 109 static const long double C2 = 1.4286068203094172321215E-6L;
 110
 111 #define SQRTH 0.70710678118654752440L
 112
 113 long double
 114 log1pl(long double xm1)
 115 {
 116 long double x, y, z;
 117 int e;
 118
 119 if( isnan(xm1) )
 120         return(xm1);
 121 if( xm1 == INFINITY )
 122         return(xm1);
 123 if(xm1 == 0.0)
 124         return(xm1);
 125
 126 x = xm1 + 1.0L;
 127
 128 /* Test for domain errors.  */
 129 if( x <= 0.0L )
 130         {
 131         if( x == 0.0L )
 132                 return( -INFINITY );
 133         else
 134                 return( NAN );
 135         }
 136
 137 /* Separate mantissa from exponent.
 138    Use frexp so that denormal numbers will be handled properly.  */
 139 x = frexpl( x, &e );
 140
 141 /* logarithm using log(x) = z + z^3 P(z)/Q(z),
 142    where z = 2(x-1)/x+1)  */
 143 if( (e > 2) || (e < -2) )
 144 {
 145 if( x < SQRTH )
 146         { /* 2( 2x-1 )/( 2x+1 ) */
 147         e -= 1;
 148         z = x - 0.5L;
 149         y = 0.5L * z + 0.5L;
 150         }
 151 else
 152         { /*  2 (x-1)/(x+1)   */
 153         z = x - 0.5L;
 154         z -= 0.5L;
 155         y = 0.5L * x  + 0.5L;
 156         }
 157 x = z / y;
 158 z = x*x;
 159 z = x * ( z * __polevll( z, (void *)R, 3 ) / __p1evll( z, (void *)S, 3 ) );
 160 z = z + e * C2;
 161 z = z + x;
 162 z = z + e * C1;
 163 return( z );
 164 }
 165
 166
 167 /* logarithm using log(1+x) = x - .5x**2 + x**3 P(x)/Q(x) */
 168
 169 if( x < SQRTH )
 170         {
 171         e -= 1;
 172         if (e != 0)
 173           x = 2.0 * x - 1.0L;
 174         else
 175           x = xm1;
 176         }
 177 else
 178         {
 179           if (e != 0)
 180             x = x - 1.0L;
 181           else
 182             x = xm1;
 183         }
 184 z = x*x;
 185 y = x * ( z * __polevll( x, (void *)P, 6 ) / __p1evll( x, (void *)Q, 6 ) );
 186 y = y + e * C2;
 187 z = y - 0.5 * z;
 188 z = z + x;
 189 z = z + e * C1;
 190 return( z );
 191 }