Computes the Riemann zeta function ζ(x) for some value of x This method may return: Double#NaN for x = 1 (would be complex infinity) NOTE: This method is not yet complete in terms of accuracy. For x < -0.5, the values
(double x)
| 314 | * @return ζ(x) |
| 315 | */ |
| 316 | public static double zeta(double x) |
| 317 | { |
| 318 | if(x == 1) |
| 319 | return Double.NaN; |
| 320 | if(x < 0 || abs(1-x) <= 0.2) |
| 321 | { |
| 322 | if(x <= 0.2 && x > -2.) |
| 323 | { |
| 324 | /* |
| 325 | * For this specific range we keep our own approximant, |
| 326 | * see below comment |
| 327 | */ |
| 328 | return hornerPolyR(zeta_p_special, x)/hornerPolyR(zeta_q_special, x); |
| 329 | } |
| 330 | /* |
| 331 | * http://dlmf.nist.gov/25.4#E2 |
| 332 | * |
| 333 | * Reflect zeta across 1 if negative to make it positive. |
| 334 | * |
| 335 | * Reflect zeta across 1 if it is just less than one so that it will |
| 336 | * be just more than 1 (much easier to compute) |
| 337 | * |
| 338 | * Reflect if just more than 1 to an area that is easier to approximate |
| 339 | */ |
| 340 | double otherPart = 2*pow(2*PI, x-1)*sin(PI/2*x); |
| 341 | if(x < 0) |
| 342 | return otherPart*exp(lnGamma(1-x)*log(zeta(1-x))); |
| 343 | else//log(zeta(1-x)) would have caused a NaN |
| 344 | return otherPart*gamma(1-x)*zeta(1-x); |
| 345 | } |
| 346 | if(x < 14) |
| 347 | return hornerPolyR(zeta_p_l14, x)/hornerPolyR(zeta_q_l14, x); |
| 348 | if(x < 50) |
| 349 | { |
| 350 | //use truncated form of http://dlmf.nist.gov/25.2#E3 |
| 351 | double mul = 1/(1-pow(2, 1-x)); |
| 352 | double sumP = 0; |
| 353 | double sumN = 0; |
| 354 | for(int i = 11; i >= 1; i-=2)//all odd values are positive |
| 355 | sumP += pow(i, -x); |
| 356 | for(int i = 10; i >= 1; i-=2)//all even values are negative |
| 357 | sumN -= pow(i, -x); |
| 358 | return mul*(sumP+sumN); |
| 359 | } |
| 360 | //else x>=50, 1 is so close we might as well use it |
| 361 | return 1; |
| 362 | } |
| 363 | |
| 364 | /** |
| 365 | * upper polynomial approximation of the zeta function between [-0.20, 0.5] |