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