123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191 |
- #include <openlibm_math.h>
- #include "math_private.h"
- static long double P[] = {
- 4.5270000862445199635215E-5L,
- 4.9854102823193375972212E-1L,
- 6.5787325942061044846969E0L,
- 2.9911919328553073277375E1L,
- 6.0949667980987787057556E1L,
- 5.7112963590585538103336E1L,
- 2.0039553499201281259648E1L,
- };
- static long double Q[] = {
- 1.5062909083469192043167E1L,
- 8.3047565967967209469434E1L,
- 2.2176239823732856465394E2L,
- 3.0909872225312059774938E2L,
- 2.1642788614495947685003E2L,
- 6.0118660497603843919306E1L,
- };
- static long double R[4] = {
- 1.9757429581415468984296E-3L,
- -7.1990767473014147232598E-1L,
- 1.0777257190312272158094E1L,
- -3.5717684488096787370998E1L,
- };
- static long double S[4] = {
- -2.6201045551331104417768E1L,
- 1.9361891836232102174846E2L,
- -4.2861221385716144629696E2L,
- };
- static const long double C1 = 6.9314575195312500000000E-1L;
- static const long double C2 = 1.4286068203094172321215E-6L;
- #define SQRTH 0.70710678118654752440L
- long double
- log1pl(long double xm1)
- {
- long double x, y, z;
- int e;
- if( isnan(xm1) )
- return(xm1);
- if( xm1 == INFINITY )
- return(xm1);
- if(xm1 == 0.0)
- return(xm1);
- x = xm1 + 1.0L;
- if( x <= 0.0L )
- {
- if( x == 0.0L )
- return( -INFINITY );
- else
- return( NAN );
- }
- x = frexpl( x, &e );
- if( (e > 2) || (e < -2) )
- {
- if( x < SQRTH )
- {
- e -= 1;
- z = x - 0.5L;
- y = 0.5L * z + 0.5L;
- }
- else
- {
- z = x - 0.5L;
- z -= 0.5L;
- y = 0.5L * x + 0.5L;
- }
- x = z / y;
- z = x*x;
- z = x * ( z * __polevll( z, R, 3 ) / __p1evll( z, S, 3 ) );
- z = z + e * C2;
- z = z + x;
- z = z + e * C1;
- return( z );
- }
- if( x < SQRTH )
- {
- e -= 1;
- if (e != 0)
- x = 2.0 * x - 1.0L;
- else
- x = xm1;
- }
- else
- {
- if (e != 0)
- x = x - 1.0L;
- else
- x = xm1;
- }
- z = x*x;
- y = x * ( z * __polevll( x, P, 6 ) / __p1evll( x, Q, 6 ) );
- y = y + e * C2;
- z = y - 0.5 * z;
- z = z + x;
- z = z + e * C1;
- return( z );
- }
|