| OLD | NEW |
| 1 /* origin: OpenBSD /usr/src/lib/libm/src/ld80/e_tgammal.c */ | 1 /* origin: OpenBSD /usr/src/lib/libm/src/ld80/e_tgammal.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 33 matching lines...) Expand 10 before | Expand all | Expand 10 after Loading... |
| 44 * IEEE -40,+40 10000 3.6e-19 7.9e-20 | 44 * IEEE -40,+40 10000 3.6e-19 7.9e-20 |
| 45 * IEEE -1755,+1755 10000 4.8e-18 6.5e-19 | 45 * IEEE -1755,+1755 10000 4.8e-18 6.5e-19 |
| 46 * | 46 * |
| 47 * Accuracy for large arguments is dominated by error in powl(). | 47 * Accuracy for large arguments is dominated by error in powl(). |
| 48 * | 48 * |
| 49 */ | 49 */ |
| 50 | 50 |
| 51 #include "libm.h" | 51 #include "libm.h" |
| 52 | 52 |
| 53 #if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024 | 53 #if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024 |
| 54 long double tgammal(long double x) | 54 long double tgammal(long double x) { |
| 55 { | 55 return tgamma(x); |
| 56 » return tgamma(x); | |
| 57 } | 56 } |
| 58 #elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384 | 57 #elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384 |
| 59 /* | 58 /* |
| 60 tgamma(x+2) = tgamma(x+2) P(x)/Q(x) | 59 tgamma(x+2) = tgamma(x+2) P(x)/Q(x) |
| 61 0 <= x <= 1 | 60 0 <= x <= 1 |
| 62 Relative error | 61 Relative error |
| 63 n=7, d=8 | 62 n=7, d=8 |
| 64 Peak error = 1.83e-20 | 63 Peak error = 1.83e-20 |
| 65 Relative error spread = 8.4e-23 | 64 Relative error spread = 8.4e-23 |
| 66 */ | 65 */ |
| 67 static const long double P[8] = { | 66 static const long double P[8] = { |
| 68 4.212760487471622013093E-5L, | 67 4.212760487471622013093E-5L, 4.542931960608009155600E-4L, |
| 69 4.542931960608009155600E-4L, | 68 4.092666828394035500949E-3L, 2.385363243461108252554E-2L, |
| 70 4.092666828394035500949E-3L, | 69 1.113062816019361559013E-1L, 3.629515436640239168939E-1L, |
| 71 2.385363243461108252554E-2L, | 70 8.378004301573126728826E-1L, 1.000000000000000000009E0L, |
| 72 1.113062816019361559013E-1L, | |
| 73 3.629515436640239168939E-1L, | |
| 74 8.378004301573126728826E-1L, | |
| 75 1.000000000000000000009E0L, | |
| 76 }; | 71 }; |
| 77 static const long double Q[9] = { | 72 static const long double Q[9] = { |
| 78 -1.397148517476170440917E-5L, | 73 -1.397148517476170440917E-5L, 2.346584059160635244282E-4L, |
| 79 2.346584059160635244282E-4L, | 74 -1.237799246653152231188E-3L, -7.955933682494738320586E-4L, |
| 80 -1.237799246653152231188E-3L, | 75 2.773706565840072979165E-2L, -4.633887671244534213831E-2L, |
| 81 -7.955933682494738320586E-4L, | 76 -2.243510905670329164562E-1L, 4.150160950588455434583E-1L, |
| 82 2.773706565840072979165E-2L, | 77 9.999999999999999999908E-1L, |
| 83 -4.633887671244534213831E-2L, | |
| 84 -2.243510905670329164562E-1L, | |
| 85 4.150160950588455434583E-1L, | |
| 86 9.999999999999999999908E-1L, | |
| 87 }; | 78 }; |
| 88 | 79 |
| 89 /* | 80 /* |
| 90 static const long double P[] = { | 81 static const long double P[] = { |
| 91 -3.01525602666895735709e0L, | 82 -3.01525602666895735709e0L, |
| 92 -3.25157411956062339893e1L, | 83 -3.25157411956062339893e1L, |
| 93 -2.92929976820724030353e2L, | 84 -2.92929976820724030353e2L, |
| 94 -1.70730828800510297666e3L, | 85 -1.70730828800510297666e3L, |
| 95 -7.96667499622741999770e3L, | 86 -7.96667499622741999770e3L, |
| 96 -2.59780216007146401957e4L, | 87 -2.59780216007146401957e4L, |
| (...skipping 18 matching lines...) Expand all Loading... |
| 115 /* Stirling's formula for the gamma function | 106 /* Stirling's formula for the gamma function |
| 116 tgamma(x) = sqrt(2 pi) x^(x-.5) exp(-x) (1 + 1/x P(1/x)) | 107 tgamma(x) = sqrt(2 pi) x^(x-.5) exp(-x) (1 + 1/x P(1/x)) |
| 117 z(x) = x | 108 z(x) = x |
| 118 13 <= x <= 1024 | 109 13 <= x <= 1024 |
| 119 Relative error | 110 Relative error |
| 120 n=8, d=0 | 111 n=8, d=0 |
| 121 Peak error = 9.44e-21 | 112 Peak error = 9.44e-21 |
| 122 Relative error spread = 8.8e-4 | 113 Relative error spread = 8.8e-4 |
| 123 */ | 114 */ |
| 124 static const long double STIR[9] = { | 115 static const long double STIR[9] = { |
| 125 7.147391378143610789273E-4L, | 116 7.147391378143610789273E-4L, -2.363848809501759061727E-5L, |
| 126 -2.363848809501759061727E-5L, | 117 -5.950237554056330156018E-4L, 6.989332260623193171870E-5L, |
| 127 -5.950237554056330156018E-4L, | 118 7.840334842744753003862E-4L, -2.294719747873185405699E-4L, |
| 128 6.989332260623193171870E-5L, | 119 -2.681327161876304418288E-3L, 3.472222222230075327854E-3L, |
| 129 7.840334842744753003862E-4L, | 120 8.333333333333331800504E-2L, |
| 130 -2.294719747873185405699E-4L, | |
| 131 -2.681327161876304418288E-3L, | |
| 132 3.472222222230075327854E-3L, | |
| 133 8.333333333333331800504E-2L, | |
| 134 }; | 121 }; |
| 135 | 122 |
| 136 #define MAXSTIR 1024.0L | 123 #define MAXSTIR 1024.0L |
| 137 static const long double SQTPI = 2.50662827463100050242E0L; | 124 static const long double SQTPI = 2.50662827463100050242E0L; |
| 138 | 125 |
| 139 /* 1/tgamma(x) = z P(z) | 126 /* 1/tgamma(x) = z P(z) |
| 140 * z(x) = 1/x | 127 * z(x) = 1/x |
| 141 * 0 < x < 0.03125 | 128 * 0 < x < 0.03125 |
| 142 * Peak relative error 4.2e-23 | 129 * Peak relative error 4.2e-23 |
| 143 */ | 130 */ |
| 144 static const long double S[9] = { | 131 static const long double S[9] = { |
| 145 -1.193945051381510095614E-3L, | 132 -1.193945051381510095614E-3L, 7.220599478036909672331E-3L, |
| 146 7.220599478036909672331E-3L, | 133 -9.622023360406271645744E-3L, -4.219773360705915470089E-2L, |
| 147 -9.622023360406271645744E-3L, | 134 1.665386113720805206758E-1L, -4.200263503403344054473E-2L, |
| 148 -4.219773360705915470089E-2L, | 135 -6.558780715202540684668E-1L, 5.772156649015328608253E-1L, |
| 149 1.665386113720805206758E-1L, | 136 1.000000000000000000000E0L, |
| 150 -4.200263503403344054473E-2L, | |
| 151 -6.558780715202540684668E-1L, | |
| 152 5.772156649015328608253E-1L, | |
| 153 1.000000000000000000000E0L, | |
| 154 }; | 137 }; |
| 155 | 138 |
| 156 /* 1/tgamma(-x) = z P(z) | 139 /* 1/tgamma(-x) = z P(z) |
| 157 * z(x) = 1/x | 140 * z(x) = 1/x |
| 158 * 0 < x < 0.03125 | 141 * 0 < x < 0.03125 |
| 159 * Peak relative error 5.16e-23 | 142 * Peak relative error 5.16e-23 |
| 160 * Relative error spread = 2.5e-24 | 143 * Relative error spread = 2.5e-24 |
| 161 */ | 144 */ |
| 162 static const long double SN[9] = { | 145 static const long double SN[9] = { |
| 163 1.133374167243894382010E-3L, | 146 1.133374167243894382010E-3L, 7.220837261893170325704E-3L, |
| 164 7.220837261893170325704E-3L, | 147 9.621911155035976733706E-3L, -4.219773343731191721664E-2L, |
| 165 9.621911155035976733706E-3L, | 148 -1.665386113944413519335E-1L, -4.200263503402112910504E-2L, |
| 166 -4.219773343731191721664E-2L, | 149 6.558780715202536547116E-1L, 5.772156649015328608727E-1L, |
| 167 -1.665386113944413519335E-1L, | 150 -1.000000000000000000000E0L, |
| 168 -4.200263503402112910504E-2L, | |
| 169 6.558780715202536547116E-1L, | |
| 170 5.772156649015328608727E-1L, | |
| 171 -1.000000000000000000000E0L, | |
| 172 }; | 151 }; |
| 173 | 152 |
| 174 static const long double PIL = 3.1415926535897932384626L; | 153 static const long double PIL = 3.1415926535897932384626L; |
| 175 | 154 |
| 176 /* Gamma function computed by Stirling's formula. | 155 /* Gamma function computed by Stirling's formula. |
| 177 */ | 156 */ |
| 178 static long double stirf(long double x) | 157 static long double stirf(long double x) { |
| 179 { | 158 long double y, w, v; |
| 180 » long double y, w, v; | |
| 181 | 159 |
| 182 » w = 1.0/x; | 160 w = 1.0 / x; |
| 183 » /* For large x, use rational coefficients from the analytical expansion.
*/ | 161 /* For large x, use rational coefficients from the analytical expansion. */ |
| 184 » if (x > 1024.0) | 162 if (x > 1024.0) |
| 185 » » w = (((((6.97281375836585777429E-5L * w | 163 w = (((((6.97281375836585777429E-5L * w + 7.84039221720066627474E-4L) * w - |
| 186 » » + 7.84039221720066627474E-4L) * w | 164 2.29472093621399176955E-4L) * |
| 187 » » - 2.29472093621399176955E-4L) * w | 165 w - |
| 188 » » - 2.68132716049382716049E-3L) * w | 166 2.68132716049382716049E-3L) * |
| 189 » » + 3.47222222222222222222E-3L) * w | 167 w + |
| 190 » » + 8.33333333333333333333E-2L) * w | 168 3.47222222222222222222E-3L) * |
| 191 » » + 1.0; | 169 w + |
| 192 » else | 170 8.33333333333333333333E-2L) * |
| 193 » » w = 1.0 + w * __polevll(w, STIR, 8); | 171 w + |
| 194 » y = expl(x); | 172 1.0; |
| 195 » if (x > MAXSTIR) { /* Avoid overflow in pow() */ | 173 else |
| 196 » » v = powl(x, 0.5L * x - 0.25L); | 174 w = 1.0 + w * __polevll(w, STIR, 8); |
| 197 » » y = v * (v / y); | 175 y = expl(x); |
| 198 » } else { | 176 if (x > MAXSTIR) { /* Avoid overflow in pow() */ |
| 199 » » y = powl(x, x - 0.5L) / y; | 177 v = powl(x, 0.5L * x - 0.25L); |
| 200 » } | 178 y = v * (v / y); |
| 201 » y = SQTPI * y * w; | 179 } else { |
| 202 » return y; | 180 y = powl(x, x - 0.5L) / y; |
| 181 } |
| 182 y = SQTPI * y * w; |
| 183 return y; |
| 203 } | 184 } |
| 204 | 185 |
| 205 long double tgammal(long double x) | 186 long double tgammal(long double x) { |
| 206 { | 187 long double p, q, z; |
| 207 » long double p, q, z; | |
| 208 | 188 |
| 209 » if (!isfinite(x)) | 189 if (!isfinite(x)) |
| 210 » » return x + INFINITY; | 190 return x + INFINITY; |
| 211 | 191 |
| 212 » q = fabsl(x); | 192 q = fabsl(x); |
| 213 » if (q > 13.0) { | 193 if (q > 13.0) { |
| 214 » » if (x < 0.0) { | 194 if (x < 0.0) { |
| 215 » » » p = floorl(q); | 195 p = floorl(q); |
| 216 » » » z = q - p; | 196 z = q - p; |
| 217 » » » if (z == 0) | 197 if (z == 0) |
| 218 » » » » return 0 / z; | 198 return 0 / z; |
| 219 » » » if (q > MAXGAML) { | 199 if (q > MAXGAML) { |
| 220 » » » » z = 0; | 200 z = 0; |
| 221 » » » } else { | 201 } else { |
| 222 » » » » if (z > 0.5) { | 202 if (z > 0.5) { |
| 223 » » » » » p += 1.0; | 203 p += 1.0; |
| 224 » » » » » z = q - p; | 204 z = q - p; |
| 225 » » » » } | 205 } |
| 226 » » » » z = q * sinl(PIL * z); | 206 z = q * sinl(PIL * z); |
| 227 » » » » z = fabsl(z) * stirf(q); | 207 z = fabsl(z) * stirf(q); |
| 228 » » » » z = PIL/z; | 208 z = PIL / z; |
| 229 » » » } | 209 } |
| 230 » » » if (0.5 * p == floorl(q * 0.5)) | 210 if (0.5 * p == floorl(q * 0.5)) |
| 231 » » » » z = -z; | 211 z = -z; |
| 232 » » } else if (x > MAXGAML) { | 212 } else if (x > MAXGAML) { |
| 233 » » » z = x * 0x1p16383L; | 213 z = x * 0x1p16383L; |
| 234 » » } else { | 214 } else { |
| 235 » » » z = stirf(x); | 215 z = stirf(x); |
| 236 » » } | 216 } |
| 237 » » return z; | 217 return z; |
| 238 » } | 218 } |
| 239 | 219 |
| 240 » z = 1.0; | 220 z = 1.0; |
| 241 » while (x >= 3.0) { | 221 while (x >= 3.0) { |
| 242 » » x -= 1.0; | 222 x -= 1.0; |
| 243 » » z *= x; | 223 z *= x; |
| 244 » } | 224 } |
| 245 » while (x < -0.03125L) { | 225 while (x < -0.03125L) { |
| 246 » » z /= x; | 226 z /= x; |
| 247 » » x += 1.0; | 227 x += 1.0; |
| 248 » } | 228 } |
| 249 » if (x <= 0.03125L) | 229 if (x <= 0.03125L) |
| 250 » » goto small; | 230 goto small; |
| 251 » while (x < 2.0) { | 231 while (x < 2.0) { |
| 252 » » z /= x; | 232 z /= x; |
| 253 » » x += 1.0; | 233 x += 1.0; |
| 254 » } | 234 } |
| 255 » if (x == 2.0) | 235 if (x == 2.0) |
| 256 » » return z; | 236 return z; |
| 257 | 237 |
| 258 » x -= 2.0; | 238 x -= 2.0; |
| 259 » p = __polevll(x, P, 7); | 239 p = __polevll(x, P, 7); |
| 260 » q = __polevll(x, Q, 8); | 240 q = __polevll(x, Q, 8); |
| 261 » z = z * p / q; | 241 z = z * p / q; |
| 262 » return z; | 242 return z; |
| 263 | 243 |
| 264 small: | 244 small: |
| 265 » /* z==1 if x was originally +-0 */ | 245 /* z==1 if x was originally +-0 */ |
| 266 » if (x == 0 && z != 1) | 246 if (x == 0 && z != 1) |
| 267 » » return x / x; | 247 return x / x; |
| 268 » if (x < 0.0) { | 248 if (x < 0.0) { |
| 269 » » x = -x; | 249 x = -x; |
| 270 » » q = z / (x * __polevll(x, SN, 8)); | 250 q = z / (x * __polevll(x, SN, 8)); |
| 271 » } else | 251 } else |
| 272 » » q = z / (x * __polevll(x, S, 8)); | 252 q = z / (x * __polevll(x, S, 8)); |
| 273 » return q; | 253 return q; |
| 274 } | 254 } |
| 275 #elif LDBL_MANT_DIG == 113 && LDBL_MAX_EXP == 16384 | 255 #elif LDBL_MANT_DIG == 113 && LDBL_MAX_EXP == 16384 |
| 276 // TODO: broken implementation to make things compile | 256 // TODO: broken implementation to make things compile |
| 277 long double tgammal(long double x) | 257 long double tgammal(long double x) { |
| 278 { | 258 return tgamma(x); |
| 279 » return tgamma(x); | |
| 280 } | 259 } |
| 281 #endif | 260 #endif |
| OLD | NEW |