mirror of
				https://github.com/saitohirga/WSJT-X.git
				synced 2025-10-29 20:10:28 -04:00 
			
		
		
		
	
		
			
				
	
	
		
			137 lines
		
	
	
		
			4.1 KiB
		
	
	
	
		
			Plaintext
		
	
	
	
	
	
			
		
		
	
	
			137 lines
		
	
	
		
			4.1 KiB
		
	
	
	
		
			Plaintext
		
	
	
	
	
	
| [section:error_inv Error Function Inverses]
 | |
| 
 | |
| [h4 Synopsis]
 | |
| 
 | |
| ``
 | |
| #include <boost/math/special_functions/erf.hpp>
 | |
| ``
 | |
| 
 | |
|    namespace boost{ namespace math{
 | |
|    
 | |
|    template <class T>
 | |
|    ``__sf_result`` erf_inv(T p);
 | |
|    
 | |
|    template <class T, class ``__Policy``>
 | |
|    ``__sf_result`` erf_inv(T p, const ``__Policy``&);
 | |
|    
 | |
|    template <class T>
 | |
|    ``__sf_result`` erfc_inv(T p);
 | |
|    
 | |
|    template <class T, class ``__Policy``>
 | |
|    ``__sf_result`` erfc_inv(T p, const ``__Policy``&);
 | |
|    
 | |
|    }} // namespaces
 | |
|    
 | |
| The return type of these functions is computed using the __arg_promotion_rules:
 | |
| the return type is `double` if T is an integer type, and T otherwise.
 | |
| 
 | |
| [optional_policy]
 | |
| 
 | |
| [h4 Description]
 | |
| 
 | |
|    template <class T>
 | |
|    ``__sf_result`` erf_inv(T z);
 | |
|    
 | |
|    template <class T, class ``__Policy``>
 | |
|    ``__sf_result`` erf_inv(T z, const ``__Policy``&);
 | |
|    
 | |
| Returns the [@http://functions.wolfram.com/GammaBetaErf/InverseErf/ inverse error function]
 | |
| of z, that is a value x such that:
 | |
| 
 | |
|    p = erf(x);
 | |
| 
 | |
| [graph erf_inv]
 | |
| 
 | |
|    template <class T>
 | |
|    ``__sf_result`` erfc_inv(T z);
 | |
|    
 | |
|    template <class T, class ``__Policy``>
 | |
|    ``__sf_result`` erfc_inv(T z, const ``__Policy``&);
 | |
|    
 | |
| Returns the inverse of the complement of the error function of z, that is a
 | |
| value x such that:
 | |
| 
 | |
|    p = erfc(x);
 | |
| 
 | |
| [graph erfc_inv]
 | |
| 
 | |
| [h4 Accuracy]
 | |
| 
 | |
| For types up to and including 80-bit long doubles the approximations used
 | |
| are accurate to less than ~ 2 epsilon.  For higher precision types these 
 | |
| functions have the same accuracy as the 
 | |
| [link math_toolkit.sf_erf.error_function forward error functions].
 | |
| 
 | |
| [table_erf_inv]
 | |
| 
 | |
| [table_erfc_inv]
 | |
| 
 | |
| [h4 Testing]
 | |
| 
 | |
| There are two sets of tests: 
 | |
| 
 | |
| * Basic sanity checks attempt to "round-trip" from
 | |
| /x/ to /p/ and back again.  These tests have quite
 | |
| generous tolerances: in general both the error functions and their
 | |
| inverses change so rapidly in some places that round tripping to more than a couple
 | |
| of significant digits isn't possible.  This is especially true when
 | |
| /p/ is very near one: in this case there isn't enough 
 | |
| "information content" in the input to the inverse function to get
 | |
| back where you started.
 | |
| * Accuracy checks using high-precision test values.  These measure
 | |
| the accuracy of the result, given /exact/ input values.
 | |
| 
 | |
| [h4 Implementation]
 | |
| 
 | |
| These functions use a rational approximation [jm_rationals] 
 | |
| to calculate an initial
 | |
| approximation to the result that is accurate to ~10[super -19], 
 | |
| then only if that has insufficient accuracy compared to the epsilon for T,
 | |
| do we clean up the result using
 | |
| [@http://en.wikipedia.org/wiki/Simple_rational_approximation Halley iteration].
 | |
| 
 | |
| Constructing rational approximations to the erf/erfc functions is actually
 | |
| surprisingly hard, especially at high precision.  For this reason no attempt
 | |
| has been made to achieve 10[super -34 ] accuracy suitable for use with 128-bit
 | |
| reals.
 | |
| 
 | |
| In the following discussion, /p/ is the value passed to erf_inv, and /q/ is
 | |
| the value passed to erfc_inv, so that /p = 1 - q/ and /q = 1 - p/ and in both
 | |
| cases we want to solve for the same result /x/.
 | |
| 
 | |
| For /p < 0.5/ the inverse erf function is reasonably smooth and the approximation:
 | |
| 
 | |
|    x = p(p + 10)(Y + R(p))
 | |
|    
 | |
| Gives a good result for a constant Y, and R(p) optimised for low absolute error
 | |
| compared to |Y|.
 | |
| 
 | |
| For q < 0.5 things get trickier, over the interval /0.5 > q > 0.25/
 | |
| the following approximation works well:
 | |
| 
 | |
|    x = sqrt(-2log(q)) / (Y + R(q))
 | |
|    
 | |
| While for q < 0.25, let 
 | |
| 
 | |
|    z = sqrt(-log(q))
 | |
| 
 | |
| Then the result is given by:
 | |
| 
 | |
|    x = z(Y + R(z - B))
 | |
| 
 | |
| As before Y is a constant and the rational function R is optimised for low
 | |
| absolute error compared to |Y|.  B is also a constant: it is the smallest value
 | |
| of /z/ for which each approximation is valid.  There are several approximations
 | |
| of this form each of which reaches a little further into the tail of the erfc 
 | |
| function (at `long double` precision the extended exponent range compared to
 | |
| `double` means that the tail goes on for a very long way indeed).
 | |
| 
 | |
| [endsect][/ :error_inv The Error Function Inverses]
 | |
| 
 | |
| [/ 
 | |
|   Copyright 2006 John Maddock and Paul A. Bristow.
 | |
|   Distributed under the Boost Software License, Version 1.0.
 | |
|   (See accompanying file LICENSE_1_0.txt or copy at
 | |
|   http://www.boost.org/LICENSE_1_0.txt).
 | |
| ]
 |