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