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