1 /* $OpenBSD: s_expm1l.c,v 1.2 2011/07/20 21:02:51 martynas Exp $ */ 2 3 /* 4 * Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net> 5 * 6 * Permission to use, copy, modify, and distribute this software for any 7 * purpose with or without fee is hereby granted, provided that the above 8 * copyright notice and this permission notice appear in all copies. 9 * 10 * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES 11 * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF 12 * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR 13 * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES 14 * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN 15 * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF 16 * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE. 17 */ 18 19 /* expm1l.c 20 * 21 * Exponential function, minus 1 22 * Long double precision 23 * 24 * 25 * SYNOPSIS: 26 * 27 * long double x, y, expm1l(); 28 * 29 * y = expm1l( x ); 30 * 31 * 32 * 33 * DESCRIPTION: 34 * 35 * Returns e (2.71828...) raised to the x power, minus 1. 36 * 37 * Range reduction is accomplished by separating the argument 38 * into an integer k and fraction f such that 39 * 40 * x k f 41 * e = 2 e. 42 * 43 * An expansion x + .5 x^2 + x^3 R(x) approximates exp(f) - 1 44 * in the basic range [-0.5 ln 2, 0.5 ln 2]. 45 * 46 * 47 * ACCURACY: 48 * 49 * Relative error: 50 * arithmetic domain # trials peak rms 51 * IEEE -45,+MAXLOG 200,000 1.2e-19 2.5e-20 52 * 53 * ERROR MESSAGES: 54 * 55 * message condition value returned 56 * expm1l overflow x > MAXLOG MAXNUM 57 * 58 */ 59 60 #include <math.h> 61 62 static const long double MAXLOGL = 1.1356523406294143949492E4L; 63 64 /* exp(x) - 1 = x + 0.5 x^2 + x^3 P(x)/Q(x) 65 -.5 ln 2 < x < .5 ln 2 66 Theoretical peak relative error = 3.4e-22 */ 67 68 static const long double 69 P0 = -1.586135578666346600772998894928250240826E4L, 70 P1 = 2.642771505685952966904660652518429479531E3L, 71 P2 = -3.423199068835684263987132888286791620673E2L, 72 P3 = 1.800826371455042224581246202420972737840E1L, 73 P4 = -5.238523121205561042771939008061958820811E-1L, 74 75 Q0 = -9.516813471998079611319047060563358064497E4L, 76 Q1 = 3.964866271411091674556850458227710004570E4L, 77 Q2 = -7.207678383830091850230366618190187434796E3L, 78 Q3 = 7.206038318724600171970199625081491823079E2L, 79 Q4 = -4.002027679107076077238836622982900945173E1L, 80 /* Q5 = 1.000000000000000000000000000000000000000E0 */ 81 82 /* C1 + C2 = ln 2 */ 83 C1 = 6.93145751953125E-1L, 84 C2 = 1.428606820309417232121458176568075500134E-6L, 85 /* ln 2^-65 */ 86 minarg = -4.5054566736396445112120088E1L; 87 static const long double huge = 0x1p10000L; 88 89 long double 90 expm1l(long double x) 91 { 92 long double px, qx, xx; 93 int k; 94 95 /* Overflow. */ 96 if (x > MAXLOGL) 97 return (huge*huge); /* overflow */ 98 99 if (x == 0.0) 100 return x; 101 102 /* Minimum value. */ 103 if (x < minarg) 104 return -1.0L; 105 106 xx = C1 + C2; 107 108 /* Express x = ln 2 (k + remainder), remainder not exceeding 1/2. */ 109 px = floorl (0.5 + x / xx); 110 k = px; 111 /* remainder times ln 2 */ 112 x -= px * C1; 113 x -= px * C2; 114 115 /* Approximate exp(remainder ln 2). */ 116 px = (((( P4 * x 117 + P3) * x 118 + P2) * x 119 + P1) * x 120 + P0) * x; 121 122 qx = (((( x 123 + Q4) * x 124 + Q3) * x 125 + Q2) * x 126 + Q1) * x 127 + Q0; 128 129 xx = x * x; 130 qx = x + (0.5 * xx + xx * px / qx); 131 132 /* exp(x) = exp(k ln 2) exp(remainder ln 2) = 2^k exp(remainder ln 2). 133 We have qx = exp(remainder ln 2) - 1, so 134 exp(x) - 1 = 2^k (qx + 1) - 1 = 2^k qx + 2^k - 1. */ 135 px = ldexpl(1.0L, k); 136 x = px * qx + (px - 1.0); 137 return x; 138 } 139