scaddins/source/analysis/bessel.cxx

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