3 * Relative error logarithm
4 * Natural logarithm of 1+x, 128-bit long double precision
10 * long double x, y, log1pl();
18 * Returns the base e (2.718...) logarithm of 1+x.
20 * The argument 1+x is separated into its exponent and fractional
21 * parts. If the exponent is between -1 and +1, the logarithm
22 * of the fraction is approximated by
24 * log(1+x) = x - 0.5 x^2 + x^3 P(x)/Q(x).
26 * Otherwise, setting z = 2(w-1)/(w+1),
28 * log(w) = z + z^3 P(z)/Q(z).
35 * arithmetic domain # trials peak rms
36 * IEEE -1, 8 100000 1.9e-34 4.3e-35
39 /* Copyright 2001 by Stephen L. Moshier */
42 #include "math_private.h"
44 /* Coefficients for log(1+x) = x - x^2 / 2 + x^3 P(x)/Q(x)
45 * 1/sqrt(2) <= 1+x < sqrt(2)
46 * Theoretical peak relative error = 5.3e-37,
47 * relative peak error spread = 2.3e-14
50 P12
= 1.538612243596254322971797716843006400388E-6L,
51 P11
= 4.998469661968096229986658302195402690910E-1L,
52 P10
= 2.321125933898420063925789532045674660756E1L
,
53 P9
= 4.114517881637811823002128927449878962058E2L
,
54 P8
= 3.824952356185897735160588078446136783779E3L
,
55 P7
= 2.128857716871515081352991964243375186031E4L
,
56 P6
= 7.594356839258970405033155585486712125861E4L
,
57 P5
= 1.797628303815655343403735250238293741397E5L
,
58 P4
= 2.854829159639697837788887080758954924001E5L
,
59 P3
= 3.007007295140399532324943111654767187848E5L
,
60 P2
= 2.014652742082537582487669938141683759923E5L
,
61 P1
= 7.771154681358524243729929227226708890930E4L
,
62 P0
= 1.313572404063446165910279910527789794488E4L
,
63 /* Q12 = 1.000000000000000000000000000000000000000E0L, */
64 Q11
= 4.839208193348159620282142911143429644326E1L
,
65 Q10
= 9.104928120962988414618126155557301584078E2L
,
66 Q9
= 9.147150349299596453976674231612674085381E3L
,
67 Q8
= 5.605842085972455027590989944010492125825E4L
,
68 Q7
= 2.248234257620569139969141618556349415120E5L
,
69 Q6
= 6.132189329546557743179177159925690841200E5L
,
70 Q5
= 1.158019977462989115839826904108208787040E6L
,
71 Q4
= 1.514882452993549494932585972882995548426E6L
,
72 Q3
= 1.347518538384329112529391120390701166528E6L
,
73 Q2
= 7.777690340007566932935753241556479363645E5L
,
74 Q1
= 2.626900195321832660448791748036714883242E5L
,
75 Q0
= 3.940717212190338497730839731583397586124E4L
;
77 /* Coefficients for log(x) = z + z^3 P(z^2)/Q(z^2),
78 * where z = 2(x-1)/(x+1)
79 * 1/sqrt(2) <= x < sqrt(2)
80 * Theoretical peak relative error = 1.1e-35,
81 * relative peak error spread 1.1e-9
84 R5
= -8.828896441624934385266096344596648080902E-1L,
85 R4
= 8.057002716646055371965756206836056074715E1L
,
86 R3
= -2.024301798136027039250415126250455056397E3L
,
87 R2
= 2.048819892795278657810231591630928516206E4L
,
88 R1
= -8.977257995689735303686582344659576526998E4L
,
89 R0
= 1.418134209872192732479751274970992665513E5L
,
90 /* S6 = 1.000000000000000000000000000000000000000E0L, */
91 S5
= -1.186359407982897997337150403816839480438E2L
,
92 S4
= 3.998526750980007367835804959888064681098E3L
,
93 S3
= -5.748542087379434595104154610899551484314E4L
,
94 S2
= 4.001557694070773974936904547424676279307E5L
,
95 S1
= -1.332535117259762928288745111081235577029E6L
,
96 S0
= 1.701761051846631278975701529965589676574E6L
;
99 static long double C1
= 6.93145751953125E-1L;
100 static long double C2
= 1.428606820309417232121458176568075500134E-6L;
102 static long double sqrth
= 0.7071067811865475244008443621048490392848L;
103 /* ln (2^16384 * (1 - 2^-113)) */
104 static long double maxlog
= 1.1356523406294143949491931077970764891253E4L
;
105 static long double big
= 2e4932L
;
106 static long double zero
= 0.0L;
109 /* Make sure these are prototyped. */
110 long double frexpl (long double, int *);
111 long double ldexpl (long double, int);
116 __log1pl (long double xm1
)
118 long double x
, y
, z
, r
, s
;
119 ieee854_long_double_shape_type u
;
123 /* Test for NaN or infinity input. */
125 ix
= u
.parts32
.w0
& 0x7fffffff;
126 if (ix
>= 0x7fff0000)
129 /* log1p(+- 0) = +- 0. */
130 if ((ix
== 0) && (u
.parts32
.w1
| u
.parts32
.w2
| u
.parts32
.w3
) == 0)
135 /* log1p(-1) = -inf */
139 return (-1.0L / zero
);
141 return (zero
/ zero
);
144 /* Separate mantissa from exponent. */
146 /* Use frexp used so that denormal numbers will be handled properly. */
149 /* Logarithm using log(x) = z + z^3 P(z^2)/Q(z^2),
150 where z = 2(x-1)/x+1). */
151 if ((e
> 2) || (e
< -2))
154 { /* 2( 2x-1 )/( 2x+1 ) */
160 { /* 2 (x-1)/(x+1) */
188 /* Logarithm using log(1+x) = x - .5x^2 + x^3 P(x)/Q(x). */
194 x
= 2.0L * x
- 1.0L; /* 2x - 1 */
206 r
= (((((((((((P12
* x
240 weak_alias (__log1pl
, log1pl
)
241 #ifdef NO_LONG_DOUBLE
242 strong_alias (__log1p
, __log1pl
)
243 weak_alias (__log1p
, log1pl
)