1 /************************************************************** 2 * 3 * Licensed to the Apache Software Foundation (ASF) under one 4 * or more contributor license agreements. See the NOTICE file 5 * distributed with this work for additional information 6 * regarding copyright ownership. The ASF licenses this file 7 * to you under the Apache License, Version 2.0 (the 8 * "License"); you may not use this file except in compliance 9 * with the License. You may obtain a copy of the License at 10 * 11 * http://www.apache.org/licenses/LICENSE-2.0 12 * 13 * Unless required by applicable law or agreed to in writing, 14 * software distributed under the License is distributed on an 15 * "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY 16 * KIND, either express or implied. See the License for the 17 * specific language governing permissions and limitations 18 * under the License. 19 * 20 *************************************************************/ 21 22 23 24 #include "bessel.hxx" 25 #include "analysishelper.hxx" 26 27 #include <rtl/math.hxx> 28 29 using ::com::sun::star::lang::IllegalArgumentException; 30 using ::com::sun::star::sheet::NoConvergenceException; 31 32 namespace sca { 33 namespace analysis { 34 35 // ============================================================================ 36 37 const double f_PI = 3.1415926535897932385; 38 const double f_2_PI = 2.0 * f_PI; 39 const double f_PI_DIV_2 = f_PI / 2.0; 40 const double f_PI_DIV_4 = f_PI / 4.0; 41 const double f_2_DIV_PI = 2.0 / f_PI; 42 43 const double THRESHOLD = 30.0; // Threshold for usage of approximation formula. 44 const double MAXEPSILON = 1e-10; // Maximum epsilon for end of iteration. 45 const sal_Int32 MAXITER = 100; // Maximum number of iterations. 46 47 // ============================================================================ 48 // BESSEL J 49 // ============================================================================ 50 51 /* The BESSEL function, first kind, unmodified: 52 The algorithm follows 53 http://www.reference-global.com/isbn/978-3-11-020354-7 54 Numerical Mathematics 1 / Numerische Mathematik 1, 55 An algorithm-based introduction / Eine algorithmisch orientierte Einfuehrung 56 Deuflhard, Peter; Hohmann, Andreas 57 Berlin, New York (Walter de Gruyter) 2008 58 4. ueberarb. u. erw. Aufl. 2008 59 eBook ISBN: 978-3-11-020355-4 60 Chapter 6.3.2 , algorithm 6.24 61 The source is in German. 62 The BesselJ-function is a special case of the adjoint summation with 63 a_k = 2*(k-1)/x for k=1,... 64 b_k = -1, for all k, directly substituted 65 m_0=1, m_k=2 for k even, and m_k=0 for k odd, calculated on the fly 66 alpha_k=1 for k=N and alpha_k=0 otherwise 67 */ 68 69 // ---------------------------------------------------------------------------- 70 71 double BesselJ( double x, sal_Int32 N ) throw (IllegalArgumentException, NoConvergenceException) 72 73 { 74 if( N < 0 ) 75 throw IllegalArgumentException(); 76 if (x==0.0) 77 return (N==0) ? 1.0 : 0.0; 78 79 /* The algorithm works only for x>0, therefore remember sign. BesselJ 80 with integer order N is an even function for even N (means J(-x)=J(x)) 81 and an odd function for odd N (means J(-x)=-J(x)).*/ 82 double fSign = (N % 2 == 1 && x < 0) ? -1.0 : 1.0; 83 double fX = fabs(x); 84 85 const double fMaxIteration = 9000000.0; //experimental, for to return in < 3 seconds 86 double fEstimateIteration = fX * 1.5 + N; 87 bool bAsymptoticPossible = pow(fX,0.4) > N; 88 if (fEstimateIteration > fMaxIteration) 89 { 90 if (bAsymptoticPossible) 91 return fSign * sqrt(f_2_DIV_PI/fX)* cos(fX-N*f_PI_DIV_2-f_PI_DIV_4); 92 else 93 throw NoConvergenceException(); 94 } 95 96 double epsilon = 1.0e-15; // relative error 97 bool bHasfound = false; 98 double k= 0.0; 99 // e_{-1} = 0; e_0 = alpha_0 / b_2 100 double u ; // u_0 = e_0/f_0 = alpha_0/m_0 = alpha_0 101 102 // first used with k=1 103 double m_bar; // m_bar_k = m_k * f_bar_{k-1} 104 double g_bar; // g_bar_k = m_bar_k - a_{k+1} + g_{k-1} 105 double g_bar_delta_u; // g_bar_delta_u_k = f_bar_{k-1} * alpha_k 106 // - g_{k-1} * delta_u_{k-1} - m_bar_k * u_{k-1} 107 // f_{-1} = 0.0; f_0 = m_0 / b_2 = 1/(-1) = -1 108 double g = 0.0; // g_0= f_{-1} / f_0 = 0/(-1) = 0 109 double delta_u = 0.0; // dummy initialize, first used with * 0 110 double f_bar = -1.0; // f_bar_k = 1/f_k, but only used for k=0 111 112 if (N==0) 113 { 114 //k=0; alpha_0 = 1.0 115 u = 1.0; // u_0 = alpha_0 116 // k = 1.0; at least one step is necessary 117 // m_bar_k = m_k * f_bar_{k-1} ==> m_bar_1 = 0.0 118 g_bar_delta_u = 0.0; // alpha_k = 0.0, m_bar = 0.0; g= 0.0 119 g_bar = - 2.0/fX; // k = 1.0, g = 0.0 120 delta_u = g_bar_delta_u / g_bar; 121 u = u + delta_u ; // u_k = u_{k-1} + delta_u_k 122 g = -1.0 / g_bar; // g_k=b_{k+2}/g_bar_k 123 f_bar = f_bar * g; // f_bar_k = f_bar_{k-1}* g_k 124 k = 2.0; 125 // From now on all alpha_k = 0.0 and k > N+1 126 } 127 else 128 { // N >= 1 and alpha_k = 0.0 for k<N 129 u=0.0; // u_0 = alpha_0 130 for (k =1.0; k<= N-1; k = k + 1.0) 131 { 132 m_bar=2.0 * fmod(k-1.0, 2.0) * f_bar; 133 g_bar_delta_u = - g * delta_u - m_bar * u; // alpha_k = 0.0 134 g_bar = m_bar - 2.0*k/fX + g; 135 delta_u = g_bar_delta_u / g_bar; 136 u = u + delta_u; 137 g = -1.0/g_bar; 138 f_bar=f_bar * g; 139 } 140 // Step alpha_N = 1.0 141 m_bar=2.0 * fmod(k-1.0, 2.0) * f_bar; 142 g_bar_delta_u = f_bar - g * delta_u - m_bar * u; // alpha_k = 1.0 143 g_bar = m_bar - 2.0*k/fX + g; 144 delta_u = g_bar_delta_u / g_bar; 145 u = u + delta_u; 146 g = -1.0/g_bar; 147 f_bar = f_bar * g; 148 k = k + 1.0; 149 } 150 // Loop until desired accuracy, always alpha_k = 0.0 151 do 152 { 153 m_bar = 2.0 * fmod(k-1.0, 2.0) * f_bar; 154 g_bar_delta_u = - g * delta_u - m_bar * u; 155 g_bar = m_bar - 2.0*k/fX + g; 156 delta_u = g_bar_delta_u / g_bar; 157 u = u + delta_u; 158 g = -1.0/g_bar; 159 f_bar = f_bar * g; 160 bHasfound = (fabs(delta_u)<=fabs(u)*epsilon); 161 k = k + 1.0; 162 } 163 while (!bHasfound && k <= fMaxIteration); 164 if (bHasfound) 165 return u * fSign; 166 else 167 throw NoConvergenceException(); // unlikely to happen 168 } 169 170 // ============================================================================ 171 // BESSEL I 172 // ============================================================================ 173 174 /* The BESSEL function, first kind, modified: 175 176 inf (x/2)^(n+2k) 177 I_n(x) = SUM TERM(n,k) with TERM(n,k) := -------------- 178 k=0 k! (n+k)! 179 180 No asymptotic approximation used, see issue 43040. 181 */ 182 183 // ---------------------------------------------------------------------------- 184 185 double BesselI( double x, sal_Int32 n ) throw( IllegalArgumentException, NoConvergenceException ) 186 { 187 const double fEpsilon = 1.0E-15; 188 const sal_Int32 nMaxIteration = 2000; 189 const double fXHalf = x / 2.0; 190 if( n < 0 ) 191 throw IllegalArgumentException(); 192 193 double fResult = 0.0; 194 195 /* Start the iteration without TERM(n,0), which is set here. 196 197 TERM(n,0) = (x/2)^n / n! 198 */ 199 sal_Int32 nK = 0; 200 double fTerm = 1.0; 201 // avoid overflow in Fak(n) 202 for( nK = 1; nK <= n; ++nK ) 203 { 204 fTerm = fTerm / static_cast< double >( nK ) * fXHalf; 205 } 206 fResult = fTerm; // Start result with TERM(n,0). 207 if( fTerm != 0.0 ) 208 { 209 nK = 1; 210 do 211 { 212 /* Calculation of TERM(n,k) from TERM(n,k-1): 213 214 (x/2)^(n+2k) 215 TERM(n,k) = -------------- 216 k! (n+k)! 217 218 (x/2)^2 (x/2)^(n+2(k-1)) 219 = -------------------------- 220 k (k-1)! (n+k) (n+k-1)! 221 222 (x/2)^2 (x/2)^(n+2(k-1)) 223 = --------- * ------------------ 224 k(n+k) (k-1)! (n+k-1)! 225 226 x^2/4 227 = -------- TERM(n,k-1) 228 k(n+k) 229 */ 230 fTerm = fTerm * fXHalf / static_cast<double>(nK) * fXHalf / static_cast<double>(nK+n); 231 fResult += fTerm; 232 nK++; 233 } 234 while( (fabs( fTerm ) > fabs(fResult) * fEpsilon) && (nK < nMaxIteration) ); 235 236 } 237 return fResult; 238 } 239 240 241 // ============================================================================ 242 243 double Besselk0( double fNum ) throw( IllegalArgumentException, NoConvergenceException ) 244 { 245 double fRet; 246 247 if( fNum <= 2.0 ) 248 { 249 double fNum2 = fNum * 0.5; 250 double y = fNum2 * fNum2; 251 252 fRet = -log( fNum2 ) * BesselI( fNum, 0 ) + 253 ( -0.57721566 + y * ( 0.42278420 + y * ( 0.23069756 + y * ( 0.3488590e-1 + 254 y * ( 0.262698e-2 + y * ( 0.10750e-3 + y * 0.74e-5 ) ) ) ) ) ); 255 } 256 else 257 { 258 double y = 2.0 / fNum; 259 260 fRet = exp( -fNum ) / sqrt( fNum ) * ( 1.25331414 + y * ( -0.7832358e-1 + 261 y * ( 0.2189568e-1 + y * ( -0.1062446e-1 + y * ( 0.587872e-2 + 262 y * ( -0.251540e-2 + y * 0.53208e-3 ) ) ) ) ) ); 263 } 264 265 return fRet; 266 } 267 268 269 double Besselk1( double fNum ) throw( IllegalArgumentException, NoConvergenceException ) 270 { 271 double fRet; 272 273 if( fNum <= 2.0 ) 274 { 275 double fNum2 = fNum * 0.5; 276 double y = fNum2 * fNum2; 277 278 fRet = log( fNum2 ) * BesselI( fNum, 1 ) + 279 ( 1.0 + y * ( 0.15443144 + y * ( -0.67278579 + y * ( -0.18156897 + y * ( -0.1919402e-1 + 280 y * ( -0.110404e-2 + y * ( -0.4686e-4 ) ) ) ) ) ) ) 281 / fNum; 282 } 283 else 284 { 285 double y = 2.0 / fNum; 286 287 fRet = exp( -fNum ) / sqrt( fNum ) * ( 1.25331414 + y * ( 0.23498619 + 288 y * ( -0.3655620e-1 + y * ( 0.1504268e-1 + y * ( -0.780353e-2 + 289 y * ( 0.325614e-2 + y * ( -0.68245e-3 ) ) ) ) ) ) ); 290 } 291 292 return fRet; 293 } 294 295 296 double BesselK( double fNum, sal_Int32 nOrder ) throw( IllegalArgumentException, NoConvergenceException ) 297 { 298 switch( nOrder ) 299 { 300 case 0: return Besselk0( fNum ); 301 case 1: return Besselk1( fNum ); 302 default: 303 { 304 double fBkp; 305 306 double fTox = 2.0 / fNum; 307 double fBkm = Besselk0( fNum ); 308 double fBk = Besselk1( fNum ); 309 310 for( sal_Int32 n = 1 ; n < nOrder ; n++ ) 311 { 312 fBkp = fBkm + double( n ) * fTox * fBk; 313 fBkm = fBk; 314 fBk = fBkp; 315 } 316 317 return fBk; 318 } 319 } 320 } 321 322 // ============================================================================ 323 // BESSEL Y 324 // ============================================================================ 325 326 /* The BESSEL function, second kind, unmodified: 327 The algorithm for order 0 and for order 1 follows 328 http://www.reference-global.com/isbn/978-3-11-020354-7 329 Numerical Mathematics 1 / Numerische Mathematik 1, 330 An algorithm-based introduction / Eine algorithmisch orientierte Einfuehrung 331 Deuflhard, Peter; Hohmann, Andreas 332 Berlin, New York (Walter de Gruyter) 2008 333 4. ueberarb. u. erw. Aufl. 2008 334 eBook ISBN: 978-3-11-020355-4 335 Chapter 6.3.2 , algorithm 6.24 336 The source is in German. 337 See #i31656# for a commented version of the implementation, attachment #desc6 338 http://www.openoffice.org/nonav/issues/showattachment.cgi/63609/Comments%20to%20the%20implementation%20of%20the%20Bessel%20functions.odt 339 */ 340 341 double Bessely0( double fX ) throw( IllegalArgumentException, NoConvergenceException ) 342 { 343 if (fX <= 0) 344 throw IllegalArgumentException(); 345 const double fMaxIteration = 9000000.0; // should not be reached 346 if (fX > 5.0e+6) // iteration is not considerable better then approximation 347 return sqrt(1/f_PI/fX) 348 *(rtl::math::sin(fX)-rtl::math::cos(fX)); 349 const double epsilon = 1.0e-15; 350 const double EulerGamma = 0.57721566490153286060; 351 double alpha = log(fX/2.0)+EulerGamma; 352 double u = alpha; 353 354 double k = 1.0; 355 double m_bar = 0.0; 356 double g_bar_delta_u = 0.0; 357 double g_bar = -2.0 / fX; 358 double delta_u = g_bar_delta_u / g_bar; 359 double g = -1.0/g_bar; 360 double f_bar = -1 * g; 361 362 double sign_alpha = 1.0; 363 double km1mod2; 364 bool bHasFound = false; 365 k = k + 1; 366 do 367 { 368 km1mod2 = fmod(k-1.0,2.0); 369 m_bar=(2.0*km1mod2) * f_bar; 370 if (km1mod2 == 0.0) 371 alpha = 0.0; 372 else 373 { 374 alpha = sign_alpha * (4.0/k); 375 sign_alpha = -sign_alpha; 376 } 377 g_bar_delta_u = f_bar * alpha - g * delta_u - m_bar * u; 378 g_bar = m_bar - (2.0*k)/fX + g; 379 delta_u = g_bar_delta_u / g_bar; 380 u = u+delta_u; 381 g = -1.0 / g_bar; 382 f_bar = f_bar*g; 383 bHasFound = (fabs(delta_u)<=fabs(u)*epsilon); 384 k=k+1; 385 } 386 while (!bHasFound && k<fMaxIteration); 387 if (bHasFound) 388 return u*f_2_DIV_PI; 389 else 390 throw NoConvergenceException(); // not likely to happen 391 } 392 393 // See #i31656# for a commented version of this implementation, attachment #desc6 394 // http://www.openoffice.org/nonav/issues/showattachment.cgi/63609/Comments%20to%20the%20implementation%20of%20the%20Bessel%20functions.odt 395 double Bessely1( double fX ) throw( IllegalArgumentException, NoConvergenceException ) 396 { 397 if (fX <= 0) 398 throw IllegalArgumentException(); 399 const double fMaxIteration = 9000000.0; // should not be reached 400 if (fX > 5.0e+6) // iteration is not considerable better then approximation 401 return - sqrt(1/f_PI/fX) 402 *(rtl::math::sin(fX)+rtl::math::cos(fX)); 403 const double epsilon = 1.0e-15; 404 const double EulerGamma = 0.57721566490153286060; 405 double alpha = 1.0/fX; 406 double f_bar = -1.0; 407 double g = 0.0; 408 double u = alpha; 409 double k = 1.0; 410 double m_bar = 0.0; 411 alpha = 1.0 - EulerGamma - log(fX/2.0); 412 double g_bar_delta_u = -alpha; 413 double g_bar = -2.0 / fX; 414 double delta_u = g_bar_delta_u / g_bar; 415 u = u + delta_u; 416 g = -1.0/g_bar; 417 f_bar = f_bar * g; 418 double sign_alpha = -1.0; 419 double km1mod2; //will be (k-1) mod 2 420 double q; // will be (k-1) div 2 421 bool bHasFound = false; 422 k = k + 1.0; 423 do 424 { 425 km1mod2 = fmod(k-1.0,2.0); 426 m_bar=(2.0*km1mod2) * f_bar; 427 q = (k-1.0)/2.0; 428 if (km1mod2 == 0.0) // k is odd 429 { 430 alpha = sign_alpha * (1.0/q + 1.0/(q+1.0)); 431 sign_alpha = -sign_alpha; 432 } 433 else 434 alpha = 0.0; 435 g_bar_delta_u = f_bar * alpha - g * delta_u - m_bar * u; 436 g_bar = m_bar - (2.0*k)/fX + g; 437 delta_u = g_bar_delta_u / g_bar; 438 u = u+delta_u; 439 g = -1.0 / g_bar; 440 f_bar = f_bar*g; 441 bHasFound = (fabs(delta_u)<=fabs(u)*epsilon); 442 k=k+1; 443 } 444 while (!bHasFound && k<fMaxIteration); 445 if (bHasFound) 446 return -u*2.0/f_PI; 447 else 448 throw NoConvergenceException(); 449 } 450 451 double BesselY( double fNum, sal_Int32 nOrder ) throw( IllegalArgumentException, NoConvergenceException ) 452 { 453 switch( nOrder ) 454 { 455 case 0: return Bessely0( fNum ); 456 case 1: return Bessely1( fNum ); 457 default: 458 { 459 double fByp; 460 461 double fTox = 2.0 / fNum; 462 double fBym = Bessely0( fNum ); 463 double fBy = Bessely1( fNum ); 464 465 for( sal_Int32 n = 1 ; n < nOrder ; n++ ) 466 { 467 fByp = double( n ) * fTox * fBy - fBym; 468 fBym = fBy; 469 fBy = fByp; 470 } 471 472 return fBy; 473 } 474 } 475 } 476 477 // ============================================================================ 478 479 } // namespace analysis 480 } // namespace sca 481 482