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