From 7cf842d1a4fd9faf11d5a56cf1eb4e3cc5777049 Mon Sep 17 00:00:00 2001 From: "Eike Rathke [er]" Date: Mon, 7 Jun 2010 16:27:28 +0200 Subject: calc55: #i40309# #i31656# #i103458# improved accuracy of Bessel functions; patch from --- scaddins/source/analysis/analysis.cxx | 8 +- scaddins/source/analysis/analysis.hxx | 9 +- scaddins/source/analysis/analysisadd.idl | 15 +- scaddins/source/analysis/analysisdefs.hxx | 2 + scaddins/source/analysis/bessel.cxx | 367 ++++++++++++++++++------------ scaddins/source/analysis/bessel.hxx | 9 +- scaddins/source/analysis/makefile.mk | 3 +- 7 files changed, 255 insertions(+), 158 deletions(-) (limited to 'scaddins/source/analysis') diff --git a/scaddins/source/analysis/analysis.cxx b/scaddins/source/analysis/analysis.cxx index bbc58ed27c9c..0d206d0286dd 100644 --- a/scaddins/source/analysis/analysis.cxx +++ b/scaddins/source/analysis/analysis.cxx @@ -915,21 +915,21 @@ double SAL_CALL AnalysisAddIn::getLcm( constREFXPS& xOpt, const SEQSEQ( double ) } -double SAL_CALL AnalysisAddIn::getBesseli( double fNum, sal_Int32 nOrder ) THROWDEF_RTE_IAE +double SAL_CALL AnalysisAddIn::getBesseli( double fNum, sal_Int32 nOrder ) THROWDEF_RTE_IAE_NCE { double fRet = sca::analysis::BesselI( fNum, nOrder ); RETURN_FINITE( fRet ); } -double SAL_CALL AnalysisAddIn::getBesselj( double fNum, sal_Int32 nOrder ) THROWDEF_RTE_IAE +double SAL_CALL AnalysisAddIn::getBesselj( double fNum, sal_Int32 nOrder ) THROWDEF_RTE_IAE_NCE { double fRet = sca::analysis::BesselJ( fNum, nOrder ); RETURN_FINITE( fRet ); } -double SAL_CALL AnalysisAddIn::getBesselk( double fNum, sal_Int32 nOrder ) THROWDEF_RTE_IAE +double SAL_CALL AnalysisAddIn::getBesselk( double fNum, sal_Int32 nOrder ) THROWDEF_RTE_IAE_NCE { if( nOrder < 0 || fNum <= 0.0 ) THROW_IAE; @@ -939,7 +939,7 @@ double SAL_CALL AnalysisAddIn::getBesselk( double fNum, sal_Int32 nOrder ) THROW } -double SAL_CALL AnalysisAddIn::getBessely( double fNum, sal_Int32 nOrder ) THROWDEF_RTE_IAE +double SAL_CALL AnalysisAddIn::getBessely( double fNum, sal_Int32 nOrder ) THROWDEF_RTE_IAE_NCE { if( nOrder < 0 || fNum <= 0.0 ) THROW_IAE; diff --git a/scaddins/source/analysis/analysis.hxx b/scaddins/source/analysis/analysis.hxx index fbf7dd012050..ff789838b299 100644 --- a/scaddins/source/analysis/analysis.hxx +++ b/scaddins/source/analysis/analysis.hxx @@ -36,6 +36,7 @@ #include #include #include +#include #include // helper for implementations @@ -143,10 +144,10 @@ public: virtual double SAL_CALL getGcd( constREFXPS& xOpt, const SEQSEQ( double )& aVLst, const SEQ( ANY )& aOptVLst ) THROWDEF_RTE_IAE; virtual double SAL_CALL getLcm( constREFXPS& xOpt, const SEQSEQ( double )& aVLst, const SEQ( ANY )& aOptVLst ) THROWDEF_RTE_IAE; - virtual double SAL_CALL getBesseli( double fNum, sal_Int32 nOrder ) THROWDEF_RTE_IAE; - virtual double SAL_CALL getBesselj( double fNum, sal_Int32 nOrder ) THROWDEF_RTE_IAE; - virtual double SAL_CALL getBesselk( double fNum, sal_Int32 nOrder ) THROWDEF_RTE_IAE; - virtual double SAL_CALL getBessely( double fNum, sal_Int32 nOrder ) THROWDEF_RTE_IAE; + virtual double SAL_CALL getBesseli( double fNum, sal_Int32 nOrder ) THROWDEF_RTE_IAE_NCE; + virtual double SAL_CALL getBesselj( double fNum, sal_Int32 nOrder ) THROWDEF_RTE_IAE_NCE; + virtual double SAL_CALL getBesselk( double fNum, sal_Int32 nOrder ) THROWDEF_RTE_IAE_NCE; + virtual double SAL_CALL getBessely( double fNum, sal_Int32 nOrder ) THROWDEF_RTE_IAE_NCE; virtual STRING SAL_CALL getBin2Oct( constREFXPS& xOpt, const STRING& aNum, const ANY& rPlaces ) THROWDEF_RTE_IAE; virtual double SAL_CALL getBin2Dec( const STRING& aNum ) THROWDEF_RTE_IAE; diff --git a/scaddins/source/analysis/analysisadd.idl b/scaddins/source/analysis/analysisadd.idl index 14b47ac256fc..6245bff9ad5a 100644 --- a/scaddins/source/analysis/analysisadd.idl +++ b/scaddins/source/analysis/analysisadd.idl @@ -25,6 +25,9 @@ * ************************************************************************/ +#ifndef __com_sun_star_sheet_NoConvergenceException_idl__ +#include +#endif #include #include @@ -145,19 +148,23 @@ module addin /// besseli. double getBesseli( [in] double Num, [in] long Order ) - raises( com::sun::star::lang::IllegalArgumentException ); + raises( com::sun::star::lang::IllegalArgumentException, + com::sun::star::sheet::NoConvergenceException ); /// besselj. double getBesselj( [in] double Num, [in] long Order ) - raises( com::sun::star::lang::IllegalArgumentException ); + raises( com::sun::star::lang::IllegalArgumentException, + com::sun::star::sheet::NoConvergenceException ); /// besselk. double getBesselk( [in] double Num, [in] long Order ) - raises( com::sun::star::lang::IllegalArgumentException ); + raises( com::sun::star::lang::IllegalArgumentException, + com::sun::star::sheet::NoConvergenceException ); /// bessely. double getBessely( [in] double Num, [in] long Order ) - raises( com::sun::star::lang::IllegalArgumentException ); + raises( com::sun::star::lang::IllegalArgumentException, + com::sun::star::sheet::NoConvergenceException ); /// bin2oct. string getBin2Oct( diff --git a/scaddins/source/analysis/analysisdefs.hxx b/scaddins/source/analysis/analysisdefs.hxx index 81f7435abc19..dae4205e2d84 100644 --- a/scaddins/source/analysis/analysisdefs.hxx +++ b/scaddins/source/analysis/analysisdefs.hxx @@ -44,6 +44,8 @@ #define THROW_RTE throw CSS::uno::RuntimeException() #define THROWDEF_RTE_IAE throw(CSS::uno::RuntimeException,CSS::lang::IllegalArgumentException) #define THROW_IAE throw CSS::lang::IllegalArgumentException() +#define THROWDEF_RTE_IAE_NCE throw(CSS::uno::RuntimeException,CSS::lang::IllegalArgumentException,CSS::sheet::NoConvergenceException) +#define THROW_NCE throw CSS::sheet::NoConvergenceException() #define CHK_Freq ( nFreq != 1 && nFreq != 2 && nFreq != 4 ) #define CHK_FINITE(d) if( !::rtl::math::isFinite( d ) ) THROW_IAE diff --git a/scaddins/source/analysis/bessel.cxx b/scaddins/source/analysis/bessel.cxx index 9b1f79b6d405..f853b49eb443 100644 --- a/scaddins/source/analysis/bessel.cxx +++ b/scaddins/source/analysis/bessel.cxx @@ -31,6 +31,7 @@ #include using ::com::sun::star::lang::IllegalArgumentException; +using ::com::sun::star::sheet::NoConvergenceException; namespace sca { namespace analysis { @@ -47,93 +48,129 @@ const double THRESHOLD = 30.0; // Threshold for usage of approximation form const double MAXEPSILON = 1e-10; // Maximum epsilon for end of iteration. const sal_Int32 MAXITER = 100; // Maximum number of iterations. - // ============================================================================ // BESSEL J // ============================================================================ /* The BESSEL function, first kind, unmodified: - - inf (-1)^k (x/2)^(n+2k) - J_n(x) = SUM TERM(n,k) with TERM(n,k) := --------------------- - k=0 k! (n+k)! - - Approximation for the BESSEL function, first kind, unmodified, for great x: - - J_n(x) ~ sqrt( 2 / (PI x) ) cos( x - n PI/2 - PI/4 ) for x>=0. - */ + The algorithm follows + http://www.reference-global.com/isbn/978-3-11-020354-7 + Numerical Mathematics 1 / Numerische Mathematik 1, + An algorithm-based introduction / Eine algorithmisch orientierte Einführung + Deuflhard, Peter; Hohmann, Andreas + Berlin, New York (Walter de Gruyter) 2008 + 4. überarb. u. erw. Aufl. 2008 + eBook ISBN: 978-3-11-020355-4 + Chapter 6.3.2 , algorithm 6.24 + The source is in German. + The BesselJ-function is a special case of the adjoint summation with + a_k = 2*(k-1)/x for k=1,... + b_k = -1, for all k, directly substituted + m_0=1, m_k=2 for k even, and m_k=0 for k odd, calculated on the fly + alpha_k=1 for k=N and alpha_k=0 otherwise +*/ // ---------------------------------------------------------------------------- -double BesselJ( double x, sal_Int32 n ) throw( IllegalArgumentException ) +double BesselJ( double x, sal_Int32 N ) throw (IllegalArgumentException, NoConvergenceException) + { - if( n < 0 ) + if( N < 0 ) throw IllegalArgumentException(); - - double fResult = 0.0; - if( fabs( x ) <= THRESHOLD ) + if (x==0.0) + return (N==0) ? 1.0 : 0.0; + + /* The algorithm works only for x>0, therefore remember sign. BesselJ + with integer order N is an even function for even N (means J(-x)=J(x)) + and an odd function for odd N (means J(-x)=-J(x)).*/ + double fSign = (N % 2 == 1 && x < 0) ? -1.0 : 1.0; + double fX = fabs(x); + + const double fMaxIteration = 9000000.0; //experimental, for to return in < 3 seconds + double fEstimateIteration = fX * 1.5 + N; + bool bAsymptoticPossible = pow(fX,0.4) > N; + if (fEstimateIteration > fMaxIteration) { - /* Start the iteration without TERM(n,0), which is set here. - - TERM(n,0) = (x/2)^n / n! - */ - double fTerm = pow( x / 2.0, (double)n ) / Fak( n ); - sal_Int32 nK = 1; // Start the iteration with k=1. - fResult = fTerm; // Start result with TERM(n,0). - - const double fSqrX = x * x / -4.0; + if (bAsymptoticPossible) + return fSign * sqrt(f_2_DIV_PI/fX)* cos(fX-N*f_PI_DIV_2-f_PI_DIV_4); + else + throw NoConvergenceException(); + } - do + double epsilon = 1.0e-15; // relative error + bool bHasfound = false; + double k= 0.0; + // e_{-1} = 0; e_0 = alpha_0 / b_2 + double u ; // u_0 = e_0/f_0 = alpha_0/m_0 = alpha_0 + + // first used with k=1 + double m_bar; // m_bar_k = m_k * f_bar_{k-1} + double g_bar; // g_bar_k = m_bar_k - a_{k+1} + g_{k-1} + double g_bar_delta_u; // g_bar_delta_u_k = f_bar_{k-1} * alpha_k + // - g_{k-1} * delta_u_{k-1} - m_bar_k * u_{k-1} + // f_{-1} = 0.0; f_0 = m_0 / b_2 = 1/(-1) = -1 + double g = 0.0; // g_0= f_{-1} / f_0 = 0/(-1) = 0 + double delta_u = 0.0; // dummy initialize, first used with * 0 + double f_bar = -1.0; // f_bar_k = 1/f_k, but only used for k=0 + + if (N==0) + { + //k=0; alpha_0 = 1.0 + u = 1.0; // u_0 = alpha_0 + // k = 1.0; at least one step is necessary + // m_bar_k = m_k * f_bar_{k-1} ==> m_bar_1 = 0.0 + g_bar_delta_u = 0.0; // alpha_k = 0.0, m_bar = 0.0; g= 0.0 + g_bar = - 2.0/fX; // k = 1.0, g = 0.0 + delta_u = g_bar_delta_u / g_bar; + u = u + delta_u ; // u_k = u_{k-1} + delta_u_k + g = -1.0 / g_bar; // g_k=b_{k+2}/g_bar_k + f_bar = f_bar * g; // f_bar_k = f_bar_{k-1}* g_k + k = 2.0; + // From now on all alpha_k = 0.0 and k > N+1 + } + else + { // N >= 1 and alpha_k = 0.0 for k MAXEPSILON) && (++nK < MAXITER) ); + // Step alpha_N = 1.0 + m_bar=2.0 * fmod(k-1.0, 2.0) * f_bar; + g_bar_delta_u = f_bar - g * delta_u - m_bar * u; // alpha_k = 1.0 + g_bar = m_bar - 2.0*k/fX + g; + delta_u = g_bar_delta_u / g_bar; + u = u + delta_u; + g = -1.0/g_bar; + f_bar = f_bar * g; + k = k + 1.0; } - else + // Loop until desired accuracy, always alpha_k = 0.0 + do { - /* Approximation for the BESSEL function, first kind, unmodified: - - J_n(x) ~ sqrt( 2 / (PI x) ) cos( x - n PI/2 - PI/4 ) for x>=0. - - The BESSEL function J_n with n IN {0,2,4,...} is axially symmetric at - x=0, means J_n(x) = J_n(-x). Therefore the approximation for x<0 is: - - J_n(x) = J_n(|x|) for x<0 and n IN {0,2,4,...}. - - The BESSEL function J_n with n IN {1,3,5,...} is point-symmetric at - x=0, means J_n(x) = -J_n(-x). Therefore the approximation for x<0 is: - - J_n(x) = -J_n(|x|) for x<0 and n IN {1,3,5,...}. - */ - double fXAbs = fabs( x ); - fResult = sqrt( f_2_DIV_PI / fXAbs ) * cos( fXAbs - n * f_PI_DIV_2 - f_PI_DIV_4 ); - if( (n & 1) && (x < 0.0) ) - fResult = -fResult; + m_bar = 2.0 * fmod(k-1.0, 2.0) * f_bar; + g_bar_delta_u = - g * delta_u - m_bar * u; + g_bar = m_bar - 2.0*k/fX + g; + delta_u = g_bar_delta_u / g_bar; + u = u + delta_u; + g = -1.0/g_bar; + f_bar = f_bar * g; + bHasfound = (fabs(delta_u)<=fabs(u)*epsilon); + k = k + 1.0; } - return fResult; + while (!bHasfound && k <= fMaxIteration); + if (bHasfound) + return u * fSign; + else + throw NoConvergenceException(); // unlikely to happen } - // ============================================================================ // BESSEL I // ============================================================================ @@ -151,7 +188,7 @@ double BesselJ( double x, sal_Int32 n ) throw( IllegalArgumentException ) // ---------------------------------------------------------------------------- -double BesselI( double x, sal_Int32 n ) throw( IllegalArgumentException ) +double BesselI( double x, sal_Int32 n ) throw( IllegalArgumentException, NoConvergenceException ) { if( n < 0 ) throw IllegalArgumentException(); @@ -222,7 +259,7 @@ double BesselI( double x, sal_Int32 n ) throw( IllegalArgumentException ) // ============================================================================ -double Besselk0( double fNum ) throw( IllegalArgumentException ) +double Besselk0( double fNum ) throw( IllegalArgumentException, NoConvergenceException ) { double fRet; @@ -248,7 +285,7 @@ double Besselk0( double fNum ) throw( IllegalArgumentException ) } -double Besselk1( double fNum ) throw( IllegalArgumentException ) +double Besselk1( double fNum ) throw( IllegalArgumentException, NoConvergenceException ) { double fRet; @@ -275,7 +312,7 @@ double Besselk1( double fNum ) throw( IllegalArgumentException ) } -double BesselK( double fNum, sal_Int32 nOrder ) throw( IllegalArgumentException ) +double BesselK( double fNum, sal_Int32 nOrder ) throw( IllegalArgumentException, NoConvergenceException ) { switch( nOrder ) { @@ -301,87 +338,136 @@ double BesselK( double fNum, sal_Int32 nOrder ) throw( IllegalArgumentException } } +// ============================================================================ +// BESSEL Y +// ============================================================================ -double Bessely0( double fNum ) throw( IllegalArgumentException ) +/* The BESSEL function, second kind, unmodified: + The algorithm for order 0 and for order 1 follows + http://www.reference-global.com/isbn/978-3-11-020354-7 + Numerical Mathematics 1 / Numerische Mathematik 1, + An algorithm-based introduction / Eine algorithmisch orientierte Einführung + Deuflhard, Peter; Hohmann, Andreas + Berlin, New York (Walter de Gruyter) 2008 + 4. überarb. u. erw. Aufl. 2008 + eBook ISBN: 978-3-11-020355-4 + Chapter 6.3.2 , algorithm 6.24 + The source is in German. + See #i31656# for a commented version of the implementation, attachment #desc6 + http://www.openoffice.org/nonav/issues/showattachment.cgi/63609/Comments%20to%20the%20implementation%20of%20the%20Bessel%20functions.odt +*/ + +double Bessely0( double fX ) throw( IllegalArgumentException, NoConvergenceException ) { - double fRet; - - if( fNum < 8.0 ) + if (fX <= 0) + throw IllegalArgumentException(); + const double fMaxIteration = 9000000.0; // should not be reached + if (fX > 5.0e+6) // iteration is not considerable better then approximation + return sqrt(1/f_PI/fX) + *(rtl::math::sin(fX)-rtl::math::cos(fX)); + const double epsilon = 1.0e-15; + const double EulerGamma = 0.57721566490153286060; + double alpha = log(fX/2.0)+EulerGamma; + double u = alpha; + + double k = 1.0; + double m_bar = 0.0; + double g_bar_delta_u = 0.0; + double g_bar = -2.0 / fX; + double delta_u = g_bar_delta_u / g_bar; + double g = -1.0/g_bar; + double f_bar = -1 * g; + + double sign_alpha = 1.0; + double km1mod2; + bool bHasFound = false; + k = k + 1; + do { - double y = fNum * fNum; - - double f1 = -2957821389.0 + y * ( 7062834065.0 + y * ( -512359803.6 + - y * ( 10879881.29 + y * ( -86327.92757 + y * 228.4622733 ) ) ) ); - - double f2 = 40076544269.0 + y * ( 745249964.8 + y * ( 7189466.438 + - y * ( 47447.26470 + y * ( 226.1030244 + y ) ) ) ); - - fRet = f1 / f2 + 0.636619772 * BesselJ( fNum, 0 ) * log( fNum ); + km1mod2 = fmod(k-1.0,2.0); + m_bar=(2.0*km1mod2) * f_bar; + if (km1mod2 == 0.0) + alpha = 0.0; + else + { + alpha = sign_alpha * (4.0/k); + sign_alpha = -sign_alpha; + } + g_bar_delta_u = f_bar * alpha - g * delta_u - m_bar * u; + g_bar = m_bar - (2.0*k)/fX + g; + delta_u = g_bar_delta_u / g_bar; + u = u+delta_u; + g = -1.0 / g_bar; + f_bar = f_bar*g; + bHasFound = (fabs(delta_u)<=fabs(u)*epsilon); + k=k+1; } + while (!bHasFound && k 5.0e+6) // iteration is not considerable better then approximation + return - sqrt(1/f_PI/fX) + *(rtl::math::sin(fX)+rtl::math::cos(fX)); + const double epsilon = 1.0e-15; + const double EulerGamma = 0.57721566490153286060; + double alpha = 1.0/fX; + double f_bar = -1.0; + double g = 0.0; + double u = alpha; + double k = 1.0; + double m_bar = 0.0; + alpha = 1.0 - EulerGamma - log(fX/2.0); + double g_bar_delta_u = -alpha; + double g_bar = -2.0 / fX; + double delta_u = g_bar_delta_u / g_bar; + u = u + delta_u; + g = -1.0/g_bar; + f_bar = f_bar * g; + double sign_alpha = -1.0; + double km1mod2; //will be (k-1) mod 2 + double q; // will be (k-1) div 2 + bool bHasFound = false; + k = k + 1.0; + do { - double y = fNum * fNum; - - double f1 = fNum * ( -0.4900604943e13 + y * ( 0.1275274390e13 + - y * ( -0.5153438139e11 + y * ( 0.7349264551e9 + - y * ( -0.4237922726e7 + y * 0.8511937935e4 ) ) ) ) ); - - double f2 = 0.2499580570e14 + y * ( 0.4244419664e12 + - y * ( 0.3733650367e10 + y * ( 0.2245904002e8 + - y * ( 0.1020426050e6 + y * ( 0.3549632885e3 + y ) ) ) ) ); - - fRet = f1 / f2 + 0.636619772 * ( BesselJ( fNum, 1 ) * log( fNum ) - 1.0 / fNum ); + km1mod2 = fmod(k-1.0,2.0); + m_bar=(2.0*km1mod2) * f_bar; + q = (k-1.0)/2.0; + if (km1mod2 == 0.0) // k is odd + { + alpha = sign_alpha * (1.0/q + 1.0/(q+1.0)); + sign_alpha = -sign_alpha; + } + else + alpha = 0.0; + g_bar_delta_u = f_bar * alpha - g * delta_u - m_bar * u; + g_bar = m_bar - (2.0*k)/fX + g; + delta_u = g_bar_delta_u / g_bar; + u = u+delta_u; + g = -1.0 / g_bar; + f_bar = f_bar*g; + bHasFound = (fabs(delta_u)<=fabs(u)*epsilon); + k=k+1; } + while (!bHasFound && k +#include #include namespace sca { @@ -39,16 +40,16 @@ namespace analysis { // ============================================================================ /** Returns the result for the unmodified BESSEL function of first kind (J), n-th order, at point x. */ -double BesselJ( double x, sal_Int32 n ) throw( ::com::sun::star::lang::IllegalArgumentException ); +double BesselJ( double x, sal_Int32 n ) throw( ::com::sun::star::lang::IllegalArgumentException, ::com::sun::star::sheet::NoConvergenceException ); /** Returns the result for the modified BESSEL function of first kind (I), n-th order, at point x. */ -double BesselI( double x, sal_Int32 n ) throw( ::com::sun::star::lang::IllegalArgumentException ); +double BesselI( double x, sal_Int32 n ) throw( ::com::sun::star::lang::IllegalArgumentException, ::com::sun::star::sheet::NoConvergenceException ); /** Returns the result for the unmodified BESSEL function of second kind (Y), n-th order, at point x. */ -double BesselY( double x, sal_Int32 n ) throw( ::com::sun::star::lang::IllegalArgumentException ); +double BesselY( double x, sal_Int32 n ) throw( ::com::sun::star::lang::IllegalArgumentException, ::com::sun::star::sheet::NoConvergenceException ); /** Returns the result for the modified BESSEL function of second kind (K), n-th order, at point x. */ -double BesselK( double x, sal_Int32 n ) throw( ::com::sun::star::lang::IllegalArgumentException ); +double BesselK( double x, sal_Int32 n ) throw( ::com::sun::star::lang::IllegalArgumentException, ::com::sun::star::sheet::NoConvergenceException ); // ============================================================================ diff --git a/scaddins/source/analysis/makefile.mk b/scaddins/source/analysis/makefile.mk index d837d644dc21..a23f9886b9b1 100644 --- a/scaddins/source/analysis/makefile.mk +++ b/scaddins/source/analysis/makefile.mk @@ -73,7 +73,8 @@ UNOTYPES=\ com.sun.star.uno.XComponentContext \ com.sun.star.util.Date \ com.sun.star.util.XNumberFormatter \ - com.sun.star.util.XNumberFormatTypes + com.sun.star.util.XNumberFormatTypes \ + com.sun.star.sheet.NoConvergenceException # --- Files ------------------------------------- -- cgit