<?xml version='1.0' encoding='UTF-8' standalone='no'?>
<doxygen xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xsi:noNamespaceSchemaLocation="compound.xsd" version="1.5.6">
  <compounddef id="R__rgamma_8c" kind="file">
    <compoundname>R_rgamma.c</compoundname>
    <includes local="yes">nmath.h</includes>
    <incdepgraph>
      <node id="318">
        <label>work/git/mixpp/bdm/math/R_rgamma.c</label>
        <link refid="R_rgamma.c"/>
        <childnode refid="319" relation="include">
        </childnode>
      </node>
      <node id="319">
        <label>nmath.h</label>
      </node>
    </incdepgraph>
      <sectiondef kind="define">
      <memberdef kind="define" id="R__rgamma_8c_188f47a79aa968bbb4c64a245c2c14437" prot="public" static="no">
        <name>repeat</name>
        <initializer>for(;;)</initializer>
        <briefdescription>
        </briefdescription>
        <detaileddescription>
        </detaileddescription>
        <inbodydescription>
        </inbodydescription>
        <location file="/home/smidl/work/git/mixpp/bdm/math/R_rgamma.c" line="52" bodyfile="/home/smidl/work/git/mixpp/bdm/math/R_rgamma.c" bodystart="52" bodyend="-1"/>
      </memberdef>
      </sectiondef>
      <sectiondef kind="func">
      <memberdef kind="function" id="R__rgamma_8c_1a593053bf473371b1f7368b592c68b95" prot="public" static="no" const="no" explicit="no" inline="no" virt="non-virtual">
        <type>double</type>
        <definition>double rgamma</definition>
        <argsstring>(double a, double scale)</argsstring>
        <name>rgamma</name>
        <param>
          <type>double</type>
          <declname>a</declname>
        </param>
        <param>
          <type>double</type>
          <declname>scale</declname>
        </param>
        <briefdescription>
        </briefdescription>
        <detaileddescription>
        </detaileddescription>
        <inbodydescription>
        </inbodydescription>
        <location file="/home/smidl/work/git/mixpp/bdm/math/R_rgamma.c" line="55" bodyfile="/home/smidl/work/git/mixpp/bdm/math/R_rgamma.c" bodystart="54" bodyend="211"/>
      </memberdef>
      </sectiondef>
    <briefdescription>
    </briefdescription>
    <detaileddescription>
    </detaileddescription>
    <programlisting>
<codeline lineno="1"><highlight class="comment">/*</highlight></codeline>
<codeline lineno="2"><highlight class="comment"><sp/>*<sp/><sp/>Mathlib<sp/>:<sp/>A<sp/>C<sp/>Library<sp/>of<sp/>Special<sp/>Functions</highlight></codeline>
<codeline lineno="3"><highlight class="comment"><sp/>*<sp/><sp/>Copyright<sp/>(C)<sp/>1998<sp/>Ross<sp/>Ihaka</highlight></codeline>
<codeline lineno="4"><highlight class="comment"><sp/>*<sp/><sp/>Copyright<sp/>(C)<sp/>2000-2005<sp/>The<sp/>R<sp/>Development<sp/>Core<sp/>Team</highlight></codeline>
<codeline lineno="5"><highlight class="comment"><sp/>*</highlight></codeline>
<codeline lineno="6"><highlight class="comment"><sp/>*<sp/><sp/>This<sp/>program<sp/>is<sp/>free<sp/>software;<sp/>you<sp/>can<sp/>redistribute<sp/>it<sp/>and/or<sp/>modify</highlight></codeline>
<codeline lineno="7"><highlight class="comment"><sp/>*<sp/><sp/>it<sp/>under<sp/>the<sp/>terms<sp/>of<sp/>the<sp/>GNU<sp/>General<sp/>Public<sp/>License<sp/>as<sp/>published<sp/>by</highlight></codeline>
<codeline lineno="8"><highlight class="comment"><sp/>*<sp/><sp/>the<sp/>Free<sp/>Software<sp/>Foundation;<sp/>either<sp/>version<sp/>2<sp/>of<sp/>the<sp/>License,<sp/>or</highlight></codeline>
<codeline lineno="9"><highlight class="comment"><sp/>*<sp/><sp/>(at<sp/>your<sp/>option)<sp/>any<sp/>later<sp/>version.</highlight></codeline>
<codeline lineno="10"><highlight class="comment"><sp/>*</highlight></codeline>
<codeline lineno="11"><highlight class="comment"><sp/>*<sp/><sp/>This<sp/>program<sp/>is<sp/>distributed<sp/>in<sp/>the<sp/>hope<sp/>that<sp/>it<sp/>will<sp/>be<sp/>useful,</highlight></codeline>
<codeline lineno="12"><highlight class="comment"><sp/>*<sp/><sp/>but<sp/>WITHOUT<sp/>ANY<sp/>WARRANTY;<sp/>without<sp/>even<sp/>the<sp/>implied<sp/>warranty<sp/>of</highlight></codeline>
<codeline lineno="13"><highlight class="comment"><sp/>*<sp/><sp/>MERCHANTABILITY<sp/>or<sp/>FITNESS<sp/>FOR<sp/>A<sp/>PARTICULAR<sp/>PURPOSE.<sp/><sp/>See<sp/>the</highlight></codeline>
<codeline lineno="14"><highlight class="comment"><sp/>*<sp/><sp/>GNU<sp/>General<sp/>Public<sp/>License<sp/>for<sp/>more<sp/>details.</highlight></codeline>
<codeline lineno="15"><highlight class="comment"><sp/>*</highlight></codeline>
<codeline lineno="16"><highlight class="comment"><sp/>*<sp/><sp/>You<sp/>should<sp/>have<sp/>received<sp/>a<sp/>copy<sp/>of<sp/>the<sp/>GNU<sp/>General<sp/>Public<sp/>License</highlight></codeline>
<codeline lineno="17"><highlight class="comment"><sp/>*<sp/><sp/>along<sp/>with<sp/>this<sp/>program;<sp/>if<sp/>not,<sp/>a<sp/>copy<sp/>is<sp/>available<sp/>at</highlight></codeline>
<codeline lineno="18"><highlight class="comment"><sp/>*<sp/><sp/>http://www.r-project.org/Licenses/</highlight></codeline>
<codeline lineno="19"><highlight class="comment"><sp/>*</highlight></codeline>
<codeline lineno="20"><highlight class="comment"><sp/>*<sp/><sp/>SYNOPSIS</highlight></codeline>
<codeline lineno="21"><highlight class="comment"><sp/>*</highlight></codeline>
<codeline lineno="22"><highlight class="comment"><sp/>*<sp/><sp/><sp/><sp/>#include<sp/>&lt;Rmath.h&gt;</highlight></codeline>
<codeline lineno="23"><highlight class="comment"><sp/>*<sp/><sp/><sp/><sp/>double<sp/>rgamma(double<sp/>a,<sp/>double<sp/>scale);</highlight></codeline>
<codeline lineno="24"><highlight class="comment"><sp/>*</highlight></codeline>
<codeline lineno="25"><highlight class="comment"><sp/>*<sp/><sp/>DESCRIPTION</highlight></codeline>
<codeline lineno="26"><highlight class="comment"><sp/>*</highlight></codeline>
<codeline lineno="27"><highlight class="comment"><sp/>*<sp/><sp/><sp/><sp/>Random<sp/>variates<sp/>from<sp/>the<sp/>gamma<sp/>distribution.</highlight></codeline>
<codeline lineno="28"><highlight class="comment"><sp/>*</highlight></codeline>
<codeline lineno="29"><highlight class="comment"><sp/>*<sp/><sp/>REFERENCES</highlight></codeline>
<codeline lineno="30"><highlight class="comment"><sp/>*</highlight></codeline>
<codeline lineno="31"><highlight class="comment"><sp/>*<sp/><sp/><sp/><sp/>[1]<sp/>Shape<sp/>parameter<sp/>a<sp/>&gt;=<sp/>1.<sp/><sp/>Algorithm<sp/>GD<sp/>in:</highlight></codeline>
<codeline lineno="32"><highlight class="comment"><sp/>*</highlight></codeline>
<codeline lineno="33"><highlight class="comment"><sp/>*<sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/>Ahrens,<sp/>J.H.<sp/>and<sp/>Dieter,<sp/>U.<sp/>(1982).</highlight></codeline>
<codeline lineno="34"><highlight class="comment"><sp/>*<sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/>Generating<sp/>gamma<sp/>variates<sp/>by<sp/>a<sp/>modified</highlight></codeline>
<codeline lineno="35"><highlight class="comment"><sp/>*<sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/>rejection<sp/>technique.</highlight></codeline>
<codeline lineno="36"><highlight class="comment"><sp/>*<sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/>Comm.<sp/>ACM,<sp/>25,<sp/>47-54.</highlight></codeline>
<codeline lineno="37"><highlight class="comment"><sp/>*</highlight></codeline>
<codeline lineno="38"><highlight class="comment"><sp/>*</highlight></codeline>
<codeline lineno="39"><highlight class="comment"><sp/>*<sp/><sp/><sp/><sp/>[2]<sp/>Shape<sp/>parameter<sp/>0<sp/>&lt;<sp/>a<sp/>&lt;<sp/>1.<sp/>Algorithm<sp/>GS<sp/>in:</highlight></codeline>
<codeline lineno="40"><highlight class="comment"><sp/>*</highlight></codeline>
<codeline lineno="41"><highlight class="comment"><sp/>*<sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/>Ahrens,<sp/>J.H.<sp/>and<sp/>Dieter,<sp/>U.<sp/>(1974).</highlight></codeline>
<codeline lineno="42"><highlight class="comment"><sp/>*<sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/>Computer<sp/>methods<sp/>for<sp/>sampling<sp/>from<sp/>gamma,<sp/>beta,</highlight></codeline>
<codeline lineno="43"><highlight class="comment"><sp/>*<sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/>poisson<sp/>and<sp/>binomial<sp/>distributions.</highlight></codeline>
<codeline lineno="44"><highlight class="comment"><sp/>*<sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/>Computing,<sp/>12,<sp/>223-246.</highlight></codeline>
<codeline lineno="45"><highlight class="comment"><sp/>*</highlight></codeline>
<codeline lineno="46"><highlight class="comment"><sp/>*<sp/><sp/><sp/><sp/>Input:<sp/>a<sp/>=<sp/>parameter<sp/>(mean)<sp/>of<sp/>the<sp/>standard<sp/>gamma<sp/>distribution.</highlight></codeline>
<codeline lineno="47"><highlight class="comment"><sp/>*<sp/><sp/><sp/><sp/>Output:<sp/>a<sp/>variate<sp/>from<sp/>the<sp/>gamma(a)-distribution</highlight></codeline>
<codeline lineno="48"><highlight class="comment"><sp/>*/</highlight><highlight class="normal"></highlight></codeline>
<codeline lineno="49"><highlight class="normal"></highlight></codeline>
<codeline lineno="50"><highlight class="normal"></highlight><highlight class="preprocessor">#include<sp/>&quot;nmath.h&quot;</highlight><highlight class="normal"></highlight></codeline>
<codeline lineno="51"><highlight class="normal"></highlight></codeline>
<codeline lineno="52"><highlight class="normal"></highlight><highlight class="preprocessor">#define<sp/>repeat<sp/>for(;;)</highlight></codeline>
<codeline lineno="53"><highlight class="preprocessor"></highlight><highlight class="normal"></highlight></codeline>
<codeline lineno="54"><highlight class="normal"></highlight><highlight class="keywordtype">double</highlight><highlight class="normal"><sp/>rgamma(</highlight><highlight class="keywordtype">double</highlight><highlight class="normal"><sp/>a,<sp/></highlight><highlight class="keywordtype">double</highlight><highlight class="normal"><sp/>scale)</highlight></codeline>
<codeline lineno="55"><highlight class="normal">{</highlight></codeline>
<codeline lineno="56"><highlight class="normal"></highlight><highlight class="comment">/*<sp/>Constants<sp/>:<sp/>*/</highlight><highlight class="normal"></highlight></codeline>
<codeline lineno="57"><highlight class="normal"><sp/><sp/><sp/><sp/></highlight><highlight class="keyword">const</highlight><highlight class="normal"><sp/></highlight><highlight class="keyword">static</highlight><highlight class="normal"><sp/></highlight><highlight class="keywordtype">double</highlight><highlight class="normal"><sp/>sqrt32<sp/>=<sp/>5.656854;</highlight></codeline>
<codeline lineno="58"><highlight class="normal"><sp/><sp/><sp/><sp/></highlight><highlight class="keyword">const</highlight><highlight class="normal"><sp/></highlight><highlight class="keyword">static</highlight><highlight class="normal"><sp/></highlight><highlight class="keywordtype">double</highlight><highlight class="normal"><sp/>exp_m1<sp/>=<sp/>0.36787944117144232159;</highlight><highlight class="comment">/*<sp/>exp(-1)<sp/>=<sp/>1/e<sp/>*/</highlight><highlight class="normal"></highlight></codeline>
<codeline lineno="59"><highlight class="normal"></highlight></codeline>
<codeline lineno="60"><highlight class="normal"><sp/><sp/><sp/><sp/></highlight><highlight class="comment">/*<sp/>Coefficients<sp/>q[k]<sp/>-<sp/>for<sp/>q0<sp/>=<sp/>sum(q[k]*a^(-k))</highlight></codeline>
<codeline lineno="61"><highlight class="comment"><sp/><sp/><sp/><sp/><sp/>*<sp/>Coefficients<sp/>a[k]<sp/>-<sp/>for<sp/>q<sp/>=<sp/>q0+(t*t/2)*sum(a[k]*v^k)</highlight></codeline>
<codeline lineno="62"><highlight class="comment"><sp/><sp/><sp/><sp/><sp/>*<sp/>Coefficients<sp/>e[k]<sp/>-<sp/>for<sp/>exp(q)-1<sp/>=<sp/>sum(e[k]*q^k)</highlight></codeline>
<codeline lineno="63"><highlight class="comment"><sp/><sp/><sp/><sp/><sp/>*/</highlight><highlight class="normal"></highlight></codeline>
<codeline lineno="64"><highlight class="normal"><sp/><sp/><sp/><sp/></highlight><highlight class="keyword">const</highlight><highlight class="normal"><sp/></highlight><highlight class="keyword">static</highlight><highlight class="normal"><sp/></highlight><highlight class="keywordtype">double</highlight><highlight class="normal"><sp/>q1<sp/>=<sp/>0.04166669;</highlight></codeline>
<codeline lineno="65"><highlight class="normal"><sp/><sp/><sp/><sp/></highlight><highlight class="keyword">const</highlight><highlight class="normal"><sp/></highlight><highlight class="keyword">static</highlight><highlight class="normal"><sp/></highlight><highlight class="keywordtype">double</highlight><highlight class="normal"><sp/>q2<sp/>=<sp/>0.02083148;</highlight></codeline>
<codeline lineno="66"><highlight class="normal"><sp/><sp/><sp/><sp/></highlight><highlight class="keyword">const</highlight><highlight class="normal"><sp/></highlight><highlight class="keyword">static</highlight><highlight class="normal"><sp/></highlight><highlight class="keywordtype">double</highlight><highlight class="normal"><sp/>q3<sp/>=<sp/>0.00801191;</highlight></codeline>
<codeline lineno="67"><highlight class="normal"><sp/><sp/><sp/><sp/></highlight><highlight class="keyword">const</highlight><highlight class="normal"><sp/></highlight><highlight class="keyword">static</highlight><highlight class="normal"><sp/></highlight><highlight class="keywordtype">double</highlight><highlight class="normal"><sp/>q4<sp/>=<sp/>0.00144121;</highlight></codeline>
<codeline lineno="68"><highlight class="normal"><sp/><sp/><sp/><sp/></highlight><highlight class="keyword">const</highlight><highlight class="normal"><sp/></highlight><highlight class="keyword">static</highlight><highlight class="normal"><sp/></highlight><highlight class="keywordtype">double</highlight><highlight class="normal"><sp/>q5<sp/>=<sp/>-7.388e-5;</highlight></codeline>
<codeline lineno="69"><highlight class="normal"><sp/><sp/><sp/><sp/></highlight><highlight class="keyword">const</highlight><highlight class="normal"><sp/></highlight><highlight class="keyword">static</highlight><highlight class="normal"><sp/></highlight><highlight class="keywordtype">double</highlight><highlight class="normal"><sp/>q6<sp/>=<sp/>2.4511e-4;</highlight></codeline>
<codeline lineno="70"><highlight class="normal"><sp/><sp/><sp/><sp/></highlight><highlight class="keyword">const</highlight><highlight class="normal"><sp/></highlight><highlight class="keyword">static</highlight><highlight class="normal"><sp/></highlight><highlight class="keywordtype">double</highlight><highlight class="normal"><sp/>q7<sp/>=<sp/>2.424e-4;</highlight></codeline>
<codeline lineno="71"><highlight class="normal"></highlight></codeline>
<codeline lineno="72"><highlight class="normal"><sp/><sp/><sp/><sp/></highlight><highlight class="keyword">const</highlight><highlight class="normal"><sp/></highlight><highlight class="keyword">static</highlight><highlight class="normal"><sp/></highlight><highlight class="keywordtype">double</highlight><highlight class="normal"><sp/>a1<sp/>=<sp/>0.3333333;</highlight></codeline>
<codeline lineno="73"><highlight class="normal"><sp/><sp/><sp/><sp/></highlight><highlight class="keyword">const</highlight><highlight class="normal"><sp/></highlight><highlight class="keyword">static</highlight><highlight class="normal"><sp/></highlight><highlight class="keywordtype">double</highlight><highlight class="normal"><sp/>a2<sp/>=<sp/>-0.250003;</highlight></codeline>
<codeline lineno="74"><highlight class="normal"><sp/><sp/><sp/><sp/></highlight><highlight class="keyword">const</highlight><highlight class="normal"><sp/></highlight><highlight class="keyword">static</highlight><highlight class="normal"><sp/></highlight><highlight class="keywordtype">double</highlight><highlight class="normal"><sp/>a3<sp/>=<sp/>0.2000062;</highlight></codeline>
<codeline lineno="75"><highlight class="normal"><sp/><sp/><sp/><sp/></highlight><highlight class="keyword">const</highlight><highlight class="normal"><sp/></highlight><highlight class="keyword">static</highlight><highlight class="normal"><sp/></highlight><highlight class="keywordtype">double</highlight><highlight class="normal"><sp/>a4<sp/>=<sp/>-0.1662921;</highlight></codeline>
<codeline lineno="76"><highlight class="normal"><sp/><sp/><sp/><sp/></highlight><highlight class="keyword">const</highlight><highlight class="normal"><sp/></highlight><highlight class="keyword">static</highlight><highlight class="normal"><sp/></highlight><highlight class="keywordtype">double</highlight><highlight class="normal"><sp/>a5<sp/>=<sp/>0.1423657;</highlight></codeline>
<codeline lineno="77"><highlight class="normal"><sp/><sp/><sp/><sp/></highlight><highlight class="keyword">const</highlight><highlight class="normal"><sp/></highlight><highlight class="keyword">static</highlight><highlight class="normal"><sp/></highlight><highlight class="keywordtype">double</highlight><highlight class="normal"><sp/>a6<sp/>=<sp/>-0.1367177;</highlight></codeline>
<codeline lineno="78"><highlight class="normal"><sp/><sp/><sp/><sp/></highlight><highlight class="keyword">const</highlight><highlight class="normal"><sp/></highlight><highlight class="keyword">static</highlight><highlight class="normal"><sp/></highlight><highlight class="keywordtype">double</highlight><highlight class="normal"><sp/>a7<sp/>=<sp/>0.1233795;</highlight></codeline>
<codeline lineno="79"><highlight class="normal"></highlight></codeline>
<codeline lineno="80"><highlight class="normal"><sp/><sp/><sp/><sp/></highlight><highlight class="comment">/*<sp/>State<sp/>variables<sp/>[FIXME<sp/>for<sp/>threading!]<sp/>:*/</highlight><highlight class="normal"></highlight></codeline>
<codeline lineno="81"><highlight class="normal"><sp/><sp/><sp/><sp/></highlight><highlight class="keyword">static</highlight><highlight class="normal"><sp/></highlight><highlight class="keywordtype">double</highlight><highlight class="normal"><sp/>aa<sp/>=<sp/>0.;</highlight></codeline>
<codeline lineno="82"><highlight class="normal"><sp/><sp/><sp/><sp/></highlight><highlight class="keyword">static</highlight><highlight class="normal"><sp/></highlight><highlight class="keywordtype">double</highlight><highlight class="normal"><sp/>aaa<sp/>=<sp/>0.;</highlight></codeline>
<codeline lineno="83"><highlight class="normal"><sp/><sp/><sp/><sp/></highlight><highlight class="keyword">static</highlight><highlight class="normal"><sp/></highlight><highlight class="keywordtype">double</highlight><highlight class="normal"><sp/>s,<sp/>s2,<sp/>d;<sp/><sp/><sp/><sp/></highlight><highlight class="comment">/*<sp/>no.<sp/>1<sp/>(step<sp/>1)<sp/>*/</highlight><highlight class="normal"></highlight></codeline>
<codeline lineno="84"><highlight class="normal"><sp/><sp/><sp/><sp/></highlight><highlight class="keyword">static</highlight><highlight class="normal"><sp/></highlight><highlight class="keywordtype">double</highlight><highlight class="normal"><sp/>q0,<sp/>b,<sp/>si,<sp/>c;</highlight><highlight class="comment">/*<sp/>no.<sp/>2<sp/>(step<sp/>4)<sp/>*/</highlight><highlight class="normal"></highlight></codeline>
<codeline lineno="85"><highlight class="normal"></highlight></codeline>
<codeline lineno="86"><highlight class="normal"><sp/><sp/><sp/><sp/></highlight><highlight class="keywordtype">double</highlight><highlight class="normal"><sp/>e,<sp/>p,<sp/>q,<sp/>r,<sp/>t,<sp/>u,<sp/>v,<sp/>w,<sp/>x,<sp/>ret_val;</highlight></codeline>
<codeline lineno="87"><highlight class="normal"></highlight></codeline>
<codeline lineno="88"><highlight class="normal"><sp/><sp/><sp/><sp/></highlight><highlight class="keywordflow">if</highlight><highlight class="normal"><sp/>(!R_FINITE(a)<sp/>||<sp/>!R_FINITE(scale)<sp/>||<sp/>a<sp/>&lt;<sp/>0.0<sp/>||<sp/>scale<sp/>&lt;=<sp/>0.0)</highlight></codeline>
<codeline lineno="89"><highlight class="normal"><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/>ML_ERR_return_NAN;</highlight></codeline>
<codeline lineno="90"><highlight class="normal"></highlight></codeline>
<codeline lineno="91"><highlight class="normal"><sp/><sp/><sp/><sp/></highlight><highlight class="keywordflow">if</highlight><highlight class="normal"><sp/>(a<sp/>&lt;<sp/>1.)<sp/>{<sp/></highlight><highlight class="comment">/*<sp/>GS<sp/>algorithm<sp/>for<sp/>parameters<sp/>a<sp/>&lt;<sp/>1<sp/>*/</highlight><highlight class="normal"></highlight></codeline>
<codeline lineno="92"><highlight class="normal"><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/></highlight><highlight class="keywordflow">if</highlight><highlight class="normal">(a<sp/>==<sp/>0)</highlight></codeline>
<codeline lineno="93"><highlight class="normal"><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/></highlight><highlight class="keywordflow">return</highlight><highlight class="normal"><sp/>0.;</highlight></codeline>
<codeline lineno="94"><highlight class="normal"><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/>e<sp/>=<sp/>1.0<sp/>+<sp/>exp_m1<sp/>*<sp/>a;</highlight></codeline>
<codeline lineno="95"><highlight class="normal"><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/>repeat<sp/>{</highlight></codeline>
<codeline lineno="96"><highlight class="normal"><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/>p<sp/>=<sp/>e<sp/>*<sp/>unif_rand();</highlight></codeline>
<codeline lineno="97"><highlight class="normal"><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/></highlight><highlight class="keywordflow">if</highlight><highlight class="normal"><sp/>(p<sp/>&gt;=<sp/>1.0)<sp/>{</highlight></codeline>
<codeline lineno="98"><highlight class="normal"><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/>x<sp/>=<sp/>-log((e<sp/>-<sp/>p)<sp/>/<sp/>a);</highlight></codeline>
<codeline lineno="99"><highlight class="normal"><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/></highlight><highlight class="keywordflow">if</highlight><highlight class="normal"><sp/>(exp_rand()<sp/>&gt;=<sp/>(1.0<sp/>-<sp/>a)<sp/>*<sp/>log(x))</highlight></codeline>
<codeline lineno="100"><highlight class="normal"><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/></highlight><highlight class="keywordflow">break</highlight><highlight class="normal">;</highlight></codeline>
<codeline lineno="101"><highlight class="normal"><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/>}<sp/></highlight><highlight class="keywordflow">else</highlight><highlight class="normal"><sp/>{</highlight></codeline>
<codeline lineno="102"><highlight class="normal"><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/>x<sp/>=<sp/>exp(log(p)<sp/>/<sp/>a);</highlight></codeline>
<codeline lineno="103"><highlight class="normal"><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/></highlight><highlight class="keywordflow">if</highlight><highlight class="normal"><sp/>(exp_rand()<sp/>&gt;=<sp/>x)</highlight></codeline>
<codeline lineno="104"><highlight class="normal"><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/></highlight><highlight class="keywordflow">break</highlight><highlight class="normal">;</highlight></codeline>
<codeline lineno="105"><highlight class="normal"><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/>}</highlight></codeline>
<codeline lineno="106"><highlight class="normal"><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/>}</highlight></codeline>
<codeline lineno="107"><highlight class="normal"><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/></highlight><highlight class="keywordflow">return</highlight><highlight class="normal"><sp/>scale<sp/>*<sp/>x;</highlight></codeline>
<codeline lineno="108"><highlight class="normal"><sp/><sp/><sp/><sp/>}</highlight></codeline>
<codeline lineno="109"><highlight class="normal"></highlight></codeline>
<codeline lineno="110"><highlight class="normal"><sp/><sp/><sp/><sp/></highlight><highlight class="comment">/*<sp/>---<sp/>a<sp/>&gt;=<sp/>1<sp/>:<sp/>GD<sp/>algorithm<sp/>---<sp/>*/</highlight><highlight class="normal"></highlight></codeline>
<codeline lineno="111"><highlight class="normal"></highlight></codeline>
<codeline lineno="112"><highlight class="normal"><sp/><sp/><sp/><sp/></highlight><highlight class="comment">/*<sp/>Step<sp/>1:<sp/>Recalculations<sp/>of<sp/>s2,<sp/>s,<sp/>d<sp/>if<sp/>a<sp/>has<sp/>changed<sp/>*/</highlight><highlight class="normal"></highlight></codeline>
<codeline lineno="113"><highlight class="normal"><sp/><sp/><sp/><sp/></highlight><highlight class="keywordflow">if</highlight><highlight class="normal"><sp/>(a<sp/>!=<sp/>aa)<sp/>{</highlight></codeline>
<codeline lineno="114"><highlight class="normal"><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/>aa<sp/>=<sp/>a;</highlight></codeline>
<codeline lineno="115"><highlight class="normal"><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/>s2<sp/>=<sp/>a<sp/>-<sp/>0.5;</highlight></codeline>
<codeline lineno="116"><highlight class="normal"><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/>s<sp/>=<sp/>sqrt(s2);</highlight></codeline>
<codeline lineno="117"><highlight class="normal"><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/>d<sp/>=<sp/>sqrt32<sp/>-<sp/>s<sp/>*<sp/>12.0;</highlight></codeline>
<codeline lineno="118"><highlight class="normal"><sp/><sp/><sp/><sp/>}</highlight></codeline>
<codeline lineno="119"><highlight class="normal"><sp/><sp/><sp/><sp/></highlight><highlight class="comment">/*<sp/>Step<sp/>2:<sp/>t<sp/>=<sp/>standard<sp/>normal<sp/>deviate,</highlight></codeline>
<codeline lineno="120"><highlight class="comment"><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/>x<sp/>=<sp/>(s,1/2)<sp/>-normal<sp/>deviate.<sp/>*/</highlight><highlight class="normal"></highlight></codeline>
<codeline lineno="121"><highlight class="normal"></highlight></codeline>
<codeline lineno="122"><highlight class="normal"><sp/><sp/><sp/><sp/></highlight><highlight class="comment">/*<sp/>immediate<sp/>acceptance<sp/>(i)<sp/>*/</highlight><highlight class="normal"></highlight></codeline>
<codeline lineno="123"><highlight class="normal"><sp/><sp/><sp/><sp/>t<sp/>=<sp/>norm_rand();</highlight></codeline>
<codeline lineno="124"><highlight class="normal"><sp/><sp/><sp/><sp/>x<sp/>=<sp/>s<sp/>+<sp/>0.5<sp/>*<sp/>t;</highlight></codeline>
<codeline lineno="125"><highlight class="normal"><sp/><sp/><sp/><sp/>ret_val<sp/>=<sp/>x<sp/>*<sp/>x;</highlight></codeline>
<codeline lineno="126"><highlight class="normal"><sp/><sp/><sp/><sp/></highlight><highlight class="keywordflow">if</highlight><highlight class="normal"><sp/>(t<sp/>&gt;=<sp/>0.0)</highlight></codeline>
<codeline lineno="127"><highlight class="normal"><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/></highlight><highlight class="keywordflow">return</highlight><highlight class="normal"><sp/>scale<sp/>*<sp/>ret_val;</highlight></codeline>
<codeline lineno="128"><highlight class="normal"></highlight></codeline>
<codeline lineno="129"><highlight class="normal"><sp/><sp/><sp/><sp/></highlight><highlight class="comment">/*<sp/>Step<sp/>3:<sp/>u<sp/>=<sp/>0,1<sp/>-<sp/>uniform<sp/>sample.<sp/>squeeze<sp/>acceptance<sp/>(s)<sp/>*/</highlight><highlight class="normal"></highlight></codeline>
<codeline lineno="130"><highlight class="normal"><sp/><sp/><sp/><sp/>u<sp/>=<sp/>unif_rand();</highlight></codeline>
<codeline lineno="131"><highlight class="normal"><sp/><sp/><sp/><sp/></highlight><highlight class="keywordflow">if</highlight><highlight class="normal"><sp/>(d<sp/>*<sp/>u<sp/>&lt;=<sp/>t<sp/>*<sp/>t<sp/>*<sp/>t)</highlight></codeline>
<codeline lineno="132"><highlight class="normal"><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/></highlight><highlight class="keywordflow">return</highlight><highlight class="normal"><sp/>scale<sp/>*<sp/>ret_val;</highlight></codeline>
<codeline lineno="133"><highlight class="normal"></highlight></codeline>
<codeline lineno="134"><highlight class="normal"><sp/><sp/><sp/><sp/></highlight><highlight class="comment">/*<sp/>Step<sp/>4:<sp/>recalculations<sp/>of<sp/>q0,<sp/>b,<sp/>si,<sp/>c<sp/>if<sp/>necessary<sp/>*/</highlight><highlight class="normal"></highlight></codeline>
<codeline lineno="135"><highlight class="normal"></highlight></codeline>
<codeline lineno="136"><highlight class="normal"><sp/><sp/><sp/><sp/></highlight><highlight class="keywordflow">if</highlight><highlight class="normal"><sp/>(a<sp/>!=<sp/>aaa)<sp/>{</highlight></codeline>
<codeline lineno="137"><highlight class="normal"><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/>aaa<sp/>=<sp/>a;</highlight></codeline>
<codeline lineno="138"><highlight class="normal"><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/>r<sp/>=<sp/>1.0<sp/>/<sp/>a;</highlight></codeline>
<codeline lineno="139"><highlight class="normal"><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/>q0<sp/>=<sp/>((((((q7<sp/>*<sp/>r<sp/>+<sp/>q6)<sp/>*<sp/>r<sp/>+<sp/>q5)<sp/>*<sp/>r<sp/>+<sp/>q4)<sp/>*<sp/>r<sp/>+<sp/>q3)<sp/>*<sp/>r</highlight></codeline>
<codeline lineno="140"><highlight class="normal"><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/>+<sp/>q2)<sp/>*<sp/>r<sp/>+<sp/>q1)<sp/>*<sp/>r;</highlight></codeline>
<codeline lineno="141"><highlight class="normal"></highlight></codeline>
<codeline lineno="142"><highlight class="normal"><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/></highlight><highlight class="comment">/*<sp/>Approximation<sp/>depending<sp/>on<sp/>size<sp/>of<sp/>parameter<sp/>a<sp/>*/</highlight><highlight class="normal"></highlight></codeline>
<codeline lineno="143"><highlight class="normal"><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/></highlight><highlight class="comment">/*<sp/>The<sp/>constants<sp/>in<sp/>the<sp/>expressions<sp/>for<sp/>b,<sp/>si<sp/>and<sp/>c<sp/>*/</highlight><highlight class="normal"></highlight></codeline>
<codeline lineno="144"><highlight class="normal"><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/></highlight><highlight class="comment">/*<sp/>were<sp/>established<sp/>by<sp/>numerical<sp/>experiments<sp/>*/</highlight><highlight class="normal"></highlight></codeline>
<codeline lineno="145"><highlight class="normal"></highlight></codeline>
<codeline lineno="146"><highlight class="normal"><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/></highlight><highlight class="keywordflow">if</highlight><highlight class="normal"><sp/>(a<sp/>&lt;=<sp/>3.686)<sp/>{</highlight></codeline>
<codeline lineno="147"><highlight class="normal"><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/>b<sp/>=<sp/>0.463<sp/>+<sp/>s<sp/>+<sp/>0.178<sp/>*<sp/>s2;</highlight></codeline>
<codeline lineno="148"><highlight class="normal"><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/>si<sp/>=<sp/>1.235;</highlight></codeline>
<codeline lineno="149"><highlight class="normal"><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/>c<sp/>=<sp/>0.195<sp/>/<sp/>s<sp/>-<sp/>0.079<sp/>+<sp/>0.16<sp/>*<sp/>s;</highlight></codeline>
<codeline lineno="150"><highlight class="normal"><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/>}<sp/></highlight><highlight class="keywordflow">else</highlight><highlight class="normal"><sp/></highlight><highlight class="keywordflow">if</highlight><highlight class="normal"><sp/>(a<sp/>&lt;=<sp/>13.022)<sp/>{</highlight></codeline>
<codeline lineno="151"><highlight class="normal"><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/>b<sp/>=<sp/>1.654<sp/>+<sp/>0.0076<sp/>*<sp/>s2;</highlight></codeline>
<codeline lineno="152"><highlight class="normal"><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/>si<sp/>=<sp/>1.68<sp/>/<sp/>s<sp/>+<sp/>0.275;</highlight></codeline>
<codeline lineno="153"><highlight class="normal"><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/>c<sp/>=<sp/>0.062<sp/>/<sp/>s<sp/>+<sp/>0.024;</highlight></codeline>
<codeline lineno="154"><highlight class="normal"><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/>}<sp/></highlight><highlight class="keywordflow">else</highlight><highlight class="normal"><sp/>{</highlight></codeline>
<codeline lineno="155"><highlight class="normal"><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/>b<sp/>=<sp/>1.77;</highlight></codeline>
<codeline lineno="156"><highlight class="normal"><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/>si<sp/>=<sp/>0.75;</highlight></codeline>
<codeline lineno="157"><highlight class="normal"><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/>c<sp/>=<sp/>0.1515<sp/>/<sp/>s;</highlight></codeline>
<codeline lineno="158"><highlight class="normal"><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/>}</highlight></codeline>
<codeline lineno="159"><highlight class="normal"><sp/><sp/><sp/><sp/>}</highlight></codeline>
<codeline lineno="160"><highlight class="normal"><sp/><sp/><sp/><sp/></highlight><highlight class="comment">/*<sp/>Step<sp/>5:<sp/>no<sp/>quotient<sp/>test<sp/>if<sp/>x<sp/>not<sp/>positive<sp/>*/</highlight><highlight class="normal"></highlight></codeline>
<codeline lineno="161"><highlight class="normal"></highlight></codeline>
<codeline lineno="162"><highlight class="normal"><sp/><sp/><sp/><sp/></highlight><highlight class="keywordflow">if</highlight><highlight class="normal"><sp/>(x<sp/>&gt;<sp/>0.0)<sp/>{</highlight></codeline>
<codeline lineno="163"><highlight class="normal"><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/></highlight><highlight class="comment">/*<sp/>Step<sp/>6:<sp/>calculation<sp/>of<sp/>v<sp/>and<sp/>quotient<sp/>q<sp/>*/</highlight><highlight class="normal"></highlight></codeline>
<codeline lineno="164"><highlight class="normal"><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/>v<sp/>=<sp/>t<sp/>/<sp/>(s<sp/>+<sp/>s);</highlight></codeline>
<codeline lineno="165"><highlight class="normal"><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/></highlight><highlight class="keywordflow">if</highlight><highlight class="normal"><sp/>(fabs(v)<sp/>&lt;=<sp/>0.25)</highlight></codeline>
<codeline lineno="166"><highlight class="normal"><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/>q<sp/>=<sp/>q0<sp/>+<sp/>0.5<sp/>*<sp/>t<sp/>*<sp/>t<sp/>*<sp/>((((((a7<sp/>*<sp/>v<sp/>+<sp/>a6)<sp/>*<sp/>v<sp/>+<sp/>a5)<sp/>*<sp/>v<sp/>+<sp/>a4)<sp/>*<sp/>v</highlight></codeline>
<codeline lineno="167"><highlight class="normal"><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/>+<sp/>a3)<sp/>*<sp/>v<sp/>+<sp/>a2)<sp/>*<sp/>v<sp/>+<sp/>a1)<sp/>*<sp/>v;</highlight></codeline>
<codeline lineno="168"><highlight class="normal"><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/></highlight><highlight class="keywordflow">else</highlight><highlight class="normal"></highlight></codeline>
<codeline lineno="169"><highlight class="normal"><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/>q<sp/>=<sp/>q0<sp/>-<sp/>s<sp/>*<sp/>t<sp/>+<sp/>0.25<sp/>*<sp/>t<sp/>*<sp/>t<sp/>+<sp/>(s2<sp/>+<sp/>s2)<sp/>*<sp/>log(1.0<sp/>+<sp/>v);</highlight></codeline>
<codeline lineno="170"><highlight class="normal"></highlight></codeline>
<codeline lineno="171"><highlight class="normal"></highlight></codeline>
<codeline lineno="172"><highlight class="normal"><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/></highlight><highlight class="comment">/*<sp/>Step<sp/>7:<sp/>quotient<sp/>acceptance<sp/>(q)<sp/>*/</highlight><highlight class="normal"></highlight></codeline>
<codeline lineno="173"><highlight class="normal"><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/></highlight><highlight class="keywordflow">if</highlight><highlight class="normal"><sp/>(log(1.0<sp/>-<sp/>u)<sp/>&lt;=<sp/>q)</highlight></codeline>
<codeline lineno="174"><highlight class="normal"><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/></highlight><highlight class="keywordflow">return</highlight><highlight class="normal"><sp/>scale<sp/>*<sp/>ret_val;</highlight></codeline>
<codeline lineno="175"><highlight class="normal"><sp/><sp/><sp/><sp/>}</highlight></codeline>
<codeline lineno="176"><highlight class="normal"></highlight></codeline>
<codeline lineno="177"><highlight class="normal"><sp/><sp/><sp/><sp/>repeat<sp/>{</highlight></codeline>
<codeline lineno="178"><highlight class="normal"><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/></highlight><highlight class="comment">/*<sp/>Step<sp/>8:<sp/>e<sp/>=<sp/>standard<sp/>exponential<sp/>deviate</highlight></codeline>
<codeline lineno="179"><highlight class="comment"><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/>*<sp/><sp/><sp/><sp/><sp/><sp/>u<sp/>=<sp/><sp/>0,1<sp/>-uniform<sp/>deviate</highlight></codeline>
<codeline lineno="180"><highlight class="comment"><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/>*<sp/><sp/><sp/><sp/><sp/><sp/>t<sp/>=<sp/>(b,si)-double<sp/>exponential<sp/>(laplace)<sp/>sample<sp/>*/</highlight><highlight class="normal"></highlight></codeline>
<codeline lineno="181"><highlight class="normal"><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/>e<sp/>=<sp/>exp_rand();</highlight></codeline>
<codeline lineno="182"><highlight class="normal"><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/>u<sp/>=<sp/>unif_rand();</highlight></codeline>
<codeline lineno="183"><highlight class="normal"><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/>u<sp/>=<sp/>u<sp/>+<sp/>u<sp/>-<sp/>1.0;</highlight></codeline>
<codeline lineno="184"><highlight class="normal"><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/></highlight><highlight class="keywordflow">if</highlight><highlight class="normal"><sp/>(u<sp/>&lt;<sp/>0.0)</highlight></codeline>
<codeline lineno="185"><highlight class="normal"><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/>t<sp/>=<sp/>b<sp/>-<sp/>si<sp/>*<sp/>e;</highlight></codeline>
<codeline lineno="186"><highlight class="normal"><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/></highlight><highlight class="keywordflow">else</highlight><highlight class="normal"></highlight></codeline>
<codeline lineno="187"><highlight class="normal"><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/>t<sp/>=<sp/>b<sp/>+<sp/>si<sp/>*<sp/>e;</highlight></codeline>
<codeline lineno="188"><highlight class="normal"><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/></highlight><highlight class="comment">/*<sp/>Step<sp/><sp/>9:<sp/><sp/>rejection<sp/>if<sp/>t<sp/>&lt;<sp/>tau(1)<sp/>=<sp/>-0.71874483771719<sp/>*/</highlight><highlight class="normal"></highlight></codeline>
<codeline lineno="189"><highlight class="normal"><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/></highlight><highlight class="keywordflow">if</highlight><highlight class="normal"><sp/>(t<sp/>&gt;=<sp/>-0.71874483771719)<sp/>{</highlight></codeline>
<codeline lineno="190"><highlight class="normal"><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/></highlight><highlight class="comment">/*<sp/>Step<sp/>10:<sp/><sp/>calculation<sp/>of<sp/>v<sp/>and<sp/>quotient<sp/>q<sp/>*/</highlight><highlight class="normal"></highlight></codeline>
<codeline lineno="191"><highlight class="normal"><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/>v<sp/>=<sp/>t<sp/>/<sp/>(s<sp/>+<sp/>s);</highlight></codeline>
<codeline lineno="192"><highlight class="normal"><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/></highlight><highlight class="keywordflow">if</highlight><highlight class="normal"><sp/>(fabs(v)<sp/>&lt;=<sp/>0.25)</highlight></codeline>
<codeline lineno="193"><highlight class="normal"><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/>q<sp/>=<sp/>q0<sp/>+<sp/>0.5<sp/>*<sp/>t<sp/>*<sp/>t<sp/>*</highlight></codeline>
<codeline lineno="194"><highlight class="normal"><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/>((((((a7<sp/>*<sp/>v<sp/>+<sp/>a6)<sp/>*<sp/>v<sp/>+<sp/>a5)<sp/>*<sp/>v<sp/>+<sp/>a4)<sp/>*<sp/>v<sp/>+<sp/>a3)<sp/>*<sp/>v</highlight></codeline>
<codeline lineno="195"><highlight class="normal"><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/>+<sp/>a2)<sp/>*<sp/>v<sp/>+<sp/>a1)<sp/>*<sp/>v;</highlight></codeline>
<codeline lineno="196"><highlight class="normal"><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/></highlight><highlight class="keywordflow">else</highlight><highlight class="normal"></highlight></codeline>
<codeline lineno="197"><highlight class="normal"><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/>q<sp/>=<sp/>q0<sp/>-<sp/>s<sp/>*<sp/>t<sp/>+<sp/>0.25<sp/>*<sp/>t<sp/>*<sp/>t<sp/>+<sp/>(s2<sp/>+<sp/>s2)<sp/>*<sp/>log(1.0<sp/>+<sp/>v);</highlight></codeline>
<codeline lineno="198"><highlight class="normal"><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/></highlight><highlight class="comment">/*<sp/>Step<sp/>11:<sp/><sp/>hat<sp/>acceptance<sp/>(h)<sp/>*/</highlight><highlight class="normal"></highlight></codeline>
<codeline lineno="199"><highlight class="normal"><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/></highlight><highlight class="comment">/*<sp/>(if<sp/>q<sp/>not<sp/>positive<sp/>go<sp/>to<sp/>step<sp/>8)<sp/>*/</highlight><highlight class="normal"></highlight></codeline>
<codeline lineno="200"><highlight class="normal"><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/></highlight><highlight class="keywordflow">if</highlight><highlight class="normal"><sp/>(q<sp/>&gt;<sp/>0.0)<sp/>{</highlight></codeline>
<codeline lineno="201"><highlight class="normal"><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/>w<sp/>=<sp/>expm1(q);</highlight></codeline>
<codeline lineno="202"><highlight class="normal"><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/></highlight><highlight class="comment">/*<sp/><sp/>^^^^^<sp/>original<sp/>code<sp/>had<sp/>approximation<sp/>with<sp/>rel.err<sp/>&lt;<sp/>2e-7<sp/>*/</highlight><highlight class="normal"></highlight></codeline>
<codeline lineno="203"><highlight class="normal"><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/></highlight><highlight class="comment">/*<sp/>if<sp/>t<sp/>is<sp/>rejected<sp/>sample<sp/>again<sp/>at<sp/>step<sp/>8<sp/>*/</highlight><highlight class="normal"></highlight></codeline>
<codeline lineno="204"><highlight class="normal"><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/></highlight><highlight class="keywordflow">if</highlight><highlight class="normal"><sp/>(c<sp/>*<sp/>fabs(u)<sp/>&lt;=<sp/>w<sp/>*<sp/>exp(e<sp/>-<sp/>0.5<sp/>*<sp/>t<sp/>*<sp/>t))</highlight></codeline>
<codeline lineno="205"><highlight class="normal"><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/></highlight><highlight class="keywordflow">break</highlight><highlight class="normal">;</highlight></codeline>
<codeline lineno="206"><highlight class="normal"><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/>}</highlight></codeline>
<codeline lineno="207"><highlight class="normal"><sp/><sp/><sp/><sp/><sp/><sp/><sp/><sp/>}</highlight></codeline>
<codeline lineno="208"><highlight class="normal"><sp/><sp/><sp/><sp/>}<sp/></highlight><highlight class="comment">/*<sp/>repeat<sp/>..<sp/>until<sp/><sp/>`t&apos;<sp/>is<sp/>accepted<sp/>*/</highlight><highlight class="normal"></highlight></codeline>
<codeline lineno="209"><highlight class="normal"><sp/><sp/><sp/><sp/>x<sp/>=<sp/>s<sp/>+<sp/>0.5<sp/>*<sp/>t;</highlight></codeline>
<codeline lineno="210"><highlight class="normal"><sp/><sp/><sp/><sp/></highlight><highlight class="keywordflow">return</highlight><highlight class="normal"><sp/>scale<sp/>*<sp/>x<sp/>*<sp/>x;</highlight></codeline>
<codeline lineno="211"><highlight class="normal">}</highlight></codeline>
    </programlisting>
    <location file="/home/smidl/work/git/mixpp/bdm/math/R_rgamma.c"/>
  </compounddef>
</doxygen>
