1 /* -*- Mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
3 * This file is part of the LibreOffice project.
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/.
9 * This file incorporates work covered by the following license notice:
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 .
22 #include <rtl/math.hxx>
24 #include <com/sun/star/lang/IllegalArgumentException.hpp>
25 #include <com/sun/star/sheet/NoConvergenceException.hpp>
27 using ::com::sun::star::lang::IllegalArgumentException
;
28 using ::com::sun::star::sheet::NoConvergenceException
;
33 const double f_PI
= 3.1415926535897932385;
34 const double f_PI_DIV_2
= f_PI
/ 2.0;
35 const double f_PI_DIV_4
= f_PI
/ 4.0;
36 const double f_2_DIV_PI
= 2.0 / f_PI
;
42 /* The BESSEL function, first kind, unmodified:
44 http://www.reference-global.com/isbn/978-3-11-020354-7
45 Numerical Mathematics 1 / Numerische Mathematik 1,
46 An algorithm-based introduction / Eine algorithmisch orientierte Einfuehrung
47 Deuflhard, Peter; Hohmann, Andreas
48 Berlin, New York (Walter de Gruyter) 2008
49 4. ueberarb. u. erw. Aufl. 2008
50 eBook ISBN: 978-3-11-020355-4
51 Chapter 6.3.2 , algorithm 6.24
52 The source is in German.
53 The BesselJ-function is a special case of the adjoint summation with
54 a_k = 2*(k-1)/x for k=1,...
55 b_k = -1, for all k, directly substituted
56 m_0=1, m_k=2 for k even, and m_k=0 for k odd, calculated on the fly
57 alpha_k=1 for k=N and alpha_k=0 otherwise
60 double BesselJ( double x
, sal_Int32 N
)
64 throw IllegalArgumentException();
66 return (N
==0) ? 1.0 : 0.0;
68 /* The algorithm works only for x>0, therefore remember sign. BesselJ
69 with integer order N is an even function for even N (means J(-x)=J(x))
70 and an odd function for odd N (means J(-x)=-J(x)).*/
71 double fSign
= (N
% 2 == 1 && x
< 0) ? -1.0 : 1.0;
74 const double fMaxIteration
= 9000000.0; //experimental, for to return in < 3 seconds
75 double fEstimateIteration
= fX
* 1.5 + N
;
76 bool bAsymptoticPossible
= pow(fX
,0.4) > N
;
77 if (fEstimateIteration
> fMaxIteration
)
79 if (!bAsymptoticPossible
)
80 throw NoConvergenceException();
81 return fSign
* sqrt(f_2_DIV_PI
/fX
)* cos(fX
-N
*f_PI_DIV_2
-f_PI_DIV_4
);
84 double const epsilon
= 1.0e-15; // relative error
85 bool bHasfound
= false;
87 // e_{-1} = 0; e_0 = alpha_0 / b_2
88 double u
; // u_0 = e_0/f_0 = alpha_0/m_0 = alpha_0
90 // first used with k=1
91 double m_bar
; // m_bar_k = m_k * f_bar_{k-1}
92 double g_bar
; // g_bar_k = m_bar_k - a_{k+1} + g_{k-1}
93 double g_bar_delta_u
; // g_bar_delta_u_k = f_bar_{k-1} * alpha_k
94 // - g_{k-1} * delta_u_{k-1} - m_bar_k * u_{k-1}
95 // f_{-1} = 0.0; f_0 = m_0 / b_2 = 1/(-1) = -1
96 double g
= 0.0; // g_0= f_{-1} / f_0 = 0/(-1) = 0
97 double delta_u
= 0.0; // dummy initialize, first used with * 0
98 double f_bar
= -1.0; // f_bar_k = 1/f_k, but only used for k=0
103 u
= 1.0; // u_0 = alpha_0
104 // k = 1.0; at least one step is necessary
105 // m_bar_k = m_k * f_bar_{k-1} ==> m_bar_1 = 0.0
106 g_bar_delta_u
= 0.0; // alpha_k = 0.0, m_bar = 0.0; g= 0.0
107 g_bar
= - 2.0/fX
; // k = 1.0, g = 0.0
108 delta_u
= g_bar_delta_u
/ g_bar
;
109 u
= u
+ delta_u
; // u_k = u_{k-1} + delta_u_k
110 g
= -1.0 / g_bar
; // g_k=b_{k+2}/g_bar_k
111 f_bar
= f_bar
* g
; // f_bar_k = f_bar_{k-1}* g_k
113 // From now on all alpha_k = 0.0 and k > N+1
116 { // N >= 1 and alpha_k = 0.0 for k<N
117 u
=0.0; // u_0 = alpha_0
118 for (k
=1.0; k
<= N
-1; k
= k
+ 1.0)
120 m_bar
=2.0 * fmod(k
-1.0, 2.0) * f_bar
;
121 g_bar_delta_u
= - g
* delta_u
- m_bar
* u
; // alpha_k = 0.0
122 g_bar
= m_bar
- 2.0*k
/fX
+ g
;
123 delta_u
= g_bar_delta_u
/ g_bar
;
128 // Step alpha_N = 1.0
129 m_bar
=2.0 * fmod(k
-1.0, 2.0) * f_bar
;
130 g_bar_delta_u
= f_bar
- g
* delta_u
- m_bar
* u
; // alpha_k = 1.0
131 g_bar
= m_bar
- 2.0*k
/fX
+ g
;
132 delta_u
= g_bar_delta_u
/ g_bar
;
138 // Loop until desired accuracy, always alpha_k = 0.0
141 m_bar
= 2.0 * fmod(k
-1.0, 2.0) * f_bar
;
142 g_bar_delta_u
= - g
* delta_u
- m_bar
* u
;
143 g_bar
= m_bar
- 2.0*k
/fX
+ g
;
144 delta_u
= g_bar_delta_u
/ g_bar
;
148 bHasfound
= (fabs(delta_u
)<=fabs(u
)*epsilon
);
151 while (!bHasfound
&& k
<= fMaxIteration
);
153 throw NoConvergenceException(); // unlikely to happen
162 /* The BESSEL function, first kind, modified:
165 I_n(x) = SUM TERM(n,k) with TERM(n,k) := --------------
168 No asymptotic approximation used, see issue 43040.
171 double BesselI( double x
, sal_Int32 n
)
173 const sal_Int32 nMaxIteration
= 2000;
174 const double fXHalf
= x
/ 2.0;
176 throw IllegalArgumentException();
178 double fResult
= 0.0;
180 /* Start the iteration without TERM(n,0), which is set here.
182 TERM(n,0) = (x/2)^n / n!
186 // avoid overflow in Fak(n)
187 for( nK
= 1; nK
<= n
; ++nK
)
189 fTerm
= fTerm
/ static_cast< double >( nK
) * fXHalf
;
191 fResult
= fTerm
; // Start result with TERM(n,0).
195 const double fEpsilon
= 1.0E-15;
198 /* Calculation of TERM(n,k) from TERM(n,k-1):
201 TERM(n,k) = --------------
204 (x/2)^2 (x/2)^(n+2(k-1))
205 = --------------------------
206 k (k-1)! (n+k) (n+k-1)!
208 (x/2)^2 (x/2)^(n+2(k-1))
209 = --------- * ------------------
210 k(n+k) (k-1)! (n+k-1)!
213 = -------- TERM(n,k-1)
216 fTerm
= fTerm
* fXHalf
/ static_cast<double>(nK
) * fXHalf
/ static_cast<double>(nK
+n
);
220 while( (fabs( fTerm
) > fabs(fResult
) * fEpsilon
) && (nK
< nMaxIteration
) );
226 /// @throws IllegalArgumentException
227 /// @throws NoConvergenceException
228 static double Besselk0( double fNum
)
234 double fNum2
= fNum
* 0.5;
235 double y
= fNum2
* fNum2
;
237 fRet
= -log( fNum2
) * BesselI( fNum
, 0 ) +
238 ( -0.57721566 + y
* ( 0.42278420 + y
* ( 0.23069756 + y
* ( 0.3488590e-1 +
239 y
* ( 0.262698e-2 + y
* ( 0.10750e-3 + y
* 0.74e-5 ) ) ) ) ) );
243 double y
= 2.0 / fNum
;
245 fRet
= exp( -fNum
) / sqrt( fNum
) * ( 1.25331414 + y
* ( -0.7832358e-1 +
246 y
* ( 0.2189568e-1 + y
* ( -0.1062446e-1 + y
* ( 0.587872e-2 +
247 y
* ( -0.251540e-2 + y
* 0.53208e-3 ) ) ) ) ) );
253 /// @throws IllegalArgumentException
254 /// @throws NoConvergenceException
255 static double Besselk1( double fNum
)
261 double fNum2
= fNum
* 0.5;
262 double y
= fNum2
* fNum2
;
264 fRet
= log( fNum2
) * BesselI( fNum
, 1 ) +
265 ( 1.0 + y
* ( 0.15443144 + y
* ( -0.67278579 + y
* ( -0.18156897 + y
* ( -0.1919402e-1 +
266 y
* ( -0.110404e-2 + y
* -0.4686e-4 ) ) ) ) ) )
271 double y
= 2.0 / fNum
;
273 fRet
= exp( -fNum
) / sqrt( fNum
) * ( 1.25331414 + y
* ( 0.23498619 +
274 y
* ( -0.3655620e-1 + y
* ( 0.1504268e-1 + y
* ( -0.780353e-2 +
275 y
* ( 0.325614e-2 + y
* -0.68245e-3 ) ) ) ) ) );
282 double BesselK( double fNum
, sal_Int32 nOrder
)
286 case 0: return Besselk0( fNum
);
287 case 1: return Besselk1( fNum
);
290 double fTox
= 2.0 / fNum
;
291 double fBkm
= Besselk0( fNum
);
292 double fBk
= Besselk1( fNum
);
294 for( sal_Int32 n
= 1 ; n
< nOrder
; n
++ )
296 const double fBkp
= fBkm
+ double( n
) * fTox
* fBk
;
310 /* The BESSEL function, second kind, unmodified:
311 The algorithm for order 0 and for order 1 follows
312 http://www.reference-global.com/isbn/978-3-11-020354-7
313 Numerical Mathematics 1 / Numerische Mathematik 1,
314 An algorithm-based introduction / Eine algorithmisch orientierte Einfuehrung
315 Deuflhard, Peter; Hohmann, Andreas
316 Berlin, New York (Walter de Gruyter) 2008
317 4. ueberarb. u. erw. Aufl. 2008
318 eBook ISBN: 978-3-11-020355-4
319 Chapter 6.3.2 , algorithm 6.24
320 The source is in German.
321 See #i31656# for a commented version of the implementation, attachment #desc6
322 https://bz.apache.org/ooo/attachment.cgi?id=63609
325 /// @throws IllegalArgumentException
326 /// @throws NoConvergenceException
327 static double Bessely0( double fX
)
330 throw IllegalArgumentException();
331 const double fMaxIteration
= 9000000.0; // should not be reached
332 if (fX
> 5.0e+6) // iteration is not considerable better then approximation
333 return sqrt(1/f_PI
/fX
)
334 *(rtl::math::sin(fX
)-rtl::math::cos(fX
));
335 const double epsilon
= 1.0e-15;
336 const double EulerGamma
= 0.57721566490153286060;
337 double alpha
= log(fX
/2.0)+EulerGamma
;
341 double g_bar_delta_u
= 0.0;
342 double g_bar
= -2.0 / fX
;
343 double delta_u
= g_bar_delta_u
/ g_bar
;
344 double g
= -1.0/g_bar
;
345 double f_bar
= -1 * g
;
347 double sign_alpha
= 1.0;
348 bool bHasFound
= false;
352 double km1mod2
= fmod(k
-1.0, 2.0);
353 double m_bar
= (2.0*km1mod2
) * f_bar
;
358 alpha
= sign_alpha
* (4.0/k
);
359 sign_alpha
= -sign_alpha
;
361 g_bar_delta_u
= f_bar
* alpha
- g
* delta_u
- m_bar
* u
;
362 g_bar
= m_bar
- (2.0*k
)/fX
+ g
;
363 delta_u
= g_bar_delta_u
/ g_bar
;
367 bHasFound
= (fabs(delta_u
)<=fabs(u
)*epsilon
);
370 while (!bHasFound
&& k
<fMaxIteration
);
372 throw NoConvergenceException(); // not likely to happen
376 // See #i31656# for a commented version of this implementation, attachment #desc6
377 // https://bz.apache.org/ooo/attachment.cgi?id=63609
378 /// @throws IllegalArgumentException
379 /// @throws NoConvergenceException
380 static double Bessely1( double fX
)
383 throw IllegalArgumentException();
384 const double fMaxIteration
= 9000000.0; // should not be reached
385 if (fX
> 5.0e+6) // iteration is not considerable better then approximation
386 return - sqrt(1/f_PI
/fX
)
387 *(rtl::math::sin(fX
)+rtl::math::cos(fX
));
388 const double epsilon
= 1.0e-15;
389 const double EulerGamma
= 0.57721566490153286060;
390 double alpha
= 1.0/fX
;
394 alpha
= 1.0 - EulerGamma
- log(fX
/2.0);
395 double g_bar_delta_u
= -alpha
;
396 double g_bar
= -2.0 / fX
;
397 double delta_u
= g_bar_delta_u
/ g_bar
;
399 double g
= -1.0/g_bar
;
401 double sign_alpha
= -1.0;
402 bool bHasFound
= false;
406 double km1mod2
= fmod(k
-1.0,2.0);
407 double m_bar
= (2.0*km1mod2
) * f_bar
;
408 double q
= (k
-1.0)/2.0;
409 if (km1mod2
== 0.0) // k is odd
411 alpha
= sign_alpha
* (1.0/q
+ 1.0/(q
+1.0));
412 sign_alpha
= -sign_alpha
;
416 g_bar_delta_u
= f_bar
* alpha
- g
* delta_u
- m_bar
* u
;
417 g_bar
= m_bar
- (2.0*k
)/fX
+ g
;
418 delta_u
= g_bar_delta_u
/ g_bar
;
422 bHasFound
= (fabs(delta_u
)<=fabs(u
)*epsilon
);
425 while (!bHasFound
&& k
<fMaxIteration
);
427 throw NoConvergenceException();
431 double BesselY( double fNum
, sal_Int32 nOrder
)
435 case 0: return Bessely0( fNum
);
436 case 1: return Bessely1( fNum
);
439 double fTox
= 2.0 / fNum
;
440 double fBym
= Bessely0( fNum
);
441 double fBy
= Bessely1( fNum
);
443 for( sal_Int32 n
= 1 ; n
< nOrder
; n
++ )
445 const double fByp
= double( n
) * fTox
* fBy
- fBym
;
455 } // namespace analysis
458 /* vim:set shiftwidth=4 softtabstop=4 expandtab: */