scaddins/source/analysis/bessel.cxx

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