bump product version to 6.3.0.0.beta1
[LibreOffice.git] / scaddins / source / analysis / bessel.cxx
blobab87ca0e10cea62aaee09d7cf335c7888d77d7dc
1 /* -*- Mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
2 /*
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 .
20 #include "bessel.hxx"
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;
30 namespace sca {
31 namespace analysis {
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;
39 // BESSEL J
42 /* The BESSEL function, first kind, unmodified:
43 The algorithm follows
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 )
63 if( N < 0 )
64 throw IllegalArgumentException();
65 if (x==0.0)
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;
72 double fX = fabs(x);
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;
86 double k= 0.0;
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
100 if (N==0)
102 //k=0; alpha_0 = 1.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
112 k = 2.0;
113 // From now on all alpha_k = 0.0 and k > N+1
115 else
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;
124 u = u + delta_u;
125 g = -1.0/g_bar;
126 f_bar=f_bar * g;
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;
133 u = u + delta_u;
134 g = -1.0/g_bar;
135 f_bar = f_bar * g;
136 k = k + 1.0;
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;
145 u = u + delta_u;
146 g = -1.0/g_bar;
147 f_bar = f_bar * g;
148 bHasfound = (fabs(delta_u)<=fabs(u)*epsilon);
149 k = k + 1.0;
151 while (!bHasfound && k <= fMaxIteration);
152 if (!bHasfound)
153 throw NoConvergenceException(); // unlikely to happen
155 return u * fSign;
159 // BESSEL I
162 /* The BESSEL function, first kind, modified:
164 inf (x/2)^(n+2k)
165 I_n(x) = SUM TERM(n,k) with TERM(n,k) := --------------
166 k=0 k! (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;
175 if( n < 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!
184 sal_Int32 nK = 0;
185 double fTerm = 1.0;
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).
192 if( fTerm != 0.0 )
194 nK = 1;
195 const double fEpsilon = 1.0E-15;
198 /* Calculation of TERM(n,k) from TERM(n,k-1):
200 (x/2)^(n+2k)
201 TERM(n,k) = --------------
202 k! (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)!
212 x^2/4
213 = -------- TERM(n,k-1)
214 k(n+k)
216 fTerm = fTerm * fXHalf / static_cast<double>(nK) * fXHalf / static_cast<double>(nK+n);
217 fResult += fTerm;
218 nK++;
220 while( (fabs( fTerm ) > fabs(fResult) * fEpsilon) && (nK < nMaxIteration) );
223 return fResult;
226 /// @throws IllegalArgumentException
227 /// @throws NoConvergenceException
228 static double Besselk0( double fNum )
230 double fRet;
232 if( fNum <= 2.0 )
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 ) ) ) ) ) );
241 else
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 ) ) ) ) ) );
250 return fRet;
253 /// @throws IllegalArgumentException
254 /// @throws NoConvergenceException
255 static double Besselk1( double fNum )
257 double fRet;
259 if( fNum <= 2.0 )
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 ) ) ) ) ) )
267 / fNum;
269 else
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 ) ) ) ) ) );
278 return fRet;
282 double BesselK( double fNum, sal_Int32 nOrder )
284 switch( nOrder )
286 case 0: return Besselk0( fNum );
287 case 1: return Besselk1( fNum );
288 default:
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;
297 fBkm = fBk;
298 fBk = fBkp;
301 return fBk;
307 // BESSEL Y
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 )
329 if (fX <= 0)
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;
338 double u = alpha;
340 double k = 1.0;
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;
349 k = k + 1;
352 double km1mod2 = fmod(k-1.0, 2.0);
353 double m_bar = (2.0*km1mod2) * f_bar;
354 if (km1mod2 == 0.0)
355 alpha = 0.0;
356 else
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;
364 u = u+delta_u;
365 g = -1.0 / g_bar;
366 f_bar = f_bar*g;
367 bHasFound = (fabs(delta_u)<=fabs(u)*epsilon);
368 k=k+1;
370 while (!bHasFound && k<fMaxIteration);
371 if (!bHasFound)
372 throw NoConvergenceException(); // not likely to happen
373 return u*f_2_DIV_PI;
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 )
382 if (fX <= 0)
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;
391 double f_bar = -1.0;
392 double u = alpha;
393 double k = 1.0;
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;
398 u = u + delta_u;
399 double g = -1.0/g_bar;
400 f_bar = f_bar * g;
401 double sign_alpha = -1.0;
402 bool bHasFound = false;
403 k = k + 1.0;
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;
414 else
415 alpha = 0.0;
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;
419 u = u+delta_u;
420 g = -1.0 / g_bar;
421 f_bar = f_bar*g;
422 bHasFound = (fabs(delta_u)<=fabs(u)*epsilon);
423 k=k+1;
425 while (!bHasFound && k<fMaxIteration);
426 if (!bHasFound)
427 throw NoConvergenceException();
428 return -u*2.0/f_PI;
431 double BesselY( double fNum, sal_Int32 nOrder )
433 switch( nOrder )
435 case 0: return Bessely0( fNum );
436 case 1: return Bessely1( fNum );
437 default:
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;
446 fBym = fBy;
447 fBy = fByp;
450 return fBy;
455 } // namespace analysis
456 } // namespace sca
458 /* vim:set shiftwidth=4 softtabstop=4 expandtab: */