Qucs-core  0.0.19
fspecial.cpp
Go to the documentation of this file.
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 }