1 | /* log10q.c |
---|---|

2 | * |

3 | * Common logarithm, 128-bit __float128 precision |

4 | * |

5 | * |

6 | * |

7 | * SYNOPSIS: |

8 | * |

9 | * __float128 x, y, log10l(); |

10 | * |

11 | * y = log10q( x ); |

12 | * |

13 | * |

14 | * |

15 | * DESCRIPTION: |

16 | * |

17 | * Returns the base 10 logarithm of x. |

18 | * |

19 | * The argument is separated into its exponent and fractional |

20 | * parts. If the exponent is between -1 and +1, the logarithm |

21 | * of the fraction is approximated by |

22 | * |

23 | * log(1+x) = x - 0.5 x^2 + x^3 P(x)/Q(x). |

24 | * |

25 | * Otherwise, setting z = 2(x-1)/x+1), |

26 | * |

27 | * log(x) = z + z^3 P(z)/Q(z). |

28 | * |

29 | * |

30 | * |

31 | * ACCURACY: |

32 | * |

33 | * Relative error: |

34 | * arithmetic domain # trials peak rms |

35 | * IEEE 0.5, 2.0 30000 2.3e-34 4.9e-35 |

36 | * IEEE exp(+-10000) 30000 1.0e-34 4.1e-35 |

37 | * |

38 | * In the tests over the interval exp(+-10000), the logarithms |

39 | * of the random arguments were uniformly distributed over |

40 | * [-10000, +10000]. |

41 | * |

42 | */ |

43 | |

44 | /* |

45 | Cephes Math Library Release 2.2: January, 1991 |

46 | Copyright 1984, 1991 by Stephen L. Moshier |

47 | Adapted for glibc November, 2001 |

48 | |

49 | This library is free software; you can redistribute it and/or |

50 | modify it under the terms of the GNU Lesser General Public |

51 | License as published by the Free Software Foundation; either |

52 | version 2.1 of the License, or (at your option) any later version. |

53 | |

54 | This library is distributed in the hope that it will be useful, |

55 | but WITHOUT ANY WARRANTY; without even the implied warranty of |

56 | MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU |

57 | Lesser General Public License for more details. |

58 | |

59 | You should have received a copy of the GNU Lesser General Public |

60 | License along with this library; if not, write to the Free Software |

61 | Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA |

62 | |

63 | */ |

64 | |

65 | #include "quadmath-imp.h" |

66 | |

67 | /* Coefficients for ln(1+x) = x - x**2/2 + x**3 P(x)/Q(x) |

68 | * 1/sqrt(2) <= x < sqrt(2) |

69 | * Theoretical peak relative error = 5.3e-37, |

70 | * relative peak error spread = 2.3e-14 |

71 | */ |

72 | static const __float128 P[13] = |

73 | { |

74 | 1.313572404063446165910279910527789794488E4Q, |

75 | 7.771154681358524243729929227226708890930E4Q, |

76 | 2.014652742082537582487669938141683759923E5Q, |

77 | 3.007007295140399532324943111654767187848E5Q, |

78 | 2.854829159639697837788887080758954924001E5Q, |

79 | 1.797628303815655343403735250238293741397E5Q, |

80 | 7.594356839258970405033155585486712125861E4Q, |

81 | 2.128857716871515081352991964243375186031E4Q, |

82 | 3.824952356185897735160588078446136783779E3Q, |

83 | 4.114517881637811823002128927449878962058E2Q, |

84 | 2.321125933898420063925789532045674660756E1Q, |

85 | 4.998469661968096229986658302195402690910E-1Q, |

86 | 1.538612243596254322971797716843006400388E-6Q |

87 | }; |

88 | static const __float128 Q[12] = |

89 | { |

90 | 3.940717212190338497730839731583397586124E4Q, |

91 | 2.626900195321832660448791748036714883242E5Q, |

92 | 7.777690340007566932935753241556479363645E5Q, |

93 | 1.347518538384329112529391120390701166528E6Q, |

94 | 1.514882452993549494932585972882995548426E6Q, |

95 | 1.158019977462989115839826904108208787040E6Q, |

96 | 6.132189329546557743179177159925690841200E5Q, |

97 | 2.248234257620569139969141618556349415120E5Q, |

98 | 5.605842085972455027590989944010492125825E4Q, |

99 | 9.147150349299596453976674231612674085381E3Q, |

100 | 9.104928120962988414618126155557301584078E2Q, |

101 | 4.839208193348159620282142911143429644326E1Q |

102 | /* 1.000000000000000000000000000000000000000E0Q, */ |

103 | }; |

104 | |

105 | /* Coefficients for log(x) = z + z^3 P(z^2)/Q(z^2), |

106 | * where z = 2(x-1)/(x+1) |

107 | * 1/sqrt(2) <= x < sqrt(2) |

108 | * Theoretical peak relative error = 1.1e-35, |

109 | * relative peak error spread 1.1e-9 |

110 | */ |

111 | static const __float128 R[6] = |

