github.com/afumu/libc@v0.0.6/musl/src/math/logl.c (about)

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