1<html> 2<head> 3<meta http-equiv="Content-Type" content="text/html; charset=UTF-8"> 4<title>Error Function Inverses</title> 5<link rel="stylesheet" href="../../math.css" type="text/css"> 6<meta name="generator" content="DocBook XSL Stylesheets V1.79.1"> 7<link rel="home" href="../../index.html" title="Math Toolkit 2.12.0"> 8<link rel="up" href="../sf_erf.html" title="Error Functions"> 9<link rel="prev" href="error_function.html" title="Error Function erf and complement erfc"> 10<link rel="next" href="../sf_poly.html" title="Polynomials"> 11</head> 12<body bgcolor="white" text="black" link="#0000FF" vlink="#840084" alink="#0000FF"> 13<table cellpadding="2" width="100%"><tr> 14<td valign="top"><img alt="Boost C++ Libraries" width="277" height="86" src="../../../../../../boost.png"></td> 15<td align="center"><a href="../../../../../../index.html">Home</a></td> 16<td align="center"><a href="../../../../../../libs/libraries.htm">Libraries</a></td> 17<td align="center"><a href="http://www.boost.org/users/people.html">People</a></td> 18<td align="center"><a href="http://www.boost.org/users/faq.html">FAQ</a></td> 19<td align="center"><a href="../../../../../../more/index.htm">More</a></td> 20</tr></table> 21<hr> 22<div class="spirit-nav"> 23<a accesskey="p" href="error_function.html"><img src="../../../../../../doc/src/images/prev.png" alt="Prev"></a><a accesskey="u" href="../sf_erf.html"><img src="../../../../../../doc/src/images/up.png" alt="Up"></a><a accesskey="h" href="../../index.html"><img src="../../../../../../doc/src/images/home.png" alt="Home"></a><a accesskey="n" href="../sf_poly.html"><img src="../../../../../../doc/src/images/next.png" alt="Next"></a> 24</div> 25<div class="section"> 26<div class="titlepage"><div><div><h3 class="title"> 27<a name="math_toolkit.sf_erf.error_inv"></a><a class="link" href="error_inv.html" title="Error Function Inverses">Error Function Inverses</a> 28</h3></div></div></div> 29<h5> 30<a name="math_toolkit.sf_erf.error_inv.h0"></a> 31 <span class="phrase"><a name="math_toolkit.sf_erf.error_inv.synopsis"></a></span><a class="link" href="error_inv.html#math_toolkit.sf_erf.error_inv.synopsis">Synopsis</a> 32 </h5> 33<pre class="programlisting"><span class="preprocessor">#include</span> <span class="special"><</span><span class="identifier">boost</span><span class="special">/</span><span class="identifier">math</span><span class="special">/</span><span class="identifier">special_functions</span><span class="special">/</span><span class="identifier">erf</span><span class="special">.</span><span class="identifier">hpp</span><span class="special">></span> 34</pre> 35<pre class="programlisting"><span class="keyword">namespace</span> <span class="identifier">boost</span><span class="special">{</span> <span class="keyword">namespace</span> <span class="identifier">math</span><span class="special">{</span> 36 37<span class="keyword">template</span> <span class="special"><</span><span class="keyword">class</span> <span class="identifier">T</span><span class="special">></span> 38<a class="link" href="../result_type.html" title="Calculation of the Type of the Result"><span class="emphasis"><em>calculated-result-type</em></span></a> <span class="identifier">erf_inv</span><span class="special">(</span><span class="identifier">T</span> <span class="identifier">p</span><span class="special">);</span> 39 40<span class="keyword">template</span> <span class="special"><</span><span class="keyword">class</span> <span class="identifier">T</span><span class="special">,</span> <span class="keyword">class</span> <a class="link" href="../../policy.html" title="Chapter 21. Policies: Controlling Precision, Error Handling etc">Policy</a><span class="special">></span> 41<a class="link" href="../result_type.html" title="Calculation of the Type of the Result"><span class="emphasis"><em>calculated-result-type</em></span></a> <span class="identifier">erf_inv</span><span class="special">(</span><span class="identifier">T</span> <span class="identifier">p</span><span class="special">,</span> <span class="keyword">const</span> <a class="link" href="../../policy.html" title="Chapter 21. Policies: Controlling Precision, Error Handling etc">Policy</a><span class="special">&);</span> 42 43<span class="keyword">template</span> <span class="special"><</span><span class="keyword">class</span> <span class="identifier">T</span><span class="special">></span> 44<a class="link" href="../result_type.html" title="Calculation of the Type of the Result"><span class="emphasis"><em>calculated-result-type</em></span></a> <span class="identifier">erfc_inv</span><span class="special">(</span><span class="identifier">T</span> <span class="identifier">p</span><span class="special">);</span> 45 46<span class="keyword">template</span> <span class="special"><</span><span class="keyword">class</span> <span class="identifier">T</span><span class="special">,</span> <span class="keyword">class</span> <a class="link" href="../../policy.html" title="Chapter 21. Policies: Controlling Precision, Error Handling etc">Policy</a><span class="special">></span> 47<a class="link" href="../result_type.html" title="Calculation of the Type of the Result"><span class="emphasis"><em>calculated-result-type</em></span></a> <span class="identifier">erfc_inv</span><span class="special">(</span><span class="identifier">T</span> <span class="identifier">p</span><span class="special">,</span> <span class="keyword">const</span> <a class="link" href="../../policy.html" title="Chapter 21. Policies: Controlling Precision, Error Handling etc">Policy</a><span class="special">&);</span> 48 49<span class="special">}}</span> <span class="comment">// namespaces</span> 50</pre> 51<p> 52 The return type of these functions is computed using the <a class="link" href="../result_type.html" title="Calculation of the Type of the Result"><span class="emphasis"><em>result 53 type calculation rules</em></span></a>: the return type is <code class="computeroutput"><span class="keyword">double</span></code> if T is an integer type, and T otherwise. 54 </p> 55<p> 56 The final <a class="link" href="../../policy.html" title="Chapter 21. Policies: Controlling Precision, Error Handling etc">Policy</a> argument is optional and can 57 be used to control the behaviour of the function: how it handles errors, 58 what level of precision to use etc. Refer to the <a class="link" href="../../policy.html" title="Chapter 21. Policies: Controlling Precision, Error Handling etc">policy 59 documentation for more details</a>. 60 </p> 61<h5> 62<a name="math_toolkit.sf_erf.error_inv.h1"></a> 63 <span class="phrase"><a name="math_toolkit.sf_erf.error_inv.description"></a></span><a class="link" href="error_inv.html#math_toolkit.sf_erf.error_inv.description">Description</a> 64 </h5> 65<pre class="programlisting"><span class="keyword">template</span> <span class="special"><</span><span class="keyword">class</span> <span class="identifier">T</span><span class="special">></span> 66<a class="link" href="../result_type.html" title="Calculation of the Type of the Result"><span class="emphasis"><em>calculated-result-type</em></span></a> <span class="identifier">erf_inv</span><span class="special">(</span><span class="identifier">T</span> <span class="identifier">z</span><span class="special">);</span> 67 68<span class="keyword">template</span> <span class="special"><</span><span class="keyword">class</span> <span class="identifier">T</span><span class="special">,</span> <span class="keyword">class</span> <a class="link" href="../../policy.html" title="Chapter 21. Policies: Controlling Precision, Error Handling etc">Policy</a><span class="special">></span> 69<a class="link" href="../result_type.html" title="Calculation of the Type of the Result"><span class="emphasis"><em>calculated-result-type</em></span></a> <span class="identifier">erf_inv</span><span class="special">(</span><span class="identifier">T</span> <span class="identifier">z</span><span class="special">,</span> <span class="keyword">const</span> <a class="link" href="../../policy.html" title="Chapter 21. Policies: Controlling Precision, Error Handling etc">Policy</a><span class="special">&);</span> 70</pre> 71<p> 72 Returns the <a href="http://functions.wolfram.com/GammaBetaErf/InverseErf/" target="_top">inverse 73 error function</a> of z, that is a value x such that: 74 </p> 75<div class="blockquote"><blockquote class="blockquote"><p> 76 <span class="serif_italic"><span class="emphasis"><em>p = erf(x);</em></span></span> 77 </p></blockquote></div> 78<div class="blockquote"><blockquote class="blockquote"><p> 79 <span class="inlinemediaobject"><img src="../../../graphs/erf_inv.svg" align="middle"></span> 80 81 </p></blockquote></div> 82<pre class="programlisting"><span class="keyword">template</span> <span class="special"><</span><span class="keyword">class</span> <span class="identifier">T</span><span class="special">></span> 83<a class="link" href="../result_type.html" title="Calculation of the Type of the Result"><span class="emphasis"><em>calculated-result-type</em></span></a> <span class="identifier">erfc_inv</span><span class="special">(</span><span class="identifier">T</span> <span class="identifier">z</span><span class="special">);</span> 84 85<span class="keyword">template</span> <span class="special"><</span><span class="keyword">class</span> <span class="identifier">T</span><span class="special">,</span> <span class="keyword">class</span> <a class="link" href="../../policy.html" title="Chapter 21. Policies: Controlling Precision, Error Handling etc">Policy</a><span class="special">></span> 86<a class="link" href="../result_type.html" title="Calculation of the Type of the Result"><span class="emphasis"><em>calculated-result-type</em></span></a> <span class="identifier">erfc_inv</span><span class="special">(</span><span class="identifier">T</span> <span class="identifier">z</span><span class="special">,</span> <span class="keyword">const</span> <a class="link" href="../../policy.html" title="Chapter 21. Policies: Controlling Precision, Error Handling etc">Policy</a><span class="special">&);</span> 87</pre> 88<p> 89 Returns the inverse of the complement of the error function of z, that is 90 a value x such that: 91 </p> 92<div class="blockquote"><blockquote class="blockquote"><p> 93 <span class="serif_italic"><span class="emphasis"><em>p = erfc(x);</em></span></span> 94 </p></blockquote></div> 95<div class="blockquote"><blockquote class="blockquote"><p> 96 <span class="inlinemediaobject"><img src="../../../graphs/erfc_inv.svg" align="middle"></span> 97 98 </p></blockquote></div> 99<h5> 100<a name="math_toolkit.sf_erf.error_inv.h2"></a> 101 <span class="phrase"><a name="math_toolkit.sf_erf.error_inv.accuracy"></a></span><a class="link" href="error_inv.html#math_toolkit.sf_erf.error_inv.accuracy">Accuracy</a> 102 </h5> 103<p> 104 For types up to and including 80-bit long doubles the approximations used 105 are accurate to less than ~ 2 epsilon. For higher precision types these functions 106 have the same accuracy as the <a class="link" href="error_function.html" title="Error Function erf and complement erfc">forward 107 error functions</a>. 108 </p> 109<div class="table"> 110<a name="math_toolkit.sf_erf.error_inv.table_erf_inv"></a><p class="title"><b>Table 8.30. Error rates for erf_inv</b></p> 111<div class="table-contents"><table class="table" summary="Error rates for erf_inv"> 112<colgroup> 113<col> 114<col> 115<col> 116<col> 117<col> 118</colgroup> 119<thead><tr> 120<th> 121 </th> 122<th> 123 <p> 124 GNU C++ version 7.1.0<br> linux<br> double 125 </p> 126 </th> 127<th> 128 <p> 129 GNU C++ version 7.1.0<br> linux<br> long double 130 </p> 131 </th> 132<th> 133 <p> 134 Sun compiler version 0x5150<br> Sun Solaris<br> long double 135 </p> 136 </th> 137<th> 138 <p> 139 Microsoft Visual C++ version 14.1<br> Win32<br> double 140 </p> 141 </th> 142</tr></thead> 143<tbody><tr> 144<td> 145 <p> 146 Inverse Erf Function 147 </p> 148 </td> 149<td> 150 <p> 151 <span class="blue">Max = 0ε (Mean = 0ε)</span> 152 </p> 153 </td> 154<td> 155 <p> 156 <span class="blue">Max = 0.996ε (Mean = 0.389ε)</span> 157 </p> 158 </td> 159<td> 160 <p> 161 <span class="blue">Max = 1.08ε (Mean = 0.395ε)</span> 162 </p> 163 </td> 164<td> 165 <p> 166 <span class="blue">Max = 1.09ε (Mean = 0.502ε)</span> 167 </p> 168 </td> 169</tr></tbody> 170</table></div> 171</div> 172<br class="table-break"><div class="table"> 173<a name="math_toolkit.sf_erf.error_inv.table_erfc_inv"></a><p class="title"><b>Table 8.31. Error rates for erfc_inv</b></p> 174<div class="table-contents"><table class="table" summary="Error rates for erfc_inv"> 175<colgroup> 176<col> 177<col> 178<col> 179<col> 180<col> 181</colgroup> 182<thead><tr> 183<th> 184 </th> 185<th> 186 <p> 187 GNU C++ version 7.1.0<br> linux<br> double 188 </p> 189 </th> 190<th> 191 <p> 192 GNU C++ version 7.1.0<br> linux<br> long double 193 </p> 194 </th> 195<th> 196 <p> 197 Sun compiler version 0x5150<br> Sun Solaris<br> long double 198 </p> 199 </th> 200<th> 201 <p> 202 Microsoft Visual C++ version 14.1<br> Win32<br> double 203 </p> 204 </th> 205</tr></thead> 206<tbody> 207<tr> 208<td> 209 <p> 210 Inverse Erfc Function 211 </p> 212 </td> 213<td> 214 <p> 215 <span class="blue">Max = 0ε (Mean = 0ε)</span> 216 </p> 217 </td> 218<td> 219 <p> 220 <span class="blue">Max = 0.996ε (Mean = 0.397ε)</span> 221 </p> 222 </td> 223<td> 224 <p> 225 <span class="blue">Max = 1.08ε (Mean = 0.403ε)</span> 226 </p> 227 </td> 228<td> 229 <p> 230 <span class="blue">Max = 1ε (Mean = 0.491ε)</span> 231 </p> 232 </td> 233</tr> 234<tr> 235<td> 236 <p> 237 Inverse Erfc Function: extreme values 238 </p> 239 </td> 240<td> 241 </td> 242<td> 243 <p> 244 <span class="blue">Max = 1.62ε (Mean = 0.383ε)</span> 245 </p> 246 </td> 247<td> 248 <p> 249 <span class="blue">Max = 1.62ε (Mean = 0.383ε)</span> 250 </p> 251 </td> 252<td> 253 </td> 254</tr> 255</tbody> 256</table></div> 257</div> 258<br class="table-break"><p> 259 The following error plot are based on an exhaustive search of the functions 260 domain, MSVC-15.5 at <code class="computeroutput"><span class="keyword">double</span></code> 261 precision, and GCC-7.1/Ubuntu for <code class="computeroutput"><span class="keyword">long</span> 262 <span class="keyword">double</span></code> and <code class="computeroutput"><span class="identifier">__float128</span></code>. 263 </p> 264<div class="blockquote"><blockquote class="blockquote"><p> 265 <span class="inlinemediaobject"><img src="../../../graphs/erfc__double.svg" align="middle"></span> 266 267 </p></blockquote></div> 268<div class="blockquote"><blockquote class="blockquote"><p> 269 <span class="inlinemediaobject"><img src="../../../graphs/erfc__80_bit_long_double.svg" align="middle"></span> 270 271 </p></blockquote></div> 272<div class="blockquote"><blockquote class="blockquote"><p> 273 <span class="inlinemediaobject"><img src="../../../graphs/erfc____float128.svg" align="middle"></span> 274 275 </p></blockquote></div> 276<h5> 277<a name="math_toolkit.sf_erf.error_inv.h3"></a> 278 <span class="phrase"><a name="math_toolkit.sf_erf.error_inv.testing"></a></span><a class="link" href="error_inv.html#math_toolkit.sf_erf.error_inv.testing">Testing</a> 279 </h5> 280<p> 281 There are two sets of tests: 282 </p> 283<div class="itemizedlist"><ul class="itemizedlist" style="list-style-type: disc; "> 284<li class="listitem"> 285 Basic sanity checks attempt to "round-trip" from <span class="emphasis"><em>x</em></span> 286 to <span class="emphasis"><em>p</em></span> and back again. These tests have quite generous 287 tolerances: in general both the error functions and their inverses change 288 so rapidly in some places that round tripping to more than a couple of 289 significant digits isn't possible. This is especially true when <span class="emphasis"><em>p</em></span> 290 is very near one: in this case there isn't enough "information content" 291 in the input to the inverse function to get back where you started. 292 </li> 293<li class="listitem"> 294 Accuracy checks using high-precision test values. These measure the accuracy 295 of the result, given <span class="emphasis"><em>exact</em></span> input values. 296 </li> 297</ul></div> 298<h5> 299<a name="math_toolkit.sf_erf.error_inv.h4"></a> 300 <span class="phrase"><a name="math_toolkit.sf_erf.error_inv.implementation"></a></span><a class="link" href="error_inv.html#math_toolkit.sf_erf.error_inv.implementation">Implementation</a> 301 </h5> 302<p> 303 These functions use a rational approximation <a class="link" href="../sf_implementation.html#math_toolkit.sf_implementation.rational_approximations_used">devised 304 by JM</a> to calculate an initial approximation to the result that is 305 accurate to ~10<sup>-19</sup>, then only if that has insufficient accuracy compared 306 to the epsilon for T, do we clean up the result using <a href="http://en.wikipedia.org/wiki/Simple_rational_approximation" target="_top">Halley 307 iteration</a>. 308 </p> 309<p> 310 Constructing rational approximations to the erf/erfc functions is actually 311 surprisingly hard, especially at high precision. For this reason no attempt 312 has been made to achieve 10<sup>-34 </sup> accuracy suitable for use with 128-bit reals. 313 </p> 314<p> 315 In the following discussion, <span class="emphasis"><em>p</em></span> is the value passed to 316 erf_inv, and <span class="emphasis"><em>q</em></span> is the value passed to erfc_inv, so that 317 <span class="emphasis"><em>p = 1 - q</em></span> and <span class="emphasis"><em>q = 1 - p</em></span> and in 318 both cases we want to solve for the same result <span class="emphasis"><em>x</em></span>. 319 </p> 320<p> 321 For <span class="emphasis"><em>p < 0.5</em></span> the inverse erf function is reasonably 322 smooth and the approximation: 323 </p> 324<div class="blockquote"><blockquote class="blockquote"><p> 325 <span class="serif_italic"><span class="emphasis"><em>x = p(p + 10)(Y + R(p))</em></span></span> 326 </p></blockquote></div> 327<p> 328 Gives a good result for a constant Y, and R(p) optimised for low absolute 329 error compared to |Y|. 330 </p> 331<p> 332 For q < 0.5 things get trickier, over the interval <span class="emphasis"><em>0.5 > 333 q > 0.25</em></span> the following approximation works well: 334 </p> 335<div class="blockquote"><blockquote class="blockquote"><p> 336 <span class="serif_italic"><span class="emphasis"><em>x = sqrt(-2log(q)) / (Y + R(q))</em></span></span> 337 </p></blockquote></div> 338<p> 339 While for q < 0.25, let 340 </p> 341<div class="blockquote"><blockquote class="blockquote"><p> 342 <span class="serif_italic"><span class="emphasis"><em>z = sqrt(-log(q))</em></span></span> 343 </p></blockquote></div> 344<p> 345 Then the result is given by: 346 </p> 347<div class="blockquote"><blockquote class="blockquote"><p> 348 <span class="serif_italic"><span class="emphasis"><em>x = z(Y + R(z - B))</em></span></span> 349 </p></blockquote></div> 350<p> 351 As before Y is a constant and the rational function R is optimised for low 352 absolute error compared to |Y|. B is also a constant: it is the smallest 353 value of <span class="emphasis"><em>z</em></span> for which each approximation is valid. There 354 are several approximations of this form each of which reaches a little further 355 into the tail of the erfc function (at <code class="computeroutput"><span class="keyword">long</span> 356 <span class="keyword">double</span></code> precision the extended exponent 357 range compared to <code class="computeroutput"><span class="keyword">double</span></code> means 358 that the tail goes on for a very long way indeed). 359 </p> 360</div> 361<table xmlns:rev="http://www.cs.rpi.edu/~gregod/boost/tools/doc/revision" width="100%"><tr> 362<td align="left"></td> 363<td align="right"><div class="copyright-footer">Copyright © 2006-2019 Nikhar 364 Agrawal, Anton Bikineev, Paul A. Bristow, Marco Guazzone, Christopher Kormanyos, 365 Hubert Holin, Bruno Lalande, John Maddock, Jeremy Murphy, Matthew Pulver, Johan 366 Råde, Gautam Sewani, Benjamin Sobotta, Nicholas Thompson, Thijs van den Berg, 367 Daryle Walker and Xiaogang Zhang<p> 368 Distributed under the Boost Software License, Version 1.0. (See accompanying 369 file LICENSE_1_0.txt or copy at <a href="http://www.boost.org/LICENSE_1_0.txt" target="_top">http://www.boost.org/LICENSE_1_0.txt</a>) 370 </p> 371</div></td> 372</tr></table> 373<hr> 374<div class="spirit-nav"> 375<a accesskey="p" href="error_function.html"><img src="../../../../../../doc/src/images/prev.png" alt="Prev"></a><a accesskey="u" href="../sf_erf.html"><img src="../../../../../../doc/src/images/up.png" alt="Up"></a><a accesskey="h" href="../../index.html"><img src="../../../../../../doc/src/images/home.png" alt="Home"></a><a accesskey="n" href="../sf_poly.html"><img src="../../../../../../doc/src/images/next.png" alt="Next"></a> 376</div> 377</body> 378</html> 379