112 | { |

113 | 1.418134209872192732479751274970992665513E5Q, |

114 | -8.977257995689735303686582344659576526998E4Q, |

115 | 2.048819892795278657810231591630928516206E4Q, |

116 | -2.024301798136027039250415126250455056397E3Q, |

117 | 8.057002716646055371965756206836056074715E1Q, |

118 | -8.828896441624934385266096344596648080902E-1Q |

119 | }; |

120 | static const __float128 S[6] = |

121 | { |

122 | 1.701761051846631278975701529965589676574E6Q, |

123 | -1.332535117259762928288745111081235577029E6Q, |

124 | 4.001557694070773974936904547424676279307E5Q, |

125 | -5.748542087379434595104154610899551484314E4Q, |

126 | 3.998526750980007367835804959888064681098E3Q, |

127 | -1.186359407982897997337150403816839480438E2Q |

128 | /* 1.000000000000000000000000000000000000000E0Q, */ |

129 | }; |

130 | |

131 | static const __float128 |

132 | /* log10(2) */ |

133 | L102A = 0.3125Q, |

134 | L102B = -1.14700043360188047862611052755069732318101185E-2Q, |

135 | /* log10(e) */ |

136 | L10EA = 0.5Q, |

137 | L10EB = -6.570551809674817234887108108339491770560299E-2Q, |

138 | /* sqrt(2)/2 */ |

139 | SQRTH = 7.071067811865475244008443621048490392848359E-1Q; |

140 | |

141 | |

142 | |

143 | /* Evaluate P[n] x^n + P[n-1] x^(n-1) + ... + P[0] */ |

144 | |

145 | static __float128 |

146 | neval (__float128 x, const __float128 *p, int n) |

147 | { |

148 | __float128 y; |

149 | |

150 | p += n; |

151 | y = *p--; |

152 | do |

153 | { |

154 | y = y * x + *p--; |

155 | } |

156 | while (--n > 0); |

157 | return y; |

158 | } |

159 | |

160 | |

161 | /* Evaluate x^n+1 + P[n] x^(n) + P[n-1] x^(n-1) + ... + P[0] */ |

162 | |

163 | static __float128 |

164 | deval (__float128 x, const __float128 *p, int n) |

165 | { |

166 | __float128 y; |

167 | |

168 | p += n; |

169 | y = x + *p--; |

170 | do |

171 | { |

172 | y = y * x + *p--; |

173 | } |

174 | while (--n > 0); |

175 | return y; |

176 | } |

177 | |

178 | |

179 | |

180 | __float128 |

181 | log10q (__float128 x) |

182 | { |

183 | __float128 z; |

184 | __float128 y; |

185 | int e; |

186 | int64_t hx, lx; |

187 | |

188 | /* Test for domain */ |

189 | GET_FLT128_WORDS64 (hx, lx, x); |

190 | if (((hx & 0x7fffffffffffffffLL) | lx) == 0) |

191 | return (-1.0Q / fabsq (x)); /* log10l(+-0)=-inf */ |

192 | if (hx < 0) |

193 | return (x - x) / (x - x); |

194 | if (hx >= 0x7fff000000000000LL) |

195 | return (x + x); |

196 | |

197 | if (x == 1.0Q) |

198 | return 0.0Q; |

199 | |

200 | /* separate mantissa from exponent */ |

201 | |

202 | /* Note, frexp is used so that denormal numbers |

203 | * will be handled properly. |

204 | */ |

205 | x = frexpq (x, &e); |

206 | |

207 | |

208 | /* logarithm using log(x) = z + z**3 P(z)/Q(z), |

209 | * where z = 2(x-1)/x+1) |

210 | */ |

211 | if ((e > 2) || (e < -2)) |

212 | { |

213 | if (x < SQRTH) |

214 | { /* 2( 2x-1 )/( 2x+1 ) */ |

215 | e -= 1; |

216 | z = x - 0.5Q; |

217 | y = 0.5Q * z + 0.5Q; |

218 | } |

219 | else |

220 | { /* 2 (x-1)/(x+1) */ |

221 | z = x - 0.5Q; |

222 | z -= 0.5Q; |

223 | y = 0.5Q * x + 0.5Q; |

224 | } |

225 | x = z / y; |

226 | z = x * x; |

227 | y = x * (z * neval (z, R, 5) / deval (z, S, 5)); |

228 | goto done; |

229 | } |

230 | |

231 | |

232 | /* logarithm using log(1+x) = x - .5x**2 + x**3 P(x)/Q(x) */ |

233 | |

234 | if (x < SQRTH) |

235 | { |

236 | e -= 1; |

237 | x = 2.0 * x - 1.0Q; /* 2x - 1 */ |

238 | } |

239 | else |

240 | { |

241 | x = x - 1.0Q; |

242 | } |

243 | z = x * x; |

244 | y = x * (z * neval (x, P, 12) / deval (x, Q, 11)); |

245 | y = y - 0.5 * z; |

246 | |

247 | done: |

248 | |

249 | /* Multiply log of fraction by log10(e) |

250 | * and base 2 exponent by log10(2). |

251 | */ |

252 | z = y * L10EB; |

253 | z += x * L10EB; |

254 | z += e * L102B; |

255 | z += y * L10EA; |

256 | z += x * L10EA; |

257 | z += e * L102A; |

258 | return (z); |

259 | } |

260 |