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

     1  /* origin: OpenBSD /usr/src/lib/libm/src/ld80/s_log1pl.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   *      Relative error logarithm
    19   *      Natural logarithm of 1+x, long double precision
    20   *
    21   *
    22   * SYNOPSIS:
    23   *
    24   * long double x, y, log1pl();
    25   *
    26   * y = log1pl( x );
    27   *
    28   *
    29   * DESCRIPTION:
    30   *
    31   * Returns the base e (2.718...) logarithm of 1+x.
    32   *
    33   * The argument 1+x is separated into its exponent and fractional
    34   * parts.  If the exponent is between -1 and +1, the logarithm
    35   * of the fraction is approximated by
    36   *
    37   *     log(1+x) = x - 0.5 x^2 + x^3 P(x)/Q(x).
    38   *
    39   * Otherwise, setting  z = 2(x-1)/x+1),
    40   *
    41   *     log(x) = z + z^3 P(z)/Q(z).
    42   *
    43   *
    44   * ACCURACY:
    45   *
    46   *                      Relative error:
    47   * arithmetic   domain     # trials      peak         rms
    48   *    IEEE     -1.0, 9.0    100000      8.2e-20    2.5e-20
    49   */
    50  
    51  #include "libm.h"
    52  
    53  #if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
    54  long double log1pl(long double x)
    55  {
    56  	return log1p(x);
    57  }
    58  #elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384
    59  /* Coefficients for log(1+x) = x - x^2 / 2 + x^3 P(x)/Q(x)
    60   * 1/sqrt(2) <= x < sqrt(2)
    61   * Theoretical peak relative error = 2.32e-20
    62   */
    63  static const long double P[] = {
    64   4.5270000862445199635215E-5L,
    65   4.9854102823193375972212E-1L,
    66   6.5787325942061044846969E0L,
    67   2.9911919328553073277375E1L,
    68   6.0949667980987787057556E1L,
    69   5.7112963590585538103336E1L,
    70   2.0039553499201281259648E1L,
    71  };
    72  static const long double Q[] = {
    73  /* 1.0000000000000000000000E0,*/
    74   1.5062909083469192043167E1L,
    75   8.3047565967967209469434E1L,
    76   2.2176239823732856465394E2L,
    77   3.0909872225312059774938E2L,
    78   2.1642788614495947685003E2L,
    79   6.0118660497603843919306E1L,
    80  };
    81  
    82  /* Coefficients for log(x) = z + z^3 P(z^2)/Q(z^2),
    83   * where z = 2(x-1)/(x+1)
    84   * 1/sqrt(2) <= x < sqrt(2)
    85   * Theoretical peak relative error = 6.16e-22
    86   */
    87  static const long double R[4] = {
    88   1.9757429581415468984296E-3L,
    89  -7.1990767473014147232598E-1L,
    90   1.0777257190312272158094E1L,
    91  -3.5717684488096787370998E1L,
    92  };
    93  static const long double S[4] = {
    94  /* 1.00000000000000000000E0L,*/
    95  -2.6201045551331104417768E1L,
    96   1.9361891836232102174846E2L,
    97  -4.2861221385716144629696E2L,
    98  };
    99  static const long double C1 = 6.9314575195312500000000E-1L;
   100  static const long double C2 = 1.4286068203094172321215E-6L;
   101  
   102  #define SQRTH 0.70710678118654752440L
   103  
   104  long double log1pl(long double xm1)
   105  {
   106  	long double x, y, z;
   107  	int e;
   108  
   109  	if (isnan(xm1))
   110  		return xm1;
   111  	if (xm1 == INFINITY)
   112  		return xm1;
   113  	if (xm1 == 0.0)
   114  		return xm1;
   115  
   116  	x = xm1 + 1.0;
   117  
   118  	/* Test for domain errors.  */
   119  	if (x <= 0.0) {
   120  		if (x == 0.0)
   121  			return -1/(x*x); /* -inf with divbyzero */
   122  		return 0/0.0f; /* nan with invalid */
   123  	}
   124  
   125  	/* Separate mantissa from exponent.
   126  	   Use frexp so that denormal numbers will be handled properly.  */
   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  	if (e > 2 || e < -2) {
   132  		if (x < SQRTH) { /* 2(2x-1)/(2x+1) */
   133  			e -= 1;
   134  			z = x - 0.5;
   135  			y = 0.5 * z + 0.5;
   136  		} else { /*  2 (x-1)/(x+1)   */
   137  			z = x - 0.5;
   138  			z -= 0.5;
   139  			y = 0.5 * x  + 0.5;
   140  		}
   141  		x = z / y;
   142  		z = x*x;
   143  		z = x * (z * __polevll(z, R, 3) / __p1evll(z, S, 3));
   144  		z = z + e * C2;
   145  		z = z + x;
   146  		z = z + e * C1;
   147  		return z;
   148  	}
   149  
   150  	/* logarithm using log(1+x) = x - .5x**2 + x**3 P(x)/Q(x) */
   151  	if (x < SQRTH) {
   152  		e -= 1;
   153  		if (e != 0)
   154  			x = 2.0 * x - 1.0;
   155  		else
   156  			x = xm1;
   157  	} else {
   158  		if (e != 0)
   159  			x = x - 1.0;
   160  		else
   161  			x = xm1;
   162  	}
   163  	z = x*x;
   164  	y = x * (z * __polevll(x, P, 6) / __p1evll(x, Q, 6));
   165  	y = y + e * C2;
   166  	z = y - 0.5 * z;
   167  	z = z + x;
   168  	z = z + e * C1;
   169  	return z;
   170  }
   171  #elif LDBL_MANT_DIG == 113 && LDBL_MAX_EXP == 16384
   172  // TODO: broken implementation to make things compile
   173  long double log1pl(long double x)
   174  {
   175  	return log1p(x);
   176  }
   177  #endif