Qucs-core
0.0.19
|
00001 /* 00002 * fspecial.cpp - special functions implementation 00003 * 00004 * Copyright (C) 2006, 2008 Stefan Jahn <stefan@lkcc.org> 00005 * 00006 * This is free software; you can redistribute it and/or modify 00007 * it under the terms of the GNU General Public License as published by 00008 * the Free Software Foundation; either version 2, or (at your option) 00009 * any later version. 00010 * 00011 * This software is distributed in the hope that it will be useful, 00012 * but WITHOUT ANY WARRANTY; without even the implied warranty of 00013 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 00014 * GNU General Public License for more details. 00015 * 00016 * You should have received a copy of the GNU General Public License 00017 * along with this package; see the file COPYING. If not, write to 00018 * the Free Software Foundation, Inc., 51 Franklin Street - Fifth Floor, 00019 * Boston, MA 02110-1301, USA. 00020 * 00021 * $Id$ 00022 * 00023 */ 00024 00025 #if HAVE_CONFIG_H 00026 # include <config.h> 00027 #endif 00028 00029 #include <stdio.h> 00030 #include <stdlib.h> 00031 #include <string.h> 00032 #include <cmath> 00033 #include <float.h> 00034 #include <algorithm> 00035 00036 #include "compat.h" 00037 #include "constants.h" 00038 #include "fspecial.h" 00039 00040 #include <limits> 00041 00042 using qucs::pi_over_2; 00043 using qucs::sqrt_pi; 00044 using qucs::sqrt2; 00045 00046 /* The function computes the complete elliptic integral of first kind 00047 K() and the second kind E() using the arithmetic-geometric mean 00048 algorithm (AGM) by Abramowitz and Stegun. */ 00049 void fspecial::ellip_ke (nr_double_t arg, nr_double_t &k, nr_double_t &e) { 00050 int iMax = 16; 00051 if (arg == 1.0) { 00052 k = std::numeric_limits<nr_double_t>::infinity(); 00053 e = 0; 00054 } 00055 else if (std::isinf (arg) && arg < 0) { 00056 k = 0; 00057 e = std::numeric_limits<nr_double_t>::infinity(); 00058 } 00059 else { 00060 nr_double_t a, b, c, f, s, fk = 1, fe = 1, t, da = arg; 00061 int i; 00062 if (arg < 0) { 00063 fk = 1 / sqrt (1 - arg); 00064 fe = sqrt (1 - arg); 00065 da = -arg / (1 - arg); 00066 } 00067 a = 1; 00068 b = sqrt (1 - da); 00069 c = sqrt (da); 00070 f = 0.5; 00071 s = f * c * c; 00072 for (i = 0; i < iMax; i++) { 00073 t = (a + b) / 2; 00074 c = (a - b) / 2; 00075 b = sqrt (a * b); 00076 a = t; 00077 f *= 2; 00078 s += f * c * c; 00079 if (c / a < std::numeric_limits<nr_double_t>::epsilon()) break; 00080 } 00081 if (i >= iMax) { 00082 k = 0; e = 0; 00083 } 00084 else { 00085 k = pi_over_2 / a; 00086 e = pi_over_2 * (1 - s) / a; 00087 if (arg < 0) { 00088 k *= fk; 00089 e *= fe; 00090 } 00091 } 00092 } 00093 } 00094 00095 const nr_double_t SN_ACC = 1e-5; // Accuracy of sn(x) is SN_ACC^2 00096 const nr_double_t K_ERR = 1e-8; // Accuracy of K(k) 00097 00098 // Computes Carlson's elliptic integral of the first kind. 00099 nr_double_t fspecial::ellip_rf (nr_double_t x, nr_double_t y, nr_double_t z) { 00100 nr_double_t al, av, dx, dy, dz, e2, e3; 00101 nr_double_t sx, sy, sz, xt, yt, zt; 00102 00103 // constants 00104 const nr_double_t c1 = 1.0 / 24.0; 00105 const nr_double_t c2 = 0.1; 00106 const nr_double_t c3 = 3.0 / 44.0; 00107 const nr_double_t c4 = 1.0 / 14.0; 00108 00109 xt = x; 00110 yt = y; 00111 zt = z; 00112 do { 00113 sx = sqrt (xt); 00114 sy = sqrt (yt); 00115 sz = sqrt (zt); 00116 al = sx * (sy + sz) + sy * sz; 00117 xt = 0.25 * (xt + al); 00118 yt = 0.25 * (yt + al); 00119 zt = 0.25 * (zt + al); 00120 av = (xt + yt + zt) / 3.0; 00121 dx = (av - xt) / av; 00122 dy = (av - yt) / av; 00123 dz = (av - zt) / av; 00124 } 00125 while (std::max (std::max (fabs (dx), fabs (dy)), fabs (dz)) > K_ERR); 00126 00127 e2 = dx * dy - dz * dz; 00128 e3 = dx * dy * dz; 00129 return (1 + (c1 * e2 - c2 - c3 * e3) * e2 + c4 * e3) / sqrt (av); 00130 } 00131 00132 // Compute the Jacobian elliptic functions sn (u,k), cn (u,k) and dn (u,k). 00133 nr_double_t fspecial::ellip_sncndn (nr_double_t u, nr_double_t k, 00134 nr_double_t& sn, nr_double_t& cn, nr_double_t& dn) { 00135 nr_double_t a, b, c, d; 00136 nr_double_t fn[14], en[14]; 00137 int l; 00138 bool bo; 00139 00140 d = 1 - k; 00141 if (k) { 00142 bo = (k < 0); 00143 if (bo) { 00144 k /= -1 / d; 00145 u *= (d = sqrt (d)); 00146 } 00147 a = 1; 00148 dn = 1; 00149 for (int i = 1; i < 14; i++) { 00150 l = i; 00151 fn[i] = a; 00152 en[i] = (k = sqrt (k)); 00153 c = 0.5 * (a + k); 00154 if (fabs (a - k) <= SN_ACC * a) 00155 break; 00156 k *= a; 00157 a = c; 00158 } 00159 u *= c; 00160 sn = sin (u); 00161 cn = cos (u); 00162 if (sn) { 00163 a = cn / sn; 00164 c *= a; 00165 for (int ii = l; ii > 0; ii--) { 00166 b = fn[ii]; 00167 a *= c; 00168 c *= dn; 00169 dn = (en[ii] + a) / (b + a); 00170 a = c / b; 00171 } 00172 a = 1 / sqrt (c * c + 1); 00173 sn = (sn >= 0 ? a : -a); 00174 cn = sn * c; 00175 } 00176 if (bo) { 00177 a = dn; 00178 dn = cn; 00179 cn = a; 00180 sn /= d; 00181 } 00182 } 00183 else { 00184 cn = 1 / cosh (u); 00185 dn = cn; 00186 sn = tanh (u); 00187 } 00188 return sn; 00189 } 00190 00191 /* data for a Chebyshev series over a given interval */ 00192 struct cheb_series_t { 00193 nr_double_t * c; /* coefficients */ 00194 int order; /* order of expansion */ 00195 nr_double_t a; /* lower interval point */ 00196 nr_double_t b; /* upper interval point */ 00197 }; 00198 typedef struct cheb_series_t cheb_series; 00199 00200 static nr_double_t cheb_eval (const cheb_series * cs, const nr_double_t x) { 00201 nr_double_t d = 0.0; 00202 nr_double_t dd = 0.0; 00203 nr_double_t y = (2.0 * x - cs->a - cs->b) / (cs->b - cs->a); 00204 nr_double_t y2 = 2.0 * y; 00205 for (int i = cs->order; i >= 1; i--) { 00206 nr_double_t t = d; 00207 d = y2 * d - dd + cs->c[i]; 00208 dd = t; 00209 } 00210 d = y * d - dd + 0.5 * cs->c[0]; 00211 return d; 00212 } 00213 00214 #if !defined (HAVE_ERF) || !defined (HAVE_ERFC) 00215 00216 /* Chebyshev fit for erfc ((t+1)/2), -1 < t < 1 */ 00217 static nr_double_t erfc_xlt1_data[20] = { 00218 1.06073416421769980345174155056, 00219 -0.42582445804381043569204735291, 00220 0.04955262679620434040357683080, 00221 0.00449293488768382749558001242, 00222 -0.00129194104658496953494224761, 00223 -0.00001836389292149396270416979, 00224 0.00002211114704099526291538556, 00225 -5.23337485234257134673693179020e-7, 00226 -2.78184788833537885382530989578e-7, 00227 1.41158092748813114560316684249e-8, 00228 2.72571296330561699984539141865e-9, 00229 -2.06343904872070629406401492476e-10, 00230 -2.14273991996785367924201401812e-11, 00231 2.22990255539358204580285098119e-12, 00232 1.36250074650698280575807934155e-13, 00233 -1.95144010922293091898995913038e-14, 00234 -6.85627169231704599442806370690e-16, 00235 1.44506492869699938239521607493e-16, 00236 2.45935306460536488037576200030e-18, 00237 -9.29599561220523396007359328540e-19 00238 }; 00239 static cheb_series erfc_xlt1_cs = { 00240 erfc_xlt1_data, 19, -1, 1 00241 }; 00242 00243 /* Chebyshev fit for erfc (x) * exp (x^2), 1 < x < 5, x = 2t + 3, -1 < t < 1 */ 00244 static nr_double_t erfc_x15_data[25] = { 00245 0.44045832024338111077637466616, 00246 -0.143958836762168335790826895326, 00247 0.044786499817939267247056666937, 00248 -0.013343124200271211203618353102, 00249 0.003824682739750469767692372556, 00250 -0.001058699227195126547306482530, 00251 0.000283859419210073742736310108, 00252 -0.000073906170662206760483959432, 00253 0.000018725312521489179015872934, 00254 -4.62530981164919445131297264430e-6, 00255 1.11558657244432857487884006422e-6, 00256 -2.63098662650834130067808832725e-7, 00257 6.07462122724551777372119408710e-8, 00258 -1.37460865539865444777251011793e-8, 00259 3.05157051905475145520096717210e-9, 00260 -6.65174789720310713757307724790e-10, 00261 1.42483346273207784489792999706e-10, 00262 -3.00141127395323902092018744545e-11, 00263 6.22171792645348091472914001250e-12, 00264 -1.26994639225668496876152836555e-12, 00265 2.55385883033257575402681845385e-13, 00266 -5.06258237507038698392265499770e-14, 00267 9.89705409478327321641264227110e-15, 00268 -1.90685978789192181051961024995e-15, 00269 3.50826648032737849245113757340e-16 00270 }; 00271 static cheb_series erfc_x15_cs = { 00272 erfc_x15_data, 24, -1, 1 00273 }; 00274 00275 /* Chebyshev fit for erfc(x) * exp(x^2), 00276 5 < x < 10, x = (5t + 15)/2, -1 < t < */ 00277 static nr_double_t erfc_x510_data[20] = { 00278 1.11684990123545698684297865808, 00279 0.003736240359381998520654927536, 00280 -0.000916623948045470238763619870, 00281 0.000199094325044940833965078819, 00282 -0.000040276384918650072591781859, 00283 7.76515264697061049477127605790e-6, 00284 -1.44464794206689070402099225301e-6, 00285 2.61311930343463958393485241947e-7, 00286 -4.61833026634844152345304095560e-8, 00287 8.00253111512943601598732144340e-9, 00288 -1.36291114862793031395712122089e-9, 00289 2.28570483090160869607683087722e-10, 00290 -3.78022521563251805044056974560e-11, 00291 6.17253683874528285729910462130e-12, 00292 -9.96019290955316888445830597430e-13, 00293 1.58953143706980770269506726000e-13, 00294 -2.51045971047162509999527428316e-14, 00295 3.92607828989125810013581287560e-15, 00296 -6.07970619384160374392535453420e-16, 00297 9.12600607264794717315507477670e-17 00298 }; 00299 static cheb_series erfc_x510_cs = { 00300 erfc_x510_data, 19, -1, 1 00301 }; 00302 00303 /* Estimates erfc (x) valid for 8 < x < 100, this is based on index 5725 00304 in Hart et al. */ 00305 static nr_double_t erfc8 (nr_double_t x) { 00306 static nr_double_t p[] = { 00307 2.97886562639399288862, 00308 7.409740605964741794425, 00309 6.1602098531096305440906, 00310 5.019049726784267463450058, 00311 1.275366644729965952479585264, 00312 0.5641895835477550741253201704 00313 }; 00314 static nr_double_t q[] = { 00315 3.3690752069827527677, 00316 9.608965327192787870698, 00317 17.08144074746600431571095, 00318 12.0489519278551290360340491, 00319 9.396034016235054150430579648, 00320 2.260528520767326969591866945, 00321 1.0 00322 }; 00323 nr_double_t n = 0.0, d = 0.0; 00324 int i; 00325 n = p[5]; 00326 for (i = 4; i >= 0; --i) n = x * n + p[i]; 00327 d = q[6]; 00328 for (i = 5; i >= 0; --i) d = x * d + q[i]; 00329 00330 return exp (-x * x) * (n / d); 00331 } 00332 00333 #endif /* !HAVE_ERF || !HAVE_ERFC */ 00334 00335 #ifndef HAVE_ERFC 00336 00337 nr_double_t fspecial::erfc (nr_double_t x) { 00338 const nr_double_t ax = fabs (x); 00339 nr_double_t val; 00340 00341 if (ax <= 1.0) { 00342 nr_double_t t = 2.0 * ax - 1.0; 00343 val = cheb_eval (&erfc_xlt1_cs, t); 00344 } 00345 else if (ax <= 5.0) { 00346 nr_double_t ex2 = exp (-x * x); 00347 nr_double_t t = 0.5 * (ax - 3.0); 00348 val = ex2 * cheb_eval (&erfc_x15_cs, t); 00349 } 00350 else if (ax < 10.0) { 00351 nr_double_t ex = exp(-x * x) / ax; 00352 nr_double_t t = (2.0 * ax - 15.0) / 5.0; 00353 val = ex * cheb_eval (&erfc_x510_cs, t); 00354 } 00355 else { 00356 val = erfc8 (ax); 00357 } 00358 return (x < 0.0) ? 2.0 - val : val; 00359 } 00360 #else 00361 nr_double_t fspecial::erfc (nr_double_t d) { 00362 return ::erfc (d); 00363 } 00364 #endif /* HAVE_ERFC */ 00365 00366 #ifndef HAVE_ERF 00367 00368 /* Abramowitz + Stegun, 7.1.5 */ 00369 static nr_double_t erfseries (nr_double_t x) { 00370 nr_double_t c = x; 00371 nr_double_t e = c; 00372 nr_double_t d; 00373 for (int k = 1; k < 30; ++k) { 00374 c *= -x * x / k; 00375 d = c / (2.0 * k + 1.0); 00376 e += d; 00377 } 00378 return 2.0 / sqrt_pi * e; 00379 } 00380 00381 nr_double_t fspecial::erf (nr_double_t x) { 00382 if (fabs (x) < 1.0) { 00383 return erfseries (x); 00384 } 00385 return 1.0 - erfc (x); 00386 } 00387 00388 #else 00389 nr_double_t fspecial::erf (nr_double_t d) { 00390 return ::erf (d); 00391 } 00392 #endif /* HAVE_ERF */ 00393 00394 // Inverse of the error function erf(). 00395 nr_double_t fspecial::erfinv (nr_double_t y) { 00396 nr_double_t x = 0.0; // The output 00397 nr_double_t z = 0.0; // Intermediate variable 00398 nr_double_t y0 = 0.7; // Central range variable 00399 00400 // Coefficients in rational approximations. 00401 nr_double_t a[4] = { 0.886226899, -1.645349621, 0.914624893, -0.140543331}; 00402 nr_double_t b[4] = {-2.118377725, 1.442710462, -0.329097515, 0.012229801}; 00403 nr_double_t c[4] = {-1.970840454, -1.624906493, 3.429567803, 1.641345311}; 00404 nr_double_t d[2] = { 3.543889200, 1.637067800}; 00405 00406 if (y < -1.0 || 1.0 < y) { 00407 x = log (-1.0); 00408 } 00409 else if (y == -1.0 || 1.0 == y) { 00410 x = -y * log(0.0); 00411 } 00412 else if (-1.0 < y && y < -y0) { 00413 z = sqrt(-log((1.0+y)/2.0)); 00414 x = -(((c[3]*z+c[2])*z+c[1])*z+c[0])/((d[1]*z+d[0])*z+1.0); 00415 } 00416 else { 00417 if (-y0 < y && y < y0) { 00418 z = y * y; 00419 x = y*(((a[3]*z+a[2])*z+a[1])*z+a[0]) / 00420 ((((b[3]*z+b[3])*z+b[1])*z+b[0])*z+1.0); 00421 } 00422 else if (y0 < y && y < 1.0) { 00423 z = sqrt(-log((1.0-y)/2.0)); 00424 x = (((c[3]*z+c[2])*z+c[1])*z+c[0])/((d[1]*z+d[0])*z+1.0); 00425 } 00426 00427 // Two steps of Newton-Raphson correction to full accuracy. 00428 x = x - (erf (x) - y) / (2.0 / sqrt_pi * exp (-x * x)); 00429 x = x - (erf (x) - y) / (2.0 / sqrt_pi * exp (-x * x)); 00430 } 00431 return x; 00432 } 00433 00434 static nr_double_t bi0_data[12] = { 00435 -.07660547252839144951, 00436 1.92733795399380827000, 00437 .22826445869203013390, 00438 .01304891466707290428, 00439 .00043442709008164874, 00440 .00000942265768600193, 00441 .00000014340062895106, 00442 .00000000161384906966, 00443 .00000000001396650044, 00444 .00000000000009579451, 00445 .00000000000000053339, 00446 .00000000000000000245 00447 }; 00448 static cheb_series bi0_cs = { 00449 bi0_data, 11, -1, 1 00450 }; 00451 00452 static nr_double_t ai0_data[21] = { 00453 .07575994494023796, 00454 .00759138081082334, 00455 .00041531313389237, 00456 .00001070076463439, 00457 -.00000790117997921, 00458 -.00000078261435014, 00459 .00000027838499429, 00460 .00000000825247260, 00461 -.00000001204463945, 00462 .00000000155964859, 00463 .00000000022925563, 00464 -.00000000011916228, 00465 .00000000001757854, 00466 .00000000000112822, 00467 -.00000000000114684, 00468 .00000000000027155, 00469 -.00000000000002415, 00470 -.00000000000000608, 00471 .00000000000000314, 00472 -.00000000000000071, 00473 .00000000000000007 00474 }; 00475 static cheb_series ai0_cs = { 00476 ai0_data, 20, -1, 1 00477 }; 00478 00479 static nr_double_t ai02_data[22] = { 00480 .05449041101410882, 00481 .00336911647825569, 00482 .00006889758346918, 00483 .00000289137052082, 00484 .00000020489185893, 00485 .00000002266668991, 00486 .00000000339623203, 00487 .00000000049406022, 00488 .00000000001188914, 00489 -.00000000003149915, 00490 -.00000000001321580, 00491 -.00000000000179419, 00492 .00000000000071801, 00493 .00000000000038529, 00494 .00000000000001539, 00495 -.00000000000004151, 00496 -.00000000000000954, 00497 .00000000000000382, 00498 .00000000000000176, 00499 -.00000000000000034, 00500 -.00000000000000027, 00501 .00000000000000003 00502 }; 00503 static cheb_series ai02_cs = { 00504 ai02_data, 21, -1, 1 00505 }; 00506 00507 // Modified Bessel function of order zero. 00508 nr_double_t fspecial::i0 (nr_double_t x) { 00509 nr_double_t y = fabs (x); 00510 nr_double_t val; 00511 00512 if (y < 2.0 * sqrt (std::numeric_limits<nr_double_t>::epsilon())) { 00513 val = 1.0; 00514 } 00515 else if (y <= 3.0) { 00516 val = 2.75 + cheb_eval (&bi0_cs, y * y / 4.5 - 1.0); 00517 } 00518 else if (y <= 8.0) { 00519 val = cheb_eval (&ai0_cs, (48.0 / y - 11.0) / 5.0); 00520 val = exp (y) * (0.375 + val) / sqrt (y); 00521 } 00522 else { 00523 val = cheb_eval (&ai02_cs, 16.0 / y - 1.0); 00524 val = exp (y) * (0.375 + val) / sqrt (y); 00525 } 00526 return val; 00527 } 00528 00529 // Lower tail quantile for the standard normal distribution function. 00530 nr_double_t fspecial::ltqnorm (nr_double_t x) { 00531 nr_double_t q, r, e, u, z = 0.0; 00532 static nr_double_t a[] = { 00533 -3.969683028665376e+01, 2.209460984245205e+02, 00534 -2.759285104469687e+02, 1.383577518672690e+02, 00535 -3.066479806614716e+01, 2.506628277459239e+00 }; 00536 static nr_double_t b[] = { 00537 -5.447609879822406e+01, 1.615858368580409e+02, 00538 -1.556989798598866e+02, 6.680131188771972e+01, 00539 -1.328068155288572e+01 }; 00540 static nr_double_t c[] = { 00541 -7.784894002430293e-03, -3.223964580411365e-01, 00542 -2.400758277161838e+00, -2.549732539343734e+00, 00543 4.374664141464968e+00, 2.938163982698783e+00 }; 00544 static nr_double_t d[] = { 00545 7.784695709041462e-03, 3.224671290700398e-01, 00546 2.445134137142996e+00, 3.754408661907416e+00 }; 00547 00548 // Define break-points. 00549 nr_double_t pl = 0.02425; 00550 nr_double_t ph = 1.0 - pl; 00551 00552 // Rational approximation for central region: 00553 if (pl <= x && x <= ph) { 00554 q = x - 0.5; 00555 r = q * q; 00556 z = (((((a[0]*r+a[1])*r+a[2])*r+a[3])*r+a[4])*r+a[5])*q/ 00557 (((((b[0]*r+b[1])*r+b[2])*r+b[3])*r+b[4])*r+1); 00558 } 00559 // Rational approximation for lower region: 00560 else if (0.0 < x && x < pl) { 00561 q = sqrt(-2*log(x)); 00562 z = (((((c[0]*q+c[1])*q+c[2])*q+c[3])*q+c[4])*q+c[5])/ 00563 ((((d[0]*q+d[1])*q+d[2])*q+d[3])*q+1); 00564 } 00565 // Rational approximation for upper region: 00566 else if (ph < x && x < 1.0) { 00567 q = sqrt(-2*log(1-x)); 00568 z = -(((((c[0]*q+c[1])*q+c[2])*q+c[3])*q+c[4])*q+c[5])/ 00569 ((((d[0]*q+d[1])*q+d[2])*q+d[3])*q+1); 00570 } 00571 // Case when X = 0: 00572 else if (x == 0.0) { 00573 z = -std::numeric_limits<nr_double_t>::infinity(); 00574 } 00575 // Case when X = 1: 00576 else if (x == 1.0) { 00577 z = +std::numeric_limits<nr_double_t>::infinity(); 00578 } 00579 // Cases when output will be NaN: 00580 else if (x < 0.0 || x > 1.0 || std::isnan (x)) { 00581 z = +std::numeric_limits<nr_double_t>::quiet_NaN(); 00582 } 00583 00584 // The relative error of the approximation has absolute value less 00585 // than 1.15e-9. One iteration of Halley's rational method (third 00586 // order) gives full machine precision. 00587 if (0.0 < x && x < 1.0) { 00588 e = 0.5 * erfc (-z / sqrt2) - x; // error 00589 u = e * sqrt2 * sqrt_pi * exp (z * z / 2); // f(z)/df(z) 00590 z = z - u / (1 + z * u / 2); // Halley's method 00591 } 00592 return z; 00593 } 00594 00595 // Inverse of the error function erfc(). 00596 nr_double_t fspecial::erfcinv (nr_double_t x) { 00597 return -ltqnorm (x / 2.0) / sqrt2; 00598 }