Alex Bikfalvi
SimStream Documentation
Psi.cpp
00001 #include "Headers.h" 00002 #include "Psi.h" 00003 #include "Chebyshev.h" 00004 #include "Zeta.h" 00005 #include "Exp.h" 00006 00007 static double psics_data[23] = { 00008 -.038057080835217922, 00009 .491415393029387130, 00010 -.056815747821244730, 00011 .008357821225914313, 00012 -.001333232857994342, 00013 .000220313287069308, 00014 -.000037040238178456, 00015 .000006283793654854, 00016 -.000001071263908506, 00017 .000000183128394654, 00018 -.000000031353509361, 00019 .000000005372808776, 00020 -.000000000921168141, 00021 .000000000157981265, 00022 -.000000000027098646, 00023 .000000000004648722, 00024 -.000000000000797527, 00025 .000000000000136827, 00026 -.000000000000023475, 00027 .000000000000004027, 00028 -.000000000000000691, 00029 .000000000000000118, 00030 -.000000000000000020 00031 }; 00032 static ChebSeries psi_cs = { 00033 psics_data, 00034 22, 00035 -1, 1, 00036 17 00037 }; 00038 00039 static double apsics_data[16] = { 00040 -.0204749044678185, 00041 -.0101801271534859, 00042 .0000559718725387, 00043 -.0000012917176570, 00044 .0000000572858606, 00045 -.0000000038213539, 00046 .0000000003397434, 00047 -.0000000000374838, 00048 .0000000000048990, 00049 -.0000000000007344, 00050 .0000000000001233, 00051 -.0000000000000228, 00052 .0000000000000045, 00053 -.0000000000000009, 00054 .0000000000000002, 00055 -.0000000000000000 00056 }; 00057 static ChebSeries apsi_cs = { 00058 apsics_data, 00059 15, 00060 -1, 1, 00061 9 00062 }; 00063 00064 #define PSI_TABLE_NMAX 100 00065 static double psi_table[PSI_TABLE_NMAX+1] = { 00066 0.0, /* Infinity */ /* psi(0) */ 00067 -M_EULER, /* psi(1) */ 00068 0.42278433509846713939348790992, /* ... */ 00069 0.92278433509846713939348790992, 00070 1.25611766843180047272682124325, 00071 1.50611766843180047272682124325, 00072 1.70611766843180047272682124325, 00073 1.87278433509846713939348790992, 00074 2.01564147795560999653634505277, 00075 2.14064147795560999653634505277, 00076 2.25175258906672110764745616389, 00077 2.35175258906672110764745616389, 00078 2.44266167997581201673836525479, 00079 2.52599501330914535007169858813, 00080 2.60291809023222227314862166505, 00081 2.67434666166079370172005023648, 00082 2.74101332832746036838671690315, 00083 2.80351332832746036838671690315, 00084 2.86233685773922507426906984432, 00085 2.91789241329478062982462539988, 00086 2.97052399224214905087725697883, 00087 3.02052399224214905087725697883, 00088 3.06814303986119666992487602645, 00089 3.11359758531574212447033057190, 00090 3.15707584618530734186163491973, 00091 3.1987425128519740085283015864, 00092 3.2387425128519740085283015864, 00093 3.2772040513135124700667631249, 00094 3.3142410883505495071038001619, 00095 3.3499553740648352213895144476, 00096 3.3844381326855248765619282407, 00097 3.4177714660188582098952615740, 00098 3.4500295305349872421533260902, 00099 3.4812795305349872421533260902, 00100 3.5115825608380175451836291205, 00101 3.5409943255438998981248055911, 00102 3.5695657541153284695533770196, 00103 3.5973435318931062473311547974, 00104 3.6243705589201332743581818244, 00105 3.6506863483938174848844976139, 00106 3.6763273740348431259101386396, 00107 3.7013273740348431259101386396, 00108 3.7257176179372821503003825420, 00109 3.7495271417468059598241920658, 00110 3.7727829557002943319172153216, 00111 3.7955102284275670591899425943, 00112 3.8177324506497892814121648166, 00113 3.8394715810845718901078169905, 00114 3.8607481768292527411716467777, 00115 3.8815815101625860745049801110, 00116 3.9019896734278921969539597029, 00117 3.9219896734278921969539597029, 00118 3.9415975165651470989147440166, 00119 3.9608282857959163296839747858, 00120 3.9796962103242182164764276160, 00121 3.9982147288427367349949461345, 00122 4.0163965470245549168131279527, 00123 4.0342536898816977739559850956, 00124 4.0517975495308205809735289552, 00125 4.0690389288411654085597358518, 00126 4.0859880813835382899156680552, 00127 4.1026547480502049565823347218, 00128 4.1190481906731557762544658694, 00129 4.1351772229312202923834981274, 00130 4.1510502388042361653993711433, 00131 4.1666752388042361653993711433, 00132 4.1820598541888515500147557587, 00133 4.1972113693403667015299072739, 00134 4.2121367424746950597388624977, 00135 4.2268426248276362362094507330, 00136 4.2413353784508246420065521823, 00137 4.2556210927365389277208378966, 00138 4.2697055997787924488475984600, 00139 4.2835944886676813377364873489, 00140 4.2972931188046676391063503626, 00141 4.3108066323181811526198638761, 00142 4.3241399656515144859531972094, 00143 4.3372978603883565912163551041, 00144 4.3502848733753695782293421171, 00145 4.3631053861958823987421626300, 00146 4.3757636140439836645649474401, 00147 4.3882636140439836645649474401, 00148 4.4006092930563293435772931191, 00149 4.4128044150075488557724150703, 00150 4.4248526077786331931218126607, 00151 4.4367573696833950978837174226, 00152 4.4485220755657480390601880108, 00153 4.4601499825424922251066996387, 00154 4.4716442354160554434975042364, 00155 4.4830078717796918071338678728, 00156 4.4942438268358715824147667492, 00157 4.5053549379469826935258778603, 00158 4.5163439489359936825368668713, 00159 4.5272135141533849868846929582, 00160 4.5379662023254279976373811303, 00161 4.5486045001977684231692960239, 00162 4.5591308159872421073798223397, 00163 4.5695474826539087740464890064, 00164 4.5798567610044242379640147796, 00165 4.5900608426370772991885045755, 00166 4.6001618527380874001986055856 00167 }; 00168 00169 00170 #define PSI_1_TABLE_NMAX 100 00171 static double psi_1_table[PSI_1_TABLE_NMAX+1] = { 00172 0.0, /* Infinity */ /* psi(1,0) */ 00173 M_PI*M_PI/6.0, /* psi(1,1) */ 00174 0.644934066848226436472415, /* ... */ 00175 0.394934066848226436472415, 00176 0.2838229557371153253613041, 00177 0.2213229557371153253613041, 00178 0.1813229557371153253613041, 00179 0.1535451779593375475835263, 00180 0.1331370146940314251345467, 00181 0.1175120146940314251345467, 00182 0.1051663356816857461222010, 00183 0.0951663356816857461222010, 00184 0.0869018728717683907503002, 00185 0.0799574284273239463058557, 00186 0.0740402686640103368384001, 00187 0.0689382278476838062261552, 00188 0.0644937834032393617817108, 00189 0.0605875334032393617817108, 00190 0.0571273257907826143768665, 00191 0.0540409060376961946237801, 00192 0.0512708229352031198315363, 00193 0.0487708229352031198315363, 00194 0.0465032492390579951149830, 00195 0.0444371335365786562720078, 00196 0.0425467743683366902984728, 00197 0.0408106632572255791873617, 00198 0.0392106632572255791873617, 00199 0.0377313733163971768204978, 00200 0.0363596312039143235969038, 00201 0.0350841209998326909438426, 00202 0.0338950603577399442137594, 00203 0.0327839492466288331026483, 00204 0.0317433665203020901265817, 00205 0.03076680402030209012658168, 00206 0.02984853037475571730748159, 00207 0.02898347847164153045627052, 00208 0.02816715194102928555831133, 00209 0.02739554700275768062003973, 00210 0.02666508681283803124093089, 00211 0.02597256603721476254286995, 00212 0.02531510384129102815759710, 00213 0.02469010384129102815759710, 00214 0.02409521984367056414807896, 00215 0.02352832641963428296894063, 00216 0.02298749353699501850166102, 00217 0.02247096461137518379091722, 00218 0.02197713745088135663042339, 00219 0.02150454765882086513703965, 00220 0.02105185413233829383780923, 00221 0.02061782635456051606003145, 00222 0.02020133322669712580597065, 00223 0.01980133322669712580597065, 00224 0.01941686571420193164987683, 00225 0.01904704322899483105816086, 00226 0.01869104465298913508094477, 00227 0.01834810912486842177504628, 00228 0.01801753061247172756017024, 00229 0.01769865306145131939690494, 00230 0.01739086605006319997554452, 00231 0.01709360088954001329302371, 00232 0.01680632711763538818529605, 00233 0.01652854933985761040751827, 00234 0.01625980437882562975715546, 00235 0.01599965869724394401313881, 00236 0.01574770606433893015574400, 00237 0.01550356543933893015574400, 00238 0.01526687904880638577704578, 00239 0.01503731063741979257227076, 00240 0.01481454387422086185273411, 00241 0.01459828089844231513993134, 00242 0.01438824099085987447620523, 00243 0.01418415935820681325171544, 00244 0.01398578601958352422176106, 00245 0.01379288478501562298719316, 00246 0.01360523231738567365335942, 00247 0.01342261726990576130858221, 00248 0.01324483949212798353080444, 00249 0.01307170929822216635628920, 00250 0.01290304679189732236910755, 00251 0.01273868124291638877278934, 00252 0.01257845051066194236996928, 00253 0.01242220051066194236996928, 00254 0.01226978472038606978956995, 00255 0.01212106372098095378719041, 00256 0.01197590477193174490346273, 00257 0.01183418141592267460867815, 00258 0.01169577311142440471248438, 00259 0.01156056489076458859566448, 00260 0.01142844704164317229232189, 00261 0.01129931481023821361463594, 00262 0.01117306812421372175754719, 00263 0.01104961133409026496742374, 00264 0.01092885297157366069257770, 00265 0.01081070552355853781923177, 00266 0.01069508522063334415522437, 00267 0.01058191183901270133041676, 00268 0.01047110851491297833872701, 00269 0.01036260157046853389428257, 00270 0.01025632035036012704977199, /* ... */ 00271 0.01015219706839427948625679, /* psi(1,99) */ 00272 0.01005016666333357139524567 /* psi(1,100) */ 00273 }; 00274 00275 00276 /* digamma for x both positive and negative; we do both 00277 * cases here because of the way we use even/odd parts 00278 * of the function 00279 */ 00280 static int PsiX(const double x, MathResult * result) 00281 { 00282 const double y = fabs(x); 00283 00284 if(x == 0.0 || x == -1.0 || x == -2.0) 00285 { 00286 DOMAIN_ERROR(result); 00287 } 00288 else if(y >= 2.0) 00289 { 00290 const double t = 8.0/(y*y)-1.0; 00291 MathResult result_c; 00292 ChebEval(&apsi_cs, t, &result_c); 00293 if(x < 0.0) 00294 { 00295 const double s = sin(M_PI*x); 00296 const double c = cos(M_PI*x); 00297 if(fabs(s) < 2.0*MATH_SQRT_DBL_MIN) 00298 { 00299 DOMAIN_ERROR(result); 00300 } 00301 else 00302 { 00303 result->val = log(y) - 0.5/x + result_c.val - M_PI * c/s; 00304 result->err = M_PI*fabs(x)*MATH_DBL_EPSILON/(s*s); 00305 result->err += result_c.err; 00306 result->err += MATH_DBL_EPSILON * fabs(result->val); 00307 return MATH_SUCCESS; 00308 } 00309 } 00310 else 00311 { 00312 result->val = log(y) - 0.5/x + result_c.val; 00313 result->err = result_c.err; 00314 result->err += MATH_DBL_EPSILON * fabs(result->val); 00315 return MATH_SUCCESS; 00316 } 00317 } 00318 else 00319 { /* -2 < x < 2 */ 00320 MathResult result_c; 00321 00322 if(x < -1.0) 00323 { /* x = -2 + v */ 00324 const double v = x + 2.0; 00325 const double t1 = 1.0/x; 00326 const double t2 = 1.0/(x+1.0); 00327 const double t3 = 1.0/v; 00328 ChebEval(&psi_cs, 2.0*v-1.0, &result_c); 00329 00330 result->val = -(t1 + t2 + t3) + result_c.val; 00331 result->err = MATH_DBL_EPSILON * (fabs(t1) + fabs(x/(t2*t2)) + fabs(x/(t3*t3))); 00332 result->err += result_c.err; 00333 result->err += MATH_DBL_EPSILON * fabs(result->val); 00334 return MATH_SUCCESS; 00335 } 00336 else if(x < 0.0) 00337 { /* x = -1 + v */ 00338 const double v = x + 1.0; 00339 const double t1 = 1.0/x; 00340 const double t2 = 1.0/v; 00341 ChebEval(&psi_cs, 2.0*v-1.0, &result_c); 00342 00343 result->val = -(t1 + t2) + result_c.val; 00344 result->err = MATH_DBL_EPSILON * (fabs(t1) + fabs(x/(t2*t2))); 00345 result->err += result_c.err; 00346 result->err += MATH_DBL_EPSILON * fabs(result->val); 00347 return MATH_SUCCESS; 00348 } 00349 else if(x < 1.0) 00350 { /* x = v */ 00351 const double t1 = 1.0/x; 00352 ChebEval(&psi_cs, 2.0*x-1.0, &result_c); 00353 00354 result->val = -t1 + result_c.val; 00355 result->err = MATH_DBL_EPSILON * t1; 00356 result->err += result_c.err; 00357 result->err += MATH_DBL_EPSILON * fabs(result->val); 00358 return MATH_SUCCESS; 00359 } 00360 else 00361 { /* x = 1 + v */ 00362 const double v = x - 1.0; 00363 return ChebEval(&psi_cs, 2.0*v-1.0, result); 00364 } 00365 } 00366 } 00367 00368 /* generic polygamma; assumes n >= 0 and x > 0 00369 */ 00370 static int PsiNXg0(const int n, const double x, MathResult * result) 00371 { 00372 if(n == 0) 00373 { 00374 return Psi(x, result); 00375 } 00376 else 00377 { 00378 /* Abramowitz + Stegun 6.4.10 */ 00379 MathResult ln_nf; 00380 MathResult hzeta; 00381 00382 int stat_hz = HZeta(n+1.0, x, &hzeta); 00383 int stat_nf = LnFact((unsigned int) n, &ln_nf); 00384 int stat_e = ExpMultErr(ln_nf.val, ln_nf.err, hzeta.val, hzeta.err, result); 00385 00386 if(MATH_IS_EVEN(n)) result->val = -result->val; 00387 00388 return MATH_ERROR_SELECT_3(stat_e, stat_nf, stat_hz); 00389 } 00390 } 00391 00392 int PsiInt(const int n, MathResult * result) 00393 { 00394 /* CHECK_POINTER(result) */ 00395 00396 if(n <= 0) 00397 { 00398 DOMAIN_ERROR(result); 00399 } 00400 else if(n <= PSI_TABLE_NMAX) 00401 { 00402 result->val = psi_table[n]; 00403 result->err = MATH_DBL_EPSILON * fabs(result->val); 00404 return MATH_SUCCESS; 00405 } 00406 else 00407 { 00408 /* Abramowitz+Stegun 6.3.18 */ 00409 const double c2 = -1.0/12.0; 00410 const double c3 = 1.0/120.0; 00411 const double c4 = -1.0/252.0; 00412 const double c5 = 1.0/240.0; 00413 const double ni2 = (1.0/n)*(1.0/n); 00414 const double ser = ni2 * (c2 + ni2 * (c3 + ni2 * (c4 + ni2*c5))); 00415 result->val = log((double)n) - 0.5/n + ser; 00416 result->err = MATH_DBL_EPSILON * (fabs(log((double)n)) + fabs(0.5/n) + fabs(ser)); 00417 result->err += MATH_DBL_EPSILON * fabs(result->val); 00418 return MATH_SUCCESS; 00419 } 00420 } 00421 00422 int Psi(const double x, MathResult * result) 00423 { 00424 /* CHECK_POINTER(result) */ 00425 return PsiX(x, result); 00426 } 00427 00428 00429 int Psi1Int(const int n, MathResult * result) 00430 { 00431 /* CHECK_POINTER(result) */ 00432 if(n <= 0) 00433 { 00434 DOMAIN_ERROR(result); 00435 } 00436 else if(n <= PSI_1_TABLE_NMAX) 00437 { 00438 result->val = psi_1_table[n]; 00439 result->err = MATH_DBL_EPSILON * result->val; 00440 return MATH_SUCCESS; 00441 } 00442 else 00443 { 00444 /* Abramowitz+Stegun 6.4.12 00445 * double-precision for n > 100 00446 */ 00447 const double c0 = -1.0/30.0; 00448 const double c1 = 1.0/42.0; 00449 const double c2 = -1.0/30.0; 00450 const double ni2 = (1.0/n)*(1.0/n); 00451 const double ser = ni2*ni2 * (c0 + ni2*(c1 + c2*ni2)); 00452 result->val = (1.0 + 0.5/n + 1.0/(6.0*n*n) + ser) / n; 00453 result->err = MATH_DBL_EPSILON * result->val; 00454 return MATH_SUCCESS; 00455 } 00456 } 00457 00458 00459 int Psi1(const double x, MathResult * result) 00460 { 00461 /* CHECK_POINTER(result) */ 00462 00463 if(x == 0.0 || x == -1.0 || x == -2.0) 00464 { 00465 DOMAIN_ERROR(result); 00466 } 00467 else if(x > 0.0) 00468 { 00469 return PsiNXg0(1, x, result); 00470 } 00471 else if(x > -5.0) 00472 { 00473 /* Abramowitz + Stegun 6.4.6 */ 00474 int M = (int)-floor(x); 00475 double fx = x + M; 00476 double sum = 0.0; 00477 int m; 00478 00479 if(fx == 0.0) DOMAIN_ERROR(result); 00480 00481 for(m = 0; m < M; ++m) 00482 sum += 1.0/((x+m)*(x+m)); 00483 00484 { 00485 int stat_psi = PsiNXg0(1, fx, result); 00486 result->val += sum; 00487 result->err += M * MATH_DBL_EPSILON * sum; 00488 return stat_psi; 00489 } 00490 } 00491 else 00492 { 00493 /* Abramowitz + Stegun 6.4.7 */ 00494 const double sin_px = sin(M_PI * x); 00495 const double d = M_PI*M_PI/(sin_px*sin_px); 00496 MathResult r; 00497 int stat_psi = PsiNXg0(1, 1.0-x, &r); 00498 result->val = d - r.val; 00499 result->err = r.err + 2.0*MATH_DBL_EPSILON*d; 00500 return stat_psi; 00501 } 00502 } 00503 00504 int PsiN(const int n, const double x, MathResult * result) 00505 { 00506 /* CHECK_POINTER(result) */ 00507 00508 if(n == 0) 00509 { 00510 return Psi(x, result); 00511 } 00512 else if(n == 1) 00513 { 00514 return Psi1(x, result); 00515 } 00516 else if(n < 0 || x <= 0.0) 00517 { 00518 DOMAIN_ERROR(result); 00519 } 00520 else 00521 { 00522 MathResult ln_nf; 00523 MathResult hzeta; 00524 int stat_hz = HZeta(n+1.0, x, &hzeta); 00525 int stat_nf = LnFact((unsigned int) n, &ln_nf); 00526 int stat_e = ExpMultErr(ln_nf.val, ln_nf.err, hzeta.val, hzeta.err, result); 00527 00528 if(MATH_IS_EVEN(n)) result->val = -result->val; 00529 00530 return MATH_ERROR_SELECT_3(stat_e, stat_nf, stat_hz); 00531 } 00532 }
Last updated: February 8, 2011