📁 File Manager Pro
v10.0.3 | PHP: 8.1.34
Server: Apache
2026-06-22 02:18:38
📂
/ (Root)
/
opt
/
alt
/
ruby19
/
share
/
doc
/
ruby
/
html
/
dd
/
dcf
📍 /opt/alt/ruby19/share/doc/ruby/html/dd/dcf
🔄 Refresh
✏️
Editing: erf_8c_source.html
Read Only
<!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Transitional//EN" "http://www.w3.org/TR/xhtml1/DTD/xhtml1-transitional.dtd"> <html xmlns="http://www.w3.org/1999/xhtml"> <head> <meta http-equiv="Content-Type" content="text/xhtml;charset=UTF-8"/> <meta http-equiv="X-UA-Compatible" content="IE=9"/> <meta name="generator" content="Doxygen 1.8.14"/> <meta name="viewport" content="width=device-width, initial-scale=1"/> <title>Ruby: missing/erf.c Source File</title> <link href="../../tabs.css" rel="stylesheet" type="text/css"/> <script type="text/javascript" src="../../jquery.js"></script> <script type="text/javascript" src="../../dynsections.js"></script> <link href="../../doxygen.css" rel="stylesheet" type="text/css" /> </head> <body> <div id="top"><!-- do not remove this div, it is closed by doxygen! --> <div id="titlearea"> <table cellspacing="0" cellpadding="0"> <tbody> <tr style="height: 56px;"> <td id="projectalign" style="padding-left: 0.5em;"> <div id="projectname">Ruby  <span id="projectnumber">1.9.3p551(2014-11-13revision48407)</span> </div> </td> </tr> </tbody> </table> </div> <!-- end header part --> <!-- Generated by Doxygen 1.8.14 --> <script type="text/javascript" src="../../menudata.js"></script> <script type="text/javascript" src="../../menu.js"></script> <script type="text/javascript"> /* @license magnet:?xt=urn:btih:cf05388f2679ee054f2beb29a391d25f4e673ac3&dn=gpl-2.0.txt GPL-v2 */ $(function() { initMenu('../../',false,false,'search.php','Search'); }); /* @license-end */</script> <div id="main-nav"></div> <div id="nav-path" class="navpath"> <ul> <li class="navelem"><a class="el" href="../../dir_f3bfeebb553c3f6ecfb19202628b4493.html">missing</a></li> </ul> </div> </div><!-- top --> <div class="header"> <div class="headertitle"> <div class="title">erf.c</div> </div> </div><!--header--> <div class="contents"> <a href="../../dd/dcf/erf_8c.html">Go to the documentation of this file.</a><div class="fragment"><div class="line"><a name="l00001"></a><span class="lineno"> 1</span> <span class="comment">/* erf.c - public domain implementation of error function erf(3m)</span></div><div class="line"><a name="l00002"></a><span class="lineno"> 2</span> <span class="comment"></span></div><div class="line"><a name="l00003"></a><span class="lineno"> 3</span> <span class="comment">reference - Haruhiko Okumura: C-gengo niyoru saishin algorithm jiten</span></div><div class="line"><a name="l00004"></a><span class="lineno"> 4</span> <span class="comment"> (New Algorithm handbook in C language) (Gijyutsu hyouron</span></div><div class="line"><a name="l00005"></a><span class="lineno"> 5</span> <span class="comment"> sha, Tokyo, 1991) p.227 [in Japanese] */</span></div><div class="line"><a name="l00006"></a><span class="lineno"> 6</span> <span class="preprocessor">#include "<a class="code" href="../../d3/d90/missing_8h.html">ruby/missing.h</a>"</span></div><div class="line"><a name="l00007"></a><span class="lineno"> 7</span> <span class="preprocessor">#include <stdio.h></span></div><div class="line"><a name="l00008"></a><span class="lineno"> 8</span> <span class="preprocessor">#include <math.h></span></div><div class="line"><a name="l00009"></a><span class="lineno"> 9</span> </div><div class="line"><a name="l00010"></a><span class="lineno"> 10</span> <span class="preprocessor">#ifdef _WIN32</span></div><div class="line"><a name="l00011"></a><span class="lineno"> 11</span> <span class="preprocessor"># include <float.h></span></div><div class="line"><a name="l00012"></a><span class="lineno"> 12</span> <span class="preprocessor"># if !defined __MINGW32__ || defined __NO_ISOCEXT</span></div><div class="line"><a name="l00013"></a><span class="lineno"> 13</span> <span class="preprocessor"># ifndef isnan</span></div><div class="line"><a name="l00014"></a><span class="lineno"> 14</span> <span class="preprocessor"># define isnan(x) _isnan(x)</span></div><div class="line"><a name="l00015"></a><span class="lineno"> 15</span> <span class="preprocessor"># endif</span></div><div class="line"><a name="l00016"></a><span class="lineno"> 16</span> <span class="preprocessor"># ifndef isinf</span></div><div class="line"><a name="l00017"></a><span class="lineno"> 17</span> <span class="preprocessor"># define isinf(x) (!_finite(x) && !_isnan(x))</span></div><div class="line"><a name="l00018"></a><span class="lineno"> 18</span> <span class="preprocessor"># endif</span></div><div class="line"><a name="l00019"></a><span class="lineno"> 19</span> <span class="preprocessor"># ifndef finite</span></div><div class="line"><a name="l00020"></a><span class="lineno"> 20</span> <span class="preprocessor"># define finite(x) _finite(x)</span></div><div class="line"><a name="l00021"></a><span class="lineno"> 21</span> <span class="preprocessor"># endif</span></div><div class="line"><a name="l00022"></a><span class="lineno"> 22</span> <span class="preprocessor"># endif</span></div><div class="line"><a name="l00023"></a><span class="lineno"> 23</span> <span class="preprocessor">#endif</span></div><div class="line"><a name="l00024"></a><span class="lineno"> 24</span> </div><div class="line"><a name="l00025"></a><span class="lineno"> 25</span> <span class="keyword">static</span> <span class="keywordtype">double</span> <a class="code" href="../../dd/dcf/erf_8c.html#aa2417338852b65efc862b7cab59d216a">q_gamma</a>(<span class="keywordtype">double</span>, <span class="keywordtype">double</span>, <span class="keywordtype">double</span>);</div><div class="line"><a name="l00026"></a><span class="lineno"> 26</span> </div><div class="line"><a name="l00027"></a><span class="lineno"> 27</span> <span class="comment">/* Incomplete gamma function</span></div><div class="line"><a name="l00028"></a><span class="lineno"> 28</span> <span class="comment"> 1 / Gamma(a) * Int_0^x exp(-t) t^(a-1) dt */</span></div><div class="line"><a name="l00029"></a><span class="lineno"><a class="line" href="../../dd/dcf/erf_8c.html#a2f5e1b011eda67f1fbd1b524b82dafb5"> 29</a></span> <span class="keyword">static</span> <span class="keywordtype">double</span> <a class="code" href="../../dd/dcf/erf_8c.html#a2f5e1b011eda67f1fbd1b524b82dafb5">p_gamma</a>(<span class="keywordtype">double</span> a, <span class="keywordtype">double</span> x, <span class="keywordtype">double</span> loggamma_a)</div><div class="line"><a name="l00030"></a><span class="lineno"> 30</span> {</div><div class="line"><a name="l00031"></a><span class="lineno"> 31</span>  <span class="keywordtype">int</span> k;</div><div class="line"><a name="l00032"></a><span class="lineno"> 32</span>  <span class="keywordtype">double</span> <a class="code" href="../../d8/d90/nkf_8c.html#a5ea5ac7abf5cce39283e422add1067d5">result</a>, term, previous;</div><div class="line"><a name="l00033"></a><span class="lineno"> 33</span> </div><div class="line"><a name="l00034"></a><span class="lineno"> 34</span>  <span class="keywordflow">if</span> (x >= 1 + a) <span class="keywordflow">return</span> 1 - <a class="code" href="../../dd/dcf/erf_8c.html#aa2417338852b65efc862b7cab59d216a">q_gamma</a>(a, x, loggamma_a);</div><div class="line"><a name="l00035"></a><span class="lineno"> 35</span>  <span class="keywordflow">if</span> (x == 0) <span class="keywordflow">return</span> 0;</div><div class="line"><a name="l00036"></a><span class="lineno"> 36</span>  <a class="code" href="../../d8/d90/nkf_8c.html#a5ea5ac7abf5cce39283e422add1067d5">result</a> = term = exp(a * log(x) - x - loggamma_a) / a;</div><div class="line"><a name="l00037"></a><span class="lineno"> 37</span>  <span class="keywordflow">for</span> (k = 1; k < 1000; k++) {</div><div class="line"><a name="l00038"></a><span class="lineno"> 38</span>  term *= x / (a + k);</div><div class="line"><a name="l00039"></a><span class="lineno"> 39</span>  previous = <a class="code" href="../../d8/d90/nkf_8c.html#a5ea5ac7abf5cce39283e422add1067d5">result</a>; <a class="code" href="../../d8/d90/nkf_8c.html#a5ea5ac7abf5cce39283e422add1067d5">result</a> += term;</div><div class="line"><a name="l00040"></a><span class="lineno"> 40</span>  <span class="keywordflow">if</span> (<a class="code" href="../../d8/d90/nkf_8c.html#a5ea5ac7abf5cce39283e422add1067d5">result</a> == previous) <span class="keywordflow">return</span> <a class="code" href="../../d8/d90/nkf_8c.html#a5ea5ac7abf5cce39283e422add1067d5">result</a>;</div><div class="line"><a name="l00041"></a><span class="lineno"> 41</span>  }</div><div class="line"><a name="l00042"></a><span class="lineno"> 42</span>  fprintf(stderr, <span class="stringliteral">"erf.c:%d:p_gamma() could not converge."</span>, __LINE__);</div><div class="line"><a name="l00043"></a><span class="lineno"> 43</span>  <span class="keywordflow">return</span> <a class="code" href="../../d8/d90/nkf_8c.html#a5ea5ac7abf5cce39283e422add1067d5">result</a>;</div><div class="line"><a name="l00044"></a><span class="lineno"> 44</span> }</div><div class="line"><a name="l00045"></a><span class="lineno"> 45</span> </div><div class="line"><a name="l00046"></a><span class="lineno"> 46</span> <span class="comment">/* Incomplete gamma function</span></div><div class="line"><a name="l00047"></a><span class="lineno"> 47</span> <span class="comment"> 1 / Gamma(a) * Int_x^inf exp(-t) t^(a-1) dt */</span></div><div class="line"><a name="l00048"></a><span class="lineno"><a class="line" href="../../dd/dcf/erf_8c.html#aa2417338852b65efc862b7cab59d216a"> 48</a></span> <span class="keyword">static</span> <span class="keywordtype">double</span> <a class="code" href="../../dd/dcf/erf_8c.html#aa2417338852b65efc862b7cab59d216a">q_gamma</a>(<span class="keywordtype">double</span> a, <span class="keywordtype">double</span> x, <span class="keywordtype">double</span> loggamma_a)</div><div class="line"><a name="l00049"></a><span class="lineno"> 49</span> {</div><div class="line"><a name="l00050"></a><span class="lineno"> 50</span>  <span class="keywordtype">int</span> k;</div><div class="line"><a name="l00051"></a><span class="lineno"> 51</span>  <span class="keywordtype">double</span> <a class="code" href="../../d8/d90/nkf_8c.html#a5ea5ac7abf5cce39283e422add1067d5">result</a>, w, temp, previous;</div><div class="line"><a name="l00052"></a><span class="lineno"> 52</span>  <span class="keywordtype">double</span> la = 1, lb = 1 + x - a; <span class="comment">/* Laguerre polynomial */</span></div><div class="line"><a name="l00053"></a><span class="lineno"> 53</span> </div><div class="line"><a name="l00054"></a><span class="lineno"> 54</span>  <span class="keywordflow">if</span> (x < 1 + a) <span class="keywordflow">return</span> 1 - <a class="code" href="../../dd/dcf/erf_8c.html#a2f5e1b011eda67f1fbd1b524b82dafb5">p_gamma</a>(a, x, loggamma_a);</div><div class="line"><a name="l00055"></a><span class="lineno"> 55</span>  w = exp(a * log(x) - x - loggamma_a);</div><div class="line"><a name="l00056"></a><span class="lineno"> 56</span>  <a class="code" href="../../d8/d90/nkf_8c.html#a5ea5ac7abf5cce39283e422add1067d5">result</a> = w / lb;</div><div class="line"><a name="l00057"></a><span class="lineno"> 57</span>  <span class="keywordflow">for</span> (k = 2; k < 1000; k++) {</div><div class="line"><a name="l00058"></a><span class="lineno"> 58</span>  temp = ((k - 1 - a) * (lb - la) + (k + x) * lb) / k;</div><div class="line"><a name="l00059"></a><span class="lineno"> 59</span>  la = lb; lb = temp;</div><div class="line"><a name="l00060"></a><span class="lineno"> 60</span>  w *= (k - 1 - a) / k;</div><div class="line"><a name="l00061"></a><span class="lineno"> 61</span>  temp = w / (la * lb);</div><div class="line"><a name="l00062"></a><span class="lineno"> 62</span>  previous = <a class="code" href="../../d8/d90/nkf_8c.html#a5ea5ac7abf5cce39283e422add1067d5">result</a>; <a class="code" href="../../d8/d90/nkf_8c.html#a5ea5ac7abf5cce39283e422add1067d5">result</a> += temp;</div><div class="line"><a name="l00063"></a><span class="lineno"> 63</span>  <span class="keywordflow">if</span> (<a class="code" href="../../d8/d90/nkf_8c.html#a5ea5ac7abf5cce39283e422add1067d5">result</a> == previous) <span class="keywordflow">return</span> <a class="code" href="../../d8/d90/nkf_8c.html#a5ea5ac7abf5cce39283e422add1067d5">result</a>;</div><div class="line"><a name="l00064"></a><span class="lineno"> 64</span>  }</div><div class="line"><a name="l00065"></a><span class="lineno"> 65</span>  fprintf(stderr, <span class="stringliteral">"erf.c:%d:q_gamma() could not converge."</span>, __LINE__);</div><div class="line"><a name="l00066"></a><span class="lineno"> 66</span>  <span class="keywordflow">return</span> <a class="code" href="../../d8/d90/nkf_8c.html#a5ea5ac7abf5cce39283e422add1067d5">result</a>;</div><div class="line"><a name="l00067"></a><span class="lineno"> 67</span> }</div><div class="line"><a name="l00068"></a><span class="lineno"> 68</span> </div><div class="line"><a name="l00069"></a><span class="lineno"><a class="line" href="../../dd/dcf/erf_8c.html#a966a7f590647f053619298a7c884c554"> 69</a></span> <span class="preprocessor">#define LOG_PI_OVER_2 0.572364942924700087071713675675 </span><span class="comment">/* log_e(PI)/2 */</span><span class="preprocessor"></span></div><div class="line"><a name="l00070"></a><span class="lineno"> 70</span> </div><div class="line"><a name="l00071"></a><span class="lineno"><a class="line" href="../../dd/dcf/erf_8c.html#aa8b35382a71885eed509a5f87bf5e266"> 71</a></span> <span class="keywordtype">double</span> <a class="code" href="../../dd/dcf/erf_8c.html#aa8b35382a71885eed509a5f87bf5e266">erf</a>(<span class="keywordtype">double</span> x)</div><div class="line"><a name="l00072"></a><span class="lineno"> 72</span> {</div><div class="line"><a name="l00073"></a><span class="lineno"> 73</span>  <span class="keywordflow">if</span> (!<a class="code" href="../../d3/d90/missing_8h.html#a92974cf7271cdbcd972c211baf210f6f">finite</a>(x)) {</div><div class="line"><a name="l00074"></a><span class="lineno"> 74</span>  <span class="keywordflow">if</span> (<a class="code" href="../../dc/db1/win32_8h.html#a2e1baae9134e580910322362dc23290e">isnan</a>(x)) <span class="keywordflow">return</span> x; <span class="comment">/* erf(NaN) = NaN */</span></div><div class="line"><a name="l00075"></a><span class="lineno"> 75</span>  <span class="keywordflow">return</span> (x>0 ? 1.0 : -1.0); <span class="comment">/* erf(+-inf) = +-1.0 */</span></div><div class="line"><a name="l00076"></a><span class="lineno"> 76</span>  }</div><div class="line"><a name="l00077"></a><span class="lineno"> 77</span>  <span class="keywordflow">if</span> (x >= 0) <span class="keywordflow">return</span> <a class="code" href="../../dd/dcf/erf_8c.html#a2f5e1b011eda67f1fbd1b524b82dafb5">p_gamma</a>(0.5, x * x, <a class="code" href="../../dd/dcf/erf_8c.html#a966a7f590647f053619298a7c884c554">LOG_PI_OVER_2</a>);</div><div class="line"><a name="l00078"></a><span class="lineno"> 78</span>  <span class="keywordflow">else</span> <span class="keywordflow">return</span> - <a class="code" href="../../dd/dcf/erf_8c.html#a2f5e1b011eda67f1fbd1b524b82dafb5">p_gamma</a>(0.5, x * x, <a class="code" href="../../dd/dcf/erf_8c.html#a966a7f590647f053619298a7c884c554">LOG_PI_OVER_2</a>);</div><div class="line"><a name="l00079"></a><span class="lineno"> 79</span> }</div><div class="line"><a name="l00080"></a><span class="lineno"> 80</span> </div><div class="line"><a name="l00081"></a><span class="lineno"><a class="line" href="../../dd/dcf/erf_8c.html#a0b6c4cfae41124da64c9a2dcc09e6045"> 81</a></span> <span class="keywordtype">double</span> <a class="code" href="../../dd/dcf/erf_8c.html#a0b6c4cfae41124da64c9a2dcc09e6045">erfc</a>(<span class="keywordtype">double</span> x)</div><div class="line"><a name="l00082"></a><span class="lineno"> 82</span> {</div><div class="line"><a name="l00083"></a><span class="lineno"> 83</span>  <span class="keywordflow">if</span> (!<a class="code" href="../../d3/d90/missing_8h.html#a92974cf7271cdbcd972c211baf210f6f">finite</a>(x)) {</div><div class="line"><a name="l00084"></a><span class="lineno"> 84</span>  <span class="keywordflow">if</span> (<a class="code" href="../../dc/db1/win32_8h.html#a2e1baae9134e580910322362dc23290e">isnan</a>(x)) <span class="keywordflow">return</span> x; <span class="comment">/* erfc(NaN) = NaN */</span></div><div class="line"><a name="l00085"></a><span class="lineno"> 85</span>  <span class="keywordflow">return</span> (x>0 ? 0.0 : 2.0); <span class="comment">/* erfc(+-inf) = 0.0, 2.0 */</span></div><div class="line"><a name="l00086"></a><span class="lineno"> 86</span>  }</div><div class="line"><a name="l00087"></a><span class="lineno"> 87</span>  <span class="keywordflow">if</span> (x >= 0) <span class="keywordflow">return</span> <a class="code" href="../../dd/dcf/erf_8c.html#aa2417338852b65efc862b7cab59d216a">q_gamma</a>(0.5, x * x, <a class="code" href="../../dd/dcf/erf_8c.html#a966a7f590647f053619298a7c884c554">LOG_PI_OVER_2</a>);</div><div class="line"><a name="l00088"></a><span class="lineno"> 88</span>  <span class="keywordflow">else</span> <span class="keywordflow">return</span> 1 + <a class="code" href="../../dd/dcf/erf_8c.html#a2f5e1b011eda67f1fbd1b524b82dafb5">p_gamma</a>(0.5, x * x, <a class="code" href="../../dd/dcf/erf_8c.html#a966a7f590647f053619298a7c884c554">LOG_PI_OVER_2</a>);</div><div class="line"><a name="l00089"></a><span class="lineno"> 89</span> }</div><div class="ttc" id="missing_8h_html"><div class="ttname"><a href="../../d3/d90/missing_8h.html">missing.h</a></div></div> <div class="ttc" id="erf_8c_html_aa8b35382a71885eed509a5f87bf5e266"><div class="ttname"><a href="../../dd/dcf/erf_8c.html#aa8b35382a71885eed509a5f87bf5e266">erf</a></div><div class="ttdeci">double erf(double x)</div><div class="ttdef"><b>Definition:</b> <a href="../../dd/dcf/erf_8c_source.html#l00071">erf.c:71</a></div></div> <div class="ttc" id="missing_8h_html_a92974cf7271cdbcd972c211baf210f6f"><div class="ttname"><a href="../../d3/d90/missing_8h.html#a92974cf7271cdbcd972c211baf210f6f">finite</a></div><div class="ttdeci">RUBY_EXTERN int finite(double)</div><div class="ttdef"><b>Definition:</b> <a href="../../d4/d21/finite_8c_source.html#l00006">finite.c:6</a></div></div> <div class="ttc" id="erf_8c_html_a966a7f590647f053619298a7c884c554"><div class="ttname"><a href="../../dd/dcf/erf_8c.html#a966a7f590647f053619298a7c884c554">LOG_PI_OVER_2</a></div><div class="ttdeci">#define LOG_PI_OVER_2</div><div class="ttdef"><b>Definition:</b> <a href="../../dd/dcf/erf_8c_source.html#l00069">erf.c:69</a></div></div> <div class="ttc" id="erf_8c_html_aa2417338852b65efc862b7cab59d216a"><div class="ttname"><a href="../../dd/dcf/erf_8c.html#aa2417338852b65efc862b7cab59d216a">q_gamma</a></div><div class="ttdeci">static double q_gamma(double, double, double)</div><div class="ttdef"><b>Definition:</b> <a href="../../dd/dcf/erf_8c_source.html#l00048">erf.c:48</a></div></div> <div class="ttc" id="nkf_8c_html_a5ea5ac7abf5cce39283e422add1067d5"><div class="ttname"><a href="../../d8/d90/nkf_8c.html#a5ea5ac7abf5cce39283e422add1067d5">result</a></div><div class="ttdeci">static VALUE result</div><div class="ttdef"><b>Definition:</b> <a href="../../d8/d90/nkf_8c_source.html#l00040">nkf.c:40</a></div></div> <div class="ttc" id="win32_8h_html_a2e1baae9134e580910322362dc23290e"><div class="ttname"><a href="../../dc/db1/win32_8h.html#a2e1baae9134e580910322362dc23290e">isnan</a></div><div class="ttdeci">#define isnan(x)</div><div class="ttdef"><b>Definition:</b> <a href="../../dc/db1/win32_8h_source.html#l00334">win32.h:334</a></div></div> <div class="ttc" id="erf_8c_html_a2f5e1b011eda67f1fbd1b524b82dafb5"><div class="ttname"><a href="../../dd/dcf/erf_8c.html#a2f5e1b011eda67f1fbd1b524b82dafb5">p_gamma</a></div><div class="ttdeci">static double p_gamma(double a, double x, double loggamma_a)</div><div class="ttdef"><b>Definition:</b> <a href="../../dd/dcf/erf_8c_source.html#l00029">erf.c:29</a></div></div> <div class="ttc" id="erf_8c_html_a0b6c4cfae41124da64c9a2dcc09e6045"><div class="ttname"><a href="../../dd/dcf/erf_8c.html#a0b6c4cfae41124da64c9a2dcc09e6045">erfc</a></div><div class="ttdeci">double erfc(double x)</div><div class="ttdef"><b>Definition:</b> <a href="../../dd/dcf/erf_8c_source.html#l00081">erf.c:81</a></div></div> </div><!-- fragment --></div><!-- contents --> <!-- start footer part --> <hr class="footer"/><address class="footer"><small> Generated by  <a href="http://www.doxygen.org/index.html"> <img class="footer" src="../../doxygen.png" alt="doxygen"/> </a> 1.8.14 </small></address> </body> </html>
💾 Save Changes
❌ Cancel