[ VIGRA Homepage | Function Index | Class Index | Namespaces | File List | Main Page ]
vigra/bessel.hxx | ![]() |
00001 /************************************************************************/ 00002 /* */ 00003 /* Copyright 2010-2011 by Ullrich Koethe */ 00004 /* */ 00005 /* This file is part of the VIGRA computer vision library. */ 00006 /* The VIGRA Website is */ 00007 /* http://hci.iwr.uni-heidelberg.de/vigra/ */ 00008 /* Please direct questions, bug reports, and contributions to */ 00009 /* ullrich.koethe@iwr.uni-heidelberg.de or */ 00010 /* vigra@informatik.uni-hamburg.de */ 00011 /* */ 00012 /* Permission is hereby granted, free of charge, to any person */ 00013 /* obtaining a copy of this software and associated documentation */ 00014 /* files (the "Software"), to deal in the Software without */ 00015 /* restriction, including without limitation the rights to use, */ 00016 /* copy, modify, merge, publish, distribute, sublicense, and/or */ 00017 /* sell copies of the Software, and to permit persons to whom the */ 00018 /* Software is furnished to do so, subject to the following */ 00019 /* conditions: */ 00020 /* */ 00021 /* The above copyright notice and this permission notice shall be */ 00022 /* included in all copies or substantial portions of the */ 00023 /* Software. */ 00024 /* */ 00025 /* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND */ 00026 /* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES */ 00027 /* OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND */ 00028 /* NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT */ 00029 /* HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, */ 00030 /* WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING */ 00031 /* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR */ 00032 /* OTHER DEALINGS IN THE SOFTWARE. */ 00033 /* */ 00034 /************************************************************************/ 00035 00036 #ifndef VIGRA_BESSEL_HXX 00037 #define VIGRA_BESSEL_HXX 00038 00039 #include "mathutil.hxx" 00040 00041 #ifdef HasBoostMath 00042 #include <boost/math/special_functions/bessel.hpp> 00043 #endif 00044 00045 namespace vigra { 00046 00047 /** \addtogroup MathFunctions 00048 */ 00049 //@{ 00050 00051 namespace detail { 00052 00053 template <class REAL> 00054 int msta1(REAL x, int mp) 00055 { 00056 double a0,f0,f1,f; 00057 int i,n0,n1,nn; 00058 00059 a0 = abs(x); 00060 n0 = (int)(1.1*a0)+1; 00061 f0 = 0.5*std::log10(6.28*n0) - n0*std::log10(1.36*a0/n0)-mp; 00062 n1 = n0+5; 00063 f1 = 0.5*std::log10(6.28*n1) - n1*std::log10(1.36*a0/n1)-mp; 00064 for(i=0;i<20;i++) 00065 { 00066 nn = int(n1-(n1-n0)/(1.0-f0/f1)); 00067 f = 0.5*std::log10(6.28*nn) - nn*std::log10(1.36*a0/nn)-mp; 00068 if(abs(nn-n1) < 1) 00069 break; 00070 n0 = n1; 00071 f0 = f1; 00072 n1 = nn; 00073 f1 = f; 00074 } 00075 return nn; 00076 } 00077 00078 template <class REAL> 00079 int msta2(REAL x, int n, int mp) 00080 { 00081 double a0,ejn,hmp,f0,f1,f,obj; 00082 int i,n0,n1,nn; 00083 00084 a0 = abs(x); 00085 hmp = 0.5*mp; 00086 ejn = 0.5*std::log10(6.28*n) - n*std::log10(1.36*a0/n); 00087 if (ejn <= hmp) 00088 { 00089 obj = mp; 00090 n0 = (int)(1.1*a0); 00091 if (n0 < 1) 00092 n0 = 1; 00093 } 00094 else 00095 { 00096 obj = hmp+ejn; 00097 n0 = n; 00098 } 00099 f0 = 0.5*std::log10(6.28*n0) - n0*std::log10(1.36*a0/n0)-obj; 00100 n1 = n0+5; 00101 f1 = 0.5*std::log10(6.28*n1) - n1*std::log10(1.36*a0/n1)-obj; 00102 for (i=0;i<20;i++) 00103 { 00104 nn = int(n1-(n1-n0)/(1.0-f0/f1)); 00105 f = 0.5*std::log10(6.28*nn) - nn*std::log10(1.36*a0/nn)-obj; 00106 if (abs(nn-n1) < 1) 00107 break; 00108 n0 = n1; 00109 f0 = f1; 00110 n1 = nn; 00111 f1 = f; 00112 } 00113 return nn+10; 00114 } 00115 00116 // 00117 // INPUT: 00118 // double x -- argument of Bessel function of 1st and 2nd kind. 00119 // int n -- order 00120 // 00121 // OUPUT: 00122 // 00123 // int nm -- highest order actually computed (nm <= n) 00124 // double jn[] -- Bessel function of 1st kind, orders from 0 to nm 00125 // double yn[] -- Bessel function of 2nd kind, orders from 0 to nm 00126 // 00127 // Computes Bessel functions of all order up to 'n' using recurrence 00128 // relations. If 'nm' < 'n' only 'nm' orders are returned. 00129 // 00130 // code has been adapted from C.R. Bond's implementation 00131 // see http://www.crbond.com/math.htm 00132 // 00133 template <class REAL> 00134 void bessjyn(int n, REAL x,int &nm, double *jn, double *yn) 00135 { 00136 double t1,t2,f,f1,f2,bj0,bj1,bjk,by0,by1,cu,s0,su,sv; 00137 double ec,bs,byk,p0,p1,q0,q1; 00138 static double a[] = { 00139 -0.7031250000000000e-1, 00140 0.1121520996093750, 00141 -0.5725014209747314, 00142 6.074042001273483}; 00143 static double b[] = { 00144 0.7324218750000000e-1, 00145 -0.2271080017089844, 00146 1.727727502584457, 00147 -2.438052969955606e1}; 00148 static double a1[] = { 00149 0.1171875, 00150 -0.1441955566406250, 00151 0.6765925884246826, 00152 -6.883914268109947}; 00153 static double b1[] = { 00154 -0.1025390625, 00155 0.2775764465332031, 00156 -1.993531733751297, 00157 2.724882731126854e1}; 00158 00159 int i,k,m; 00160 nm = n; 00161 if (x < 1e-15) 00162 { 00163 for (i=0;i<=n;i++) 00164 { 00165 jn[i] = 0.0; 00166 yn[i] = -1e308; 00167 } 00168 jn[0] = 1.0; 00169 return; 00170 } 00171 if (x <= 300.0 || n > (int)(0.9*x)) 00172 { 00173 if (n == 0) 00174 nm = 1; 00175 m = msta1(x,200); 00176 if (m < nm) 00177 nm = m; 00178 else 00179 m = msta2(x,nm,15); 00180 bs = 0.0; 00181 su = 0.0; 00182 sv = 0.0; 00183 f2 = 0.0; 00184 f1 = 1.0e-100; 00185 for (k = m;k>=0;k--) 00186 { 00187 f = 2.0*(k+1.0)/x*f1 - f2; 00188 if (k <= nm) 00189 jn[k] = f; 00190 if ((k == 2*(int)(k/2)) && (k != 0)) 00191 { 00192 bs += 2.0*f; 00193 su += (-1)*((k & 2)-1)*f/(double)k; 00194 } 00195 else if (k > 1) 00196 { 00197 sv += (-1)*((k & 2)-1)*(double)k*f/(k*k-1.0); 00198 } 00199 f2 = f1; 00200 f1 = f; 00201 } 00202 s0 = bs+f; 00203 for (k=0;k<=nm;k++) 00204 { 00205 jn[k] /= s0; 00206 } 00207 ec = std::log(0.5*x) + M_EULER_GAMMA; 00208 by0 = M_2_PI*(ec*jn[0]-4.0*su/s0); 00209 yn[0] = by0; 00210 by1 = M_2_PI*((ec-1.0)*jn[1]-jn[0]/x-4.0*sv/s0); 00211 yn[1] = by1; 00212 } 00213 else 00214 { 00215 t1 = x-M_PI_4; 00216 p0 = 1.0; 00217 q0 = -0.125/x; 00218 for (k=0;k<4;k++) 00219 { 00220 p0 += a[k]*std::pow(x,-2*k-2); 00221 q0 += b[k]*std::pow(x,-2*k-3); 00222 } 00223 cu = std::sqrt(M_2_PI/x); 00224 bj0 = cu*(p0*std::cos(t1)-q0*std::sin(t1)); 00225 by0 = cu*(p0*std::sin(t1)+q0*std::cos(t1)); 00226 jn[0] = bj0; 00227 yn[0] = by0; 00228 t2 = x-0.75*M_PI; 00229 p1 = 1.0; 00230 q1 = 0.375/x; 00231 for (k=0;k<4;k++) 00232 { 00233 p1 += a1[k]*std::pow(x,-2*k-2); 00234 q1 += b1[k]*std::pow(x,-2*k-3); 00235 } 00236 bj1 = cu*(p1*std::cos(t2)-q1*std::sin(t2)); 00237 by1 = cu*(p1*std::sin(t2)+q1*std::cos(t2)); 00238 jn[1] = bj1; 00239 yn[1] = by1; 00240 for (k=2;k<=nm;k++) 00241 { 00242 bjk = 2.0*(k-1.0)*bj1/x-bj0; 00243 jn[k] = bjk; 00244 bj0 = bj1; 00245 bj1 = bjk; 00246 } 00247 } 00248 for (k=2;k<=nm;k++) 00249 { 00250 byk = 2.0*(k-1.0)*by1/x-by0; 00251 yn[k] = byk; 00252 by0 = by1; 00253 by1 = byk; 00254 } 00255 } 00256 00257 00258 00259 } // namespace detail 00260 00261 /*! Bessel function of the first kind. 00262 00263 Computes the value of BesselJ of integer order <tt>n</tt> and argument <tt>x</tt>. 00264 Negative <tt>x</tt> are unsupported and will result in a <tt>std::domain_error</tt>. 00265 00266 This function wraps a number of existing implementations and falls back to 00267 a rather slow algorithm if none of them is available. In particular, 00268 it uses boost::math when <tt>HasBoostMath</tt> is \#defined, or native 00269 implementations on gcc and MSVC otherwise. 00270 00271 <b>\#include</b> <vigra/bessel.hxx><br> 00272 Namespace: vigra 00273 */ 00274 inline double besselJ(int n, double x) 00275 { 00276 if(x < 0.0) 00277 throw std::domain_error("besselJ(n, x): x cannot be negative"); 00278 if(x < 1e-15) 00279 return n == 0 ? 1.0 : 0.0; 00280 #if defined(HasBoostMath) 00281 return boost::math::cyl_bessel_j((double)n, x); 00282 #elif defined(__GNUC__) 00283 return ::jn(n, x); 00284 #elif defined(_MSC_VER) 00285 return _jn(n, x); 00286 #else 00287 int an = abs(n), nr = n, s = an+2; 00288 ArrayVector<double> t(2*s); 00289 detail::bessjyn(an, x, nr, &t[0], &t[s]); 00290 if(n < 0 && odd(an)) 00291 return -t[an]; 00292 else 00293 return t[an]; 00294 #endif 00295 } 00296 00297 /*! Bessel function of the second kind. 00298 00299 Computes the value of BesselY of integer order <tt>n</tt> and argument <tt>x</tt>. 00300 Negative <tt>x</tt> are unsupported and will result in a <tt>std::domain_error</tt>. 00301 00302 This function wraps a number of existing implementations and falls back to 00303 a rather slow algorithm if none of them is available. In particular, 00304 it uses boost::math when <tt>HasBoostMath</tt> is \#defined, or native 00305 implementations on gcc and MSVC otherwise. 00306 00307 <b>\#include</b> <vigra/bessel.hxx><br> 00308 Namespace: vigra 00309 */ 00310 inline double besselY(int n, double x) 00311 { 00312 if(x < 0.0) 00313 throw std::domain_error("besselY(n, x): x cannot be negative"); 00314 if(x == 0.0 ) 00315 return -std::numeric_limits<double>::infinity(); 00316 #if defined(HasBoostMath) 00317 return boost::math::cyl_neumann((double)n, x); 00318 #elif defined(__GNUC__) 00319 return ::yn(n, x); 00320 #elif defined(_MSC_VER) 00321 return _yn(n, x); 00322 #else 00323 int an = abs(n), nr = n, s = an+2; 00324 ArrayVector<double> t(2*s); 00325 detail::bessjyn(an, x, nr, &t[0], &t[s]); 00326 if(n < 0.0 && odd(n)) 00327 return -t[an+s]; 00328 else 00329 return t[an+s]; 00330 #endif 00331 } 00332 00333 //@} 00334 00335 } // namespace vigra 00336 00337 #endif // VIGRA_BESSEL_HXX
© Ullrich Köthe (ullrich.koethe@iwr.uni-heidelberg.de) |
html generated using doxygen and Python
|