scaddins/source/analysis/bessel.cxx

   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 Einführung
  60     Deuflhard, Peter; Hohmann, Andreas
  61     Berlin, New York (Walter de Gruyter) 2008
  62     4. überarb. 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     Approximation for the BESSEL function, first kind, modified, for great x:
 185
 186         I_n(x)  ~  e^x / sqrt( 2 PI x )  for  x>=0.
 187  */
 188
 189 // ----------------------------------------------------------------------------
 190
 191 double BesselI( double x, sal_Int32 n ) throw( IllegalArgumentException, NoConvergenceException )
 192 {
 193     if( n < 0 )
 194         throw IllegalArgumentException();
 195
 196     double fResult = 0.0;
 197     if( fabs( x ) <= THRESHOLD )
 198     {
 199         /*  Start the iteration without TERM(n,0), which is set here.
 200
 201                 TERM(n,0) = (x/2)^n / n!
 202          */
 203         double fTerm = pow( x / 2.0, (double)n ) / Fak( n );
 204         sal_Int32 nK = 1;   // Start the iteration with k=1.
 205         fResult = fTerm;    // Start result with TERM(n,0).
 206
 207         const double fSqrX = x * x / 4.0;
 208
 209         do
 210         {
 211             /*  Calculation of TERM(n,k) from TERM(n,k-1):
 212
 213                                    (x/2)^(n+2k)
 214                     TERM(n,k)  =  --------------
 215                                     k! (n+k)!
 216
 217                                    (x/2)^2 (x/2)^(n+2(k-1))
 218                                =  --------------------------
 219                                    k (k-1)! (n+k) (n+k-1)!
 220
 221                                    (x/2)^2     (x/2)^(n+2(k-1))
 222                                =  --------- * ------------------
 223                                    k(n+k)      (k-1)! (n+k-1)!
 224
 225                                    x^2/4
 226                                =  -------- TERM(n,k-1)
 227                                    k(n+k)
 228              */
 229             fTerm *= fSqrX;    // defined above as x^2/4
 230             fTerm /= (nK * (nK + n));
 231             fResult += fTerm;
 232         }
 233         while( (fabs( fTerm ) > MAXEPSILON) && (++nK < MAXITER) );
 234     }
 235     else
 236     {
 237         /*  Approximation for the BESSEL function, first kind, modified:
 238
 239                 I_n(x)  ~  e^x / sqrt( 2 PI x )  for  x>=0.
 240
 241             The BESSEL function I_n with n IN {0,2,4,...} is axially symmetric at
 242             x=0, means I_n(x) = I_n(-x). Therefore the approximation for x<0 is:
 243
 244                 I_n(x)  =  I_n(|x|)  for  x<0  and  n IN {0,2,4,...}.
 245
 246             The BESSEL function I_n with n IN {1,3,5,...} is point-symmetric at
 247             x=0, means I_n(x) = -I_n(-x). Therefore the approximation for x<0 is:
 248
 249                 I_n(x)  =  -I_n(|x|)  for  x<0  and  n IN {1,3,5,...}.
 250          */
 251         double fXAbs = fabs( x );
 252         fResult = exp( fXAbs ) / sqrt( f_2_PI * fXAbs );
 253         if( (n & 1) && (x < 0.0) )
 254             fResult = -fResult;
 255     }
 256     return fResult;
 257 }
 258
 259
 260 // ============================================================================
 261
 262 double Besselk0( double fNum ) throw( IllegalArgumentException, NoConvergenceException )
 263 {
 264     double  fRet;
 265
 266     if( fNum <= 2.0 )
 267     {
 268         double  fNum2 = fNum * 0.5;
 269         double  y = fNum2 * fNum2;
 270
 271         fRet = -log( fNum2 ) * BesselI( fNum, 0 ) +
 272                 ( -0.57721566 + y * ( 0.42278420 + y * ( 0.23069756 + y * ( 0.3488590e-1 +
 273                     y * ( 0.262698e-2 + y * ( 0.10750e-3 + y * 0.74e-5 ) ) ) ) ) );
 274     }
 275     else
 276     {
 277         double  y = 2.0 / fNum;
 278
 279         fRet = exp( -fNum ) / sqrt( fNum ) * ( 1.25331414 + y * ( -0.7832358e-1 +
 280                 y * ( 0.2189568e-1 + y * ( -0.1062446e-1 + y * ( 0.587872e-2 +
 281                 y * ( -0.251540e-2 + y * 0.53208e-3 ) ) ) ) ) );
 282     }
 283
 284     return fRet;
 285 }
 286
 287
 288 double Besselk1( double fNum ) throw( IllegalArgumentException, NoConvergenceException )
 289 {
 290     double  fRet;
 291
 292     if( fNum <= 2.0 )
 293     {
 294         double  fNum2 = fNum * 0.5;
 295         double  y = fNum2 * fNum2;
 296
 297         fRet = log( fNum2 ) * BesselI( fNum, 1 ) +
 298                 ( 1.0 + y * ( 0.15443144 + y * ( -0.67278579 + y * ( -0.18156897 + y * ( -0.1919402e-1 +
 299                     y * ( -0.110404e-2 + y * ( -0.4686e-4 ) ) ) ) ) ) )
 300                 / fNum;
 301     }
 302     else
 303     {
 304         double  y = 2.0 / fNum;
 305
 306         fRet = exp( -fNum ) / sqrt( fNum ) * ( 1.25331414 + y * ( 0.23498619 +
 307                 y * ( -0.3655620e-1 + y * ( 0.1504268e-1 + y * ( -0.780353e-2 +
 308                 y * ( 0.325614e-2 + y * ( -0.68245e-3 ) ) ) ) ) ) );
 309     }
 310
 311     return fRet;
 312 }
 313
 314
 315 double BesselK( double fNum, sal_Int32 nOrder ) throw( IllegalArgumentException, NoConvergenceException )
 316 {
 317     switch( nOrder )
 318     {
 319         case 0:     return Besselk0( fNum );
 320         case 1:     return Besselk1( fNum );
 321         default:
 322         {
 323             double      fBkp;
 324
 325             double      fTox = 2.0 / fNum;
 326             double      fBkm = Besselk0( fNum );
 327             double      fBk = Besselk1( fNum );
 328
 329             for( sal_Int32 n = 1 ; n < nOrder ; n++ )
 330             {
 331                 fBkp = fBkm + double( n ) * fTox * fBk;
 332                 fBkm = fBk;
 333                 fBk = fBkp;
 334             }
 335
 336             return fBk;
 337         }
 338     }
 339 }
 340
 341 // ============================================================================
 342 // BESSEL Y
 343 // ============================================================================
 344
 345 /*  The BESSEL function, second kind, unmodified:
 346     The algorithm for order 0 and for order 1 follows
 347     http://www.reference-global.com/isbn/978-3-11-020354-7
 348     Numerical Mathematics 1 / Numerische Mathematik 1,
 349     An algorithm-based introduction / Eine algorithmisch orientierte Einführung
 350     Deuflhard, Peter; Hohmann, Andreas
 351     Berlin, New York (Walter de Gruyter) 2008
 352     4. überarb. u. erw. Aufl. 2008
 353     eBook ISBN: 978-3-11-020355-4
 354     Chapter 6.3.2 , algorithm 6.24
 355     The source is in German.
 356     See #i31656# for a commented version of the implementation, attachment #desc6
 357     http://www.openoffice.org/nonav/issues/showattachment.cgi/63609/Comments%20to%20the%20implementation%20of%20the%20Bessel%20functions.odt
 358 */
 359
 360 double Bessely0( double fX ) throw( IllegalArgumentException, NoConvergenceException )
 361 {
 362     if (fX <= 0)
 363         throw IllegalArgumentException();
 364     const double fMaxIteration = 9000000.0; // should not be reached
 365     if (fX > 5.0e+6) // iteration is not considerable better then approximation
 366         return sqrt(1/f_PI/fX)
 367                 *(rtl::math::sin(fX)-rtl::math::cos(fX));
 368     const double epsilon = 1.0e-15;
 369     const double EulerGamma = 0.57721566490153286060;
 370     double alpha = log(fX/2.0)+EulerGamma;
 371     double u = alpha;
 372
 373     double k = 1.0;
 374     double m_bar = 0.0;
 375     double g_bar_delta_u = 0.0;
 376     double g_bar = -2.0 / fX;
 377     double delta_u = g_bar_delta_u / g_bar;
 378     double g = -1.0/g_bar;
 379     double f_bar = -1 * g;
 380
 381     double sign_alpha = 1.0;
 382     double km1mod2;
 383     bool bHasFound = false;
 384     k = k + 1;
 385     do
 386     {
 387         km1mod2 = fmod(k-1.0,2.0);
 388         m_bar=(2.0*km1mod2) * f_bar;
 389         if (km1mod2 == 0.0)
 390             alpha = 0.0;
 391         else
 392         {
 393            alpha = sign_alpha * (4.0/k);
 394            sign_alpha = -sign_alpha;
 395         }
 396         g_bar_delta_u = f_bar * alpha - g * delta_u - m_bar * u;
 397         g_bar = m_bar - (2.0*k)/fX + g;
 398         delta_u = g_bar_delta_u / g_bar;
 399         u = u+delta_u;
 400         g = -1.0 / g_bar;
 401         f_bar = f_bar*g;
 402         bHasFound = (fabs(delta_u)<=fabs(u)*epsilon);
 403         k=k+1;
 404     }
 405     while (!bHasFound && k<fMaxIteration);
 406     if (bHasFound)
 407         return u*f_2_DIV_PI;
 408     else
 409         throw NoConvergenceException(); // not likely to happen
 410 }
 411
 412 // See #i31656# for a commented version of this implementation, attachment #desc6
 413 // http://www.openoffice.org/nonav/issues/showattachment.cgi/63609/Comments%20to%20the%20implementation%20of%20the%20Bessel%20functions.odt
 414 double Bessely1( double fX ) throw( IllegalArgumentException, NoConvergenceException )
 415 {
 416     if (fX <= 0)
 417         throw IllegalArgumentException();
 418     const double fMaxIteration = 9000000.0; // should not be reached
 419     if (fX > 5.0e+6) // iteration is not considerable better then approximation
 420         return - sqrt(1/f_PI/fX)
 421                 *(rtl::math::sin(fX)+rtl::math::cos(fX));
 422     const double epsilon = 1.0e-15;
 423     const double EulerGamma = 0.57721566490153286060;
 424     double alpha = 1.0/fX;
 425     double f_bar = -1.0;
 426     double g = 0.0;
 427     double u = alpha;
 428     double k = 1.0;
 429     double m_bar = 0.0;
 430     alpha = 1.0 - EulerGamma - log(fX/2.0);
 431     double g_bar_delta_u = -alpha;
 432     double g_bar = -2.0 / fX;
 433     double delta_u = g_bar_delta_u / g_bar;
 434     u = u + delta_u;
 435     g = -1.0/g_bar;
 436     f_bar = f_bar * g;
 437     double sign_alpha = -1.0;
 438     double km1mod2; //will be (k-1) mod 2
 439     double q; // will be (k-1) div 2
 440     bool bHasFound = false;
 441     k = k + 1.0;
 442     do
 443     {
 444         km1mod2 = fmod(k-1.0,2.0);
 445         m_bar=(2.0*km1mod2) * f_bar;
 446         q = (k-1.0)/2.0;
 447         if (km1mod2 == 0.0) // k is odd
 448         {
 449            alpha = sign_alpha * (1.0/q + 1.0/(q+1.0));
 450            sign_alpha = -sign_alpha;
 451         }
 452         else
 453             alpha = 0.0;
 454         g_bar_delta_u = f_bar * alpha - g * delta_u - m_bar * u;
 455         g_bar = m_bar - (2.0*k)/fX + g;
 456         delta_u = g_bar_delta_u / g_bar;
 457         u = u+delta_u;
 458         g = -1.0 / g_bar;
 459         f_bar = f_bar*g;
 460         bHasFound = (fabs(delta_u)<=fabs(u)*epsilon);
 461         k=k+1;
 462     }
 463     while (!bHasFound && k<fMaxIteration);
 464     if (bHasFound)
 465         return -u*2.0/f_PI;
 466     else
 467         throw NoConvergenceException();
 468 }
 469
 470 double BesselY( double fNum, sal_Int32 nOrder ) throw( IllegalArgumentException, NoConvergenceException )
 471 {
 472     switch( nOrder )
 473     {
 474         case 0:     return Bessely0( fNum );
 475         case 1:     return Bessely1( fNum );
 476         default:
 477         {
 478             double      fByp;
 479
 480             double      fTox = 2.0 / fNum;
 481             double      fBym = Bessely0( fNum );
 482             double      fBy = Bessely1( fNum );
 483
 484             for( sal_Int32 n = 1 ; n < nOrder ; n++ )
 485             {
 486                 fByp = double( n ) * fTox * fBy - fBym;
 487                 fBym = fBy;
 488                 fBy = fByp;
 489             }
 490
 491             return fBy;
 492         }
 493     }
 494 }
 495
 496 // ============================================================================
 497
 498 } // namespace analysis
 499 } // namespace sca
 500