| Rich Felker | b69f695 | 2012-03-13 01:17:53 -0400 | [diff] [blame] | 1 | #include "libm.h" |
| 2 | |
| nsz | df8b3e5 | 2012-03-13 16:38:22 +0100 | [diff] [blame] | 3 | /* k is such that k*ln2 has minimal relative error and x - kln2 > log(DBL_MIN) */ |
| 4 | static const int k = 2043; |
| 5 | static const double kln2 = 0x1.62066151add8bp+10; |
| Rich Felker | b69f695 | 2012-03-13 01:17:53 -0400 | [diff] [blame] | 6 | |
| nsz | df8b3e5 | 2012-03-13 16:38:22 +0100 | [diff] [blame] | 7 | /* exp(x)/2 for x >= log(DBL_MAX), slightly better than 0.5*exp(x/2)*exp(x/2) */ |
| Rich Felker | b69f695 | 2012-03-13 01:17:53 -0400 | [diff] [blame] | 8 | double __expo2(double x) |
| 9 | { |
| 10 | double scale; |
| Rich Felker | b69f695 | 2012-03-13 01:17:53 -0400 | [diff] [blame] | 11 | |
| nsz | df8b3e5 | 2012-03-13 16:38:22 +0100 | [diff] [blame] | 12 | /* note that k is odd and scale*scale overflows */ |
| 13 | INSERT_WORDS(scale, (uint32_t)(0x3ff + k/2) << 20, 0); |
| 14 | /* exp(x - k ln2) * 2**(k-1) */ |
| Rich Felker | b69f695 | 2012-03-13 01:17:53 -0400 | [diff] [blame] | 15 | return exp(x - kln2) * scale * scale; |
| 16 | } |