OLD | NEW |
1 /* origin: OpenBSD /usr/src/lib/libm/src/ld80/e_expm1l.c */ | 1 /* origin: OpenBSD /usr/src/lib/libm/src/ld80/e_expm1l.c */ |
2 /* | 2 /* |
3 * Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net> | 3 * Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net> |
4 * | 4 * |
5 * Permission to use, copy, modify, and distribute this software for any | 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 | 6 * purpose with or without fee is hereby granted, provided that the above |
7 * copyright notice and this permission notice appear in all copies. | 7 * copyright notice and this permission notice appear in all copies. |
8 * | 8 * |
9 * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES | 9 * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES |
10 * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF | 10 * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF |
(...skipping 32 matching lines...) Expand 10 before | Expand all | Expand 10 after Loading... |
43 * ACCURACY: | 43 * ACCURACY: |
44 * | 44 * |
45 * Relative error: | 45 * Relative error: |
46 * arithmetic domain # trials peak rms | 46 * arithmetic domain # trials peak rms |
47 * IEEE -45,+maxarg 200,000 1.2e-19 2.5e-20 | 47 * IEEE -45,+maxarg 200,000 1.2e-19 2.5e-20 |
48 */ | 48 */ |
49 | 49 |
50 #include "libm.h" | 50 #include "libm.h" |
51 | 51 |
52 #if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024 | 52 #if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024 |
53 long double expm1l(long double x) | 53 long double expm1l(long double x) { |
54 { | 54 return expm1(x); |
55 » return expm1(x); | |
56 } | 55 } |
57 #elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384 | 56 #elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384 |
58 | 57 |
59 /* exp(x) - 1 = x + 0.5 x^2 + x^3 P(x)/Q(x) | 58 /* exp(x) - 1 = x + 0.5 x^2 + x^3 P(x)/Q(x) |
60 -.5 ln 2 < x < .5 ln 2 | 59 -.5 ln 2 < x < .5 ln 2 |
61 Theoretical peak relative error = 3.4e-22 */ | 60 Theoretical peak relative error = 3.4e-22 */ |
62 static const long double | 61 static const long double P0 = -1.586135578666346600772998894928250240826E4L, |
63 P0 = -1.586135578666346600772998894928250240826E4L, | 62 P1 = 2.642771505685952966904660652518429479531E3L, |
64 P1 = 2.642771505685952966904660652518429479531E3L, | 63 P2 = -3.423199068835684263987132888286791620673E2L, |
65 P2 = -3.423199068835684263987132888286791620673E2L, | 64 P3 = 1.800826371455042224581246202420972737840E1L, |
66 P3 = 1.800826371455042224581246202420972737840E1L, | 65 P4 = -5.238523121205561042771939008061958820811E-1L, |
67 P4 = -5.238523121205561042771939008061958820811E-1L, | 66 Q0 = -9.516813471998079611319047060563358064497E4L, |
68 Q0 = -9.516813471998079611319047060563358064497E4L, | 67 Q1 = 3.964866271411091674556850458227710004570E4L, |
69 Q1 = 3.964866271411091674556850458227710004570E4L, | 68 Q2 = -7.207678383830091850230366618190187434796E3L, |
70 Q2 = -7.207678383830091850230366618190187434796E3L, | 69 Q3 = 7.206038318724600171970199625081491823079E2L, |
71 Q3 = 7.206038318724600171970199625081491823079E2L, | 70 Q4 = -4.002027679107076077238836622982900945173E1L, |
72 Q4 = -4.002027679107076077238836622982900945173E1L, | 71 /* Q5 = 1.000000000000000000000000000000000000000E0 */ |
73 /* Q5 = 1.000000000000000000000000000000000000000E0 */ | 72 /* C1 + C2 = ln 2 */ |
74 /* C1 + C2 = ln 2 */ | 73 C1 = 6.93145751953125E-1L, |
75 C1 = 6.93145751953125E-1L, | 74 C2 = 1.428606820309417232121458176568075500134E-6L, |
76 C2 = 1.428606820309417232121458176568075500134E-6L, | 75 /* ln 2^-65 */ |
77 /* ln 2^-65 */ | 76 minarg = -4.5054566736396445112120088E1L, |
78 minarg = -4.5054566736396445112120088E1L, | 77 /* ln 2^16384 */ |
79 /* ln 2^16384 */ | 78 maxarg = 1.1356523406294143949492E4L; |
80 maxarg = 1.1356523406294143949492E4L; | |
81 | 79 |
82 long double expm1l(long double x) | 80 long double expm1l(long double x) { |
83 { | 81 long double px, qx, xx; |
84 » long double px, qx, xx; | 82 int k; |
85 » int k; | |
86 | 83 |
87 » if (isnan(x)) | 84 if (isnan(x)) |
88 » » return x; | 85 return x; |
89 » if (x > maxarg) | 86 if (x > maxarg) |
90 » » return x*0x1p16383L; /* overflow, unless x==inf */ | 87 return x * 0x1p16383L; /* overflow, unless x==inf */ |
91 » if (x == 0.0) | 88 if (x == 0.0) |
92 » » return x; | 89 return x; |
93 » if (x < minarg) | 90 if (x < minarg) |
94 » » return -1.0; | 91 return -1.0; |
95 | 92 |
96 » xx = C1 + C2; | 93 xx = C1 + C2; |
97 » /* Express x = ln 2 (k + remainder), remainder not exceeding 1/2. */ | 94 /* Express x = ln 2 (k + remainder), remainder not exceeding 1/2. */ |
98 » px = floorl(0.5 + x / xx); | 95 px = floorl(0.5 + x / xx); |
99 » k = px; | 96 k = px; |
100 » /* remainder times ln 2 */ | 97 /* remainder times ln 2 */ |
101 » x -= px * C1; | 98 x -= px * C1; |
102 » x -= px * C2; | 99 x -= px * C2; |
103 | 100 |
104 » /* Approximate exp(remainder ln 2).*/ | 101 /* Approximate exp(remainder ln 2).*/ |
105 » px = (((( P4 * x + P3) * x + P2) * x + P1) * x + P0) * x; | 102 px = ((((P4 * x + P3) * x + P2) * x + P1) * x + P0) * x; |
106 » qx = (((( x + Q4) * x + Q3) * x + Q2) * x + Q1) * x + Q0; | 103 qx = ((((x + Q4) * x + Q3) * x + Q2) * x + Q1) * x + Q0; |
107 » xx = x * x; | 104 xx = x * x; |
108 » qx = x + (0.5 * xx + xx * px / qx); | 105 qx = x + (0.5 * xx + xx * px / qx); |
109 | 106 |
110 » /* exp(x) = exp(k ln 2) exp(remainder ln 2) = 2^k exp(remainder ln 2). | 107 /* exp(x) = exp(k ln 2) exp(remainder ln 2) = 2^k exp(remainder ln 2). |
111 » We have qx = exp(remainder ln 2) - 1, so | 108 We have qx = exp(remainder ln 2) - 1, so |
112 » exp(x) - 1 = 2^k (qx + 1) - 1 = 2^k qx + 2^k - 1. */ | 109 exp(x) - 1 = 2^k (qx + 1) - 1 = 2^k qx + 2^k - 1. */ |
113 » px = scalbnl(1.0, k); | 110 px = scalbnl(1.0, k); |
114 » x = px * qx + (px - 1.0); | 111 x = px * qx + (px - 1.0); |
115 » return x; | 112 return x; |
116 } | 113 } |
117 #elif LDBL_MANT_DIG == 113 && LDBL_MAX_EXP == 16384 | 114 #elif LDBL_MANT_DIG == 113 && LDBL_MAX_EXP == 16384 |
118 // TODO: broken implementation to make things compile | 115 // TODO: broken implementation to make things compile |
119 long double expm1l(long double x) | 116 long double expm1l(long double x) { |
120 { | 117 return expm1(x); |
121 » return expm1(x); | |
122 } | 118 } |
123 #endif | 119 #endif |
OLD | NEW |