compiler/stdc/math/ld80/e_log10l.c

   1 /*      $OpenBSD: e_log10l.c,v 1.2 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 /*                                                      log10l.c
  20  *
  21  *      Common logarithm, long double precision
  22  *
  23  *
  24  *
  25  * SYNOPSIS:
  26  *
  27  * long double x, y, log10l();
  28  *
  29  * y = log10l( x );
  30  *
  31  *
  32  *
  33  * DESCRIPTION:
  34  *
  35  * Returns the base 10 logarithm of x.
  36  *
  37  * The argument is separated into its exponent and fractional
  38  * parts.  If the exponent is between -1 and +1, the logarithm
  39  * of the fraction is approximated by
  40  *
  41  *     log(1+x) = x - 0.5 x**2 + x**3 P(x)/Q(x).
  42  *
  43  * Otherwise, setting  z = 2(x-1)/x+1),
  44  *
  45  *     log(x) = z + z**3 P(z)/Q(z).
  46  *
  47  *
  48  *
  49  * ACCURACY:
  50  *
  51  *                      Relative error:
  52  * arithmetic   domain     # trials      peak         rms
  53  *    IEEE      0.5, 2.0     30000      9.0e-20     2.6e-20
  54  *    IEEE     exp(+-10000)  30000      6.0e-20     2.3e-20
  55  *
  56  * In the tests over the interval exp(+-10000), the logarithms
  57  * of the random arguments were uniformly distributed over
  58  * [-10000, +10000].
  59  *
  60  * ERROR MESSAGES:
  61  *
  62  * log singularity:  x = 0; returns MINLOG
  63  * log domain:       x < 0; returns MINLOG
  64  */
  65
  66 #include "math.h"
  67
  68 #include "math_private.h"
  69
  70 /* Coefficients for log(1+x) = x - x**2/2 + x**3 P(x)/Q(x)
  71  * 1/sqrt(2) <= x < sqrt(2)
  72  * Theoretical peak relative error = 6.2e-22
  73  */
  74 static const long double P[] = {
  75  4.9962495940332550844739E-1L,
  76  1.0767376367209449010438E1L,
  77  7.7671073698359539859595E1L,
  78  2.5620629828144409632571E2L,
  79  4.2401812743503691187826E2L,
  80  3.4258224542413922935104E2L,
  81  1.0747524399916215149070E2L,
  82 };
  83 static const long double Q[] = {
  84 /* 1.0000000000000000000000E0,*/
  85  2.3479774160285863271658E1L,
  86  1.9444210022760132894510E2L,
  87  7.7952888181207260646090E2L,
  88  1.6911722418503949084863E3L,
  89  2.0307734695595183428202E3L,
  90  1.2695660352705325274404E3L,
  91  3.2242573199748645407652E2L,
  92 };
  93
  94 /* Coefficients for log(x) = z + z^3 P(z^2)/Q(z^2),
  95  * where z = 2(x-1)/(x+1)
  96  * 1/sqrt(2) <= x < sqrt(2)
  97  * Theoretical peak relative error = 6.16e-22
  98  */
  99
 100 static const long double R[4] = {
 101  1.9757429581415468984296E-3L,
 102 -7.1990767473014147232598E-1L,
 103  1.0777257190312272158094E1L,
 104 -3.5717684488096787370998E1L,
 105 };
 106 static const long double S[4] = {
 107 /* 1.00000000000000000000E0L,*/
 108 -2.6201045551331104417768E1L,
 109  1.9361891836232102174846E2L,
 110 -4.2861221385716144629696E2L,
 111 };
 112 /* log10(2) */
 113 #define L102A 0.3125L
 114 #define L102B -1.1470004336018804786261e-2L
 115 /* log10(e) */
 116 #define L10EA 0.5L
 117 #define L10EB -6.5705518096748172348871e-2L
 118
 119 #define SQRTH 0.70710678118654752440L
 120
 121 long double
 122 log10l(long double x)
 123 {
 124 long double y;
 125 volatile long double z;
 126 int e;
 127
 128 if( isnan(x) )
 129         return(x);
 130 /* Test for domain */
 131 if( x <= 0.0L )
 132         {
 133         if( x == 0.0L )
 134                 return (-1.0L / (x - x));
 135         else
 136                 return (x - x) / (x - x);
 137         }
 138 if( x == INFINITY )
 139         return(INFINITY);
 140 /* separate mantissa from exponent */
 141
 142 /* Note, frexp is used so that denormal numbers
 143  * will be handled properly.
 144  */
 145 x = frexpl( x, &e );
 146
 147
 148 /* logarithm using log(x) = z + z**3 P(z)/Q(z),
 149  * where z = 2(x-1)/x+1)
 150  */
 151 if( (e > 2) || (e < -2) )
 152 {
 153 if( x < SQRTH )
 154         { /* 2( 2x-1 )/( 2x+1 ) */
 155         e -= 1;
 156         z = x - 0.5L;
 157         y = 0.5L * z + 0.5L;
 158         }
 159 else
 160         { /*  2 (x-1)/(x+1)   */
 161         z = x - 0.5L;
 162         z -= 0.5L;
 163         y = 0.5L * x  + 0.5L;
 164         }
 165 x = z / y;
 166 z = x*x;
 167 y = x * ( z * __polevll( z, (void *)R, 3 ) / __p1evll( z, (void *)S, 3 ) );
 168 goto done;
 169 }
 170
 171
 172 /* logarithm using log(1+x) = x - .5x**2 + x**3 P(x)/Q(x) */
 173
 174 if( x < SQRTH )
 175         {
 176         e -= 1;
 177         x = ldexpl( x, 1 ) - 1.0L; /*  2x - 1  */
 178         }
 179 else
 180         {
 181         x = x - 1.0L;
 182         }
 183 z = x*x;
 184 y = x * ( z * __polevll( x, (void *)P, 6 ) / __p1evll( x, (void *)Q, 7 ) );
 185 y = y - ldexpl( z, -1 );   /* -0.5x^2 + ... */
 186
 187 done:
 188
 189 /* Multiply log of fraction by log10(e)
 190  * and base 2 exponent by log10(2).
 191  *
 192  * ***CAUTION***
 193  *
 194  * This sequence of operations is critical and it may
 195  * be horribly defeated by some compiler optimizers.
 196  */
 197 z = y * (L10EB);
 198 z += x * (L10EB);
 199 z += e * (L102B);
 200 z += y * (L10EA);
 201 z += x * (L10EA);
 202 z += e * (L102A);
 203
 204 return( z );
 205 }