Added: dev/commons/statistics/1.0-RC1/site/xref/org/apache/commons/statistics/distribution/ExponentialDistribution.html ============================================================================== --- dev/commons/statistics/1.0-RC1/site/xref/org/apache/commons/statistics/distribution/ExponentialDistribution.html (added) +++ dev/commons/statistics/1.0-RC1/site/xref/org/apache/commons/statistics/distribution/ExponentialDistribution.html Thu Dec 1 16:47:12 2022 @@ -0,0 +1,209 @@ +<!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"; xml:lang="en" lang="en"> +<head><meta http-equiv="content-type" content="text/html; charset=UTF-8" /> +<title>ExponentialDistribution xref</title> +<link type="text/css" rel="stylesheet" href="../../../../../stylesheet.css" /> +</head> +<body> +<div id="overview"><a href="../../../../../../apidocs/org/apache/commons/statistics/distribution/ExponentialDistribution.html">View Javadoc</a></div><pre> +<a class="jxr_linenumber" name="L1" href="#L1">1</a> <em class="jxr_comment">/*</em> +<a class="jxr_linenumber" name="L2" href="#L2">2</a> <em class="jxr_comment"> * Licensed to the Apache Software Foundation (ASF) under one or more</em> +<a class="jxr_linenumber" name="L3" href="#L3">3</a> <em class="jxr_comment"> * contributor license agreements. See the NOTICE file distributed with</em> +<a class="jxr_linenumber" name="L4" href="#L4">4</a> <em class="jxr_comment"> * this work for additional information regarding copyright ownership.</em> +<a class="jxr_linenumber" name="L5" href="#L5">5</a> <em class="jxr_comment"> * The ASF licenses this file to You under the Apache License, Version 2.0</em> +<a class="jxr_linenumber" name="L6" href="#L6">6</a> <em class="jxr_comment"> * (the "License"); you may not use this file except in compliance with</em> +<a class="jxr_linenumber" name="L7" href="#L7">7</a> <em class="jxr_comment"> * the License. You may obtain a copy of the License at</em> +<a class="jxr_linenumber" name="L8" href="#L8">8</a> <em class="jxr_comment"> *</em> +<a class="jxr_linenumber" name="L9" href="#L9">9</a> <em class="jxr_comment"> * <a href="http://www.apache.org/licenses/LICENSE-2.0"; target="alexandria_uri">http://www.apache.org/licenses/LICENSE-2.0</a></em> +<a class="jxr_linenumber" name="L10" href="#L10">10</a> <em class="jxr_comment"> *</em> +<a class="jxr_linenumber" name="L11" href="#L11">11</a> <em class="jxr_comment"> * Unless required by applicable law or agreed to in writing, software</em> +<a class="jxr_linenumber" name="L12" href="#L12">12</a> <em class="jxr_comment"> * distributed under the License is distributed on an "AS IS" BASIS,</em> +<a class="jxr_linenumber" name="L13" href="#L13">13</a> <em class="jxr_comment"> * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.</em> +<a class="jxr_linenumber" name="L14" href="#L14">14</a> <em class="jxr_comment"> * See the License for the specific language governing permissions and</em> +<a class="jxr_linenumber" name="L15" href="#L15">15</a> <em class="jxr_comment"> * limitations under the License.</em> +<a class="jxr_linenumber" name="L16" href="#L16">16</a> <em class="jxr_comment"> */</em> +<a class="jxr_linenumber" name="L17" href="#L17">17</a> <strong class="jxr_keyword">package</strong> org.apache.commons.statistics.distribution; +<a class="jxr_linenumber" name="L18" href="#L18">18</a> +<a class="jxr_linenumber" name="L19" href="#L19">19</a> <strong class="jxr_keyword">import</strong> org.apache.commons.rng.UniformRandomProvider; +<a class="jxr_linenumber" name="L20" href="#L20">20</a> <strong class="jxr_keyword">import</strong> org.apache.commons.rng.sampling.distribution.ZigguratSampler; +<a class="jxr_linenumber" name="L21" href="#L21">21</a> +<a class="jxr_linenumber" name="L22" href="#L22">22</a> <em class="jxr_javadoccomment">/**</em> +<a class="jxr_linenumber" name="L23" href="#L23">23</a> <em class="jxr_javadoccomment"> * Implementation of the exponential distribution.</em> +<a class="jxr_linenumber" name="L24" href="#L24">24</a> <em class="jxr_javadoccomment"> *</em> +<a class="jxr_linenumber" name="L25" href="#L25">25</a> <em class="jxr_javadoccomment"> * <p>The probability density function of \( X \) is:</em> +<a class="jxr_linenumber" name="L26" href="#L26">26</a> <em class="jxr_javadoccomment"> *</em> +<a class="jxr_linenumber" name="L27" href="#L27">27</a> <em class="jxr_javadoccomment"> * <p>\[ f(x; \mu) = \frac{1}{\mu} e^{-x / \mu} \]</em> +<a class="jxr_linenumber" name="L28" href="#L28">28</a> <em class="jxr_javadoccomment"> *</em> +<a class="jxr_linenumber" name="L29" href="#L29">29</a> <em class="jxr_javadoccomment"> * <p>for \( \mu &gt; 0 \) the mean and</em> +<a class="jxr_linenumber" name="L30" href="#L30">30</a> <em class="jxr_javadoccomment"> * \( x \in [0, \infty) \).</em> +<a class="jxr_linenumber" name="L31" href="#L31">31</a> <em class="jxr_javadoccomment"> *</em> +<a class="jxr_linenumber" name="L32" href="#L32">32</a> <em class="jxr_javadoccomment"> * <p>This implementation uses the scale parameter \( \mu \) which is the mean of the distribution.</em> +<a class="jxr_linenumber" name="L33" href="#L33">33</a> <em class="jxr_javadoccomment"> * A common alternative parameterization uses the rate parameter \( \lambda \) which is the reciprocal</em> +<a class="jxr_linenumber" name="L34" href="#L34">34</a> <em class="jxr_javadoccomment"> * of the mean. The distribution can be be created using \( \mu = \frac{1}{\lambda} \).</em> +<a class="jxr_linenumber" name="L35" href="#L35">35</a> <em class="jxr_javadoccomment"> *</em> +<a class="jxr_linenumber" name="L36" href="#L36">36</a> <em class="jxr_javadoccomment"> * @see <a href="<a href="https://en.wikipedia.org/wiki/Exponential_distribution"; target="alexandria_uri">https://en.wikipedia.org/wiki/Exponential_distribution</a>">Exponential distribution (Wikipedia)</a></em> +<a class="jxr_linenumber" name="L37" href="#L37">37</a> <em class="jxr_javadoccomment"> * @see <a href="<a href="https://mathworld.wolfram.com/ExponentialDistribution.html"; target="alexandria_uri">https://mathworld.wolfram.com/ExponentialDistribution.html</a>">Exponential distribution (MathWorld)</a></em> +<a class="jxr_linenumber" name="L38" href="#L38">38</a> <em class="jxr_javadoccomment"> */</em> +<a class="jxr_linenumber" name="L39" href="#L39">39</a> <strong class="jxr_keyword">public</strong> <strong class="jxr_keyword">final</strong> <strong class="jxr_keyword">class</strong> <a name="ExponentialDistribution" href="../../../../../org/apache/commons/statistics/distribution/ExponentialDistribution.html#ExponentialDistribution">ExponentialDistribution</a> <strong class="jxr_keyword">extends</strong> <a name="AbstractContinuousDistribution" href="../../../../../org/apache/commons/statistics/distribution/AbstractContinuousDistribution.html#AbstractContinuousDistribution">AbstractContinuousDistribution</a> { +<a class="jxr_linenumber" name="L40" href="#L40">40</a> <em class="jxr_javadoccomment">/** Support lower bound. */</em> +<a class="jxr_linenumber" name="L41" href="#L41">41</a> <strong class="jxr_keyword">private</strong> <strong class="jxr_keyword">static</strong> <strong class="jxr_keyword">final</strong> <strong class="jxr_keyword">double</strong> SUPPORT_LO = 0; +<a class="jxr_linenumber" name="L42" href="#L42">42</a> <em class="jxr_javadoccomment">/** Support upper bound. */</em> +<a class="jxr_linenumber" name="L43" href="#L43">43</a> <strong class="jxr_keyword">private</strong> <strong class="jxr_keyword">static</strong> <strong class="jxr_keyword">final</strong> <strong class="jxr_keyword">double</strong> SUPPORT_HI = Double.POSITIVE_INFINITY; +<a class="jxr_linenumber" name="L44" href="#L44">44</a> <em class="jxr_javadoccomment">/** ln(2). */</em> +<a class="jxr_linenumber" name="L45" href="#L45">45</a> <strong class="jxr_keyword">private</strong> <strong class="jxr_keyword">static</strong> <strong class="jxr_keyword">final</strong> <strong class="jxr_keyword">double</strong> LN_2 = 0.6931471805599453094172; +<a class="jxr_linenumber" name="L46" href="#L46">46</a> <em class="jxr_javadoccomment">/** The mean of this distribution. */</em> +<a class="jxr_linenumber" name="L47" href="#L47">47</a> <strong class="jxr_keyword">private</strong> <strong class="jxr_keyword">final</strong> <strong class="jxr_keyword">double</strong> mean; +<a class="jxr_linenumber" name="L48" href="#L48">48</a> <em class="jxr_javadoccomment">/** The logarithm of the mean, stored to reduce computing time. */</em> +<a class="jxr_linenumber" name="L49" href="#L49">49</a> <strong class="jxr_keyword">private</strong> <strong class="jxr_keyword">final</strong> <strong class="jxr_keyword">double</strong> logMean; +<a class="jxr_linenumber" name="L50" href="#L50">50</a> +<a class="jxr_linenumber" name="L51" href="#L51">51</a> <em class="jxr_javadoccomment">/**</em> +<a class="jxr_linenumber" name="L52" href="#L52">52</a> <em class="jxr_javadoccomment"> * @param mean Mean of this distribution.</em> +<a class="jxr_linenumber" name="L53" href="#L53">53</a> <em class="jxr_javadoccomment"> */</em> +<a class="jxr_linenumber" name="L54" href="#L54">54</a> <strong class="jxr_keyword">private</strong> <a name="ExponentialDistribution" href="../../../../../org/apache/commons/statistics/distribution/ExponentialDistribution.html#ExponentialDistribution">ExponentialDistribution</a>(<strong class="jxr_keyword">double</strong> mean) { +<a class="jxr_linenumber" name="L55" href="#L55">55</a> <strong class="jxr_keyword">this</strong>.mean = mean; +<a class="jxr_linenumber" name="L56" href="#L56">56</a> logMean = Math.log(mean); +<a class="jxr_linenumber" name="L57" href="#L57">57</a> } +<a class="jxr_linenumber" name="L58" href="#L58">58</a> +<a class="jxr_linenumber" name="L59" href="#L59">59</a> <em class="jxr_javadoccomment">/**</em> +<a class="jxr_linenumber" name="L60" href="#L60">60</a> <em class="jxr_javadoccomment"> * Creates an exponential distribution.</em> +<a class="jxr_linenumber" name="L61" href="#L61">61</a> <em class="jxr_javadoccomment"> *</em> +<a class="jxr_linenumber" name="L62" href="#L62">62</a> <em class="jxr_javadoccomment"> * @param mean Mean of this distribution. This is a scale parameter.</em> +<a class="jxr_linenumber" name="L63" href="#L63">63</a> <em class="jxr_javadoccomment"> * @return the distribution</em> +<a class="jxr_linenumber" name="L64" href="#L64">64</a> <em class="jxr_javadoccomment"> * @throws IllegalArgumentException if {@code mean <= 0}.</em> +<a class="jxr_linenumber" name="L65" href="#L65">65</a> <em class="jxr_javadoccomment"> */</em> +<a class="jxr_linenumber" name="L66" href="#L66">66</a> <strong class="jxr_keyword">public</strong> <strong class="jxr_keyword">static</strong> <a name="ExponentialDistribution" href="../../../../../org/apache/commons/statistics/distribution/ExponentialDistribution.html#ExponentialDistribution">ExponentialDistribution</a> of(<strong class="jxr_keyword">double</strong> mean) { +<a class="jxr_linenumber" name="L67" href="#L67">67</a> <strong class="jxr_keyword">if</strong> (mean <= 0) { +<a class="jxr_linenumber" name="L68" href="#L68">68</a> <strong class="jxr_keyword">throw</strong> <strong class="jxr_keyword">new</strong> <a name="DistributionException" href="../../../../../org/apache/commons/statistics/distribution/DistributionException.html#DistributionException">DistributionException</a>(DistributionException.NOT_STRICTLY_POSITIVE, mean); +<a class="jxr_linenumber" name="L69" href="#L69">69</a> } +<a class="jxr_linenumber" name="L70" href="#L70">70</a> <strong class="jxr_keyword">return</strong> <strong class="jxr_keyword">new</strong> <a name="ExponentialDistribution" href="../../../../../org/apache/commons/statistics/distribution/ExponentialDistribution.html#ExponentialDistribution">ExponentialDistribution</a>(mean); +<a class="jxr_linenumber" name="L71" href="#L71">71</a> } +<a class="jxr_linenumber" name="L72" href="#L72">72</a> +<a class="jxr_linenumber" name="L73" href="#L73">73</a> <em class="jxr_javadoccomment">/** {@inheritDoc} */</em> +<a class="jxr_linenumber" name="L74" href="#L74">74</a> @Override +<a class="jxr_linenumber" name="L75" href="#L75">75</a> <strong class="jxr_keyword">public</strong> <strong class="jxr_keyword">double</strong> density(<strong class="jxr_keyword">double</strong> x) { +<a class="jxr_linenumber" name="L76" href="#L76">76</a> <strong class="jxr_keyword">if</strong> (x < SUPPORT_LO) { +<a class="jxr_linenumber" name="L77" href="#L77">77</a> <strong class="jxr_keyword">return</strong> 0; +<a class="jxr_linenumber" name="L78" href="#L78">78</a> } +<a class="jxr_linenumber" name="L79" href="#L79">79</a> <strong class="jxr_keyword">return</strong> Math.exp(-x / mean) / mean; +<a class="jxr_linenumber" name="L80" href="#L80">80</a> } +<a class="jxr_linenumber" name="L81" href="#L81">81</a> +<a class="jxr_linenumber" name="L82" href="#L82">82</a> <em class="jxr_javadoccomment">/** {@inheritDoc} **/</em> +<a class="jxr_linenumber" name="L83" href="#L83">83</a> @Override +<a class="jxr_linenumber" name="L84" href="#L84">84</a> <strong class="jxr_keyword">public</strong> <strong class="jxr_keyword">double</strong> logDensity(<strong class="jxr_keyword">double</strong> x) { +<a class="jxr_linenumber" name="L85" href="#L85">85</a> <strong class="jxr_keyword">if</strong> (x < SUPPORT_LO) { +<a class="jxr_linenumber" name="L86" href="#L86">86</a> <strong class="jxr_keyword">return</strong> Double.NEGATIVE_INFINITY; +<a class="jxr_linenumber" name="L87" href="#L87">87</a> } +<a class="jxr_linenumber" name="L88" href="#L88">88</a> <strong class="jxr_keyword">return</strong> -x / mean - logMean; +<a class="jxr_linenumber" name="L89" href="#L89">89</a> } +<a class="jxr_linenumber" name="L90" href="#L90">90</a> +<a class="jxr_linenumber" name="L91" href="#L91">91</a> <em class="jxr_javadoccomment">/** {@inheritDoc} */</em> +<a class="jxr_linenumber" name="L92" href="#L92">92</a> @Override +<a class="jxr_linenumber" name="L93" href="#L93">93</a> <strong class="jxr_keyword">public</strong> <strong class="jxr_keyword">double</strong> cumulativeProbability(<strong class="jxr_keyword">double</strong> x) { +<a class="jxr_linenumber" name="L94" href="#L94">94</a> <strong class="jxr_keyword">if</strong> (x <= SUPPORT_LO) { +<a class="jxr_linenumber" name="L95" href="#L95">95</a> <strong class="jxr_keyword">return</strong> 0; +<a class="jxr_linenumber" name="L96" href="#L96">96</a> } +<a class="jxr_linenumber" name="L97" href="#L97">97</a> <strong class="jxr_keyword">return</strong> -Math.expm1(-x / mean); +<a class="jxr_linenumber" name="L98" href="#L98">98</a> } +<a class="jxr_linenumber" name="L99" href="#L99">99</a> +<a class="jxr_linenumber" name="L100" href="#L100">100</a> <em class="jxr_javadoccomment">/** {@inheritDoc} */</em> +<a class="jxr_linenumber" name="L101" href="#L101">101</a> @Override +<a class="jxr_linenumber" name="L102" href="#L102">102</a> <strong class="jxr_keyword">public</strong> <strong class="jxr_keyword">double</strong> survivalProbability(<strong class="jxr_keyword">double</strong> x) { +<a class="jxr_linenumber" name="L103" href="#L103">103</a> <strong class="jxr_keyword">if</strong> (x <= SUPPORT_LO) { +<a class="jxr_linenumber" name="L104" href="#L104">104</a> <strong class="jxr_keyword">return</strong> 1; +<a class="jxr_linenumber" name="L105" href="#L105">105</a> } +<a class="jxr_linenumber" name="L106" href="#L106">106</a> <strong class="jxr_keyword">return</strong> Math.exp(-x / mean); +<a class="jxr_linenumber" name="L107" href="#L107">107</a> } +<a class="jxr_linenumber" name="L108" href="#L108">108</a> +<a class="jxr_linenumber" name="L109" href="#L109">109</a> <em class="jxr_javadoccomment">/**</em> +<a class="jxr_linenumber" name="L110" href="#L110">110</a> <em class="jxr_javadoccomment"> * {@inheritDoc}</em> +<a class="jxr_linenumber" name="L111" href="#L111">111</a> <em class="jxr_javadoccomment"> *</em> +<a class="jxr_linenumber" name="L112" href="#L112">112</a> <em class="jxr_javadoccomment"> * <p>Returns {@code 0} when {@code p == 0} and</em> +<a class="jxr_linenumber" name="L113" href="#L113">113</a> <em class="jxr_javadoccomment"> * {@link Double#POSITIVE_INFINITY} when {@code p == 1}.</em> +<a class="jxr_linenumber" name="L114" href="#L114">114</a> <em class="jxr_javadoccomment"> */</em> +<a class="jxr_linenumber" name="L115" href="#L115">115</a> @Override +<a class="jxr_linenumber" name="L116" href="#L116">116</a> <strong class="jxr_keyword">public</strong> <strong class="jxr_keyword">double</strong> inverseCumulativeProbability(<strong class="jxr_keyword">double</strong> p) { +<a class="jxr_linenumber" name="L117" href="#L117">117</a> ArgumentUtils.checkProbability(p); +<a class="jxr_linenumber" name="L118" href="#L118">118</a> <strong class="jxr_keyword">if</strong> (p == 1) { +<a class="jxr_linenumber" name="L119" href="#L119">119</a> <strong class="jxr_keyword">return</strong> Double.POSITIVE_INFINITY; +<a class="jxr_linenumber" name="L120" href="#L120">120</a> } +<a class="jxr_linenumber" name="L121" href="#L121">121</a> <em class="jxr_comment">// Subtract from zero to prevent returning -0.0 for p=-0.0</em> +<a class="jxr_linenumber" name="L122" href="#L122">122</a> <strong class="jxr_keyword">return</strong> 0 - mean * Math.log1p(-p); +<a class="jxr_linenumber" name="L123" href="#L123">123</a> } +<a class="jxr_linenumber" name="L124" href="#L124">124</a> +<a class="jxr_linenumber" name="L125" href="#L125">125</a> <em class="jxr_javadoccomment">/**</em> +<a class="jxr_linenumber" name="L126" href="#L126">126</a> <em class="jxr_javadoccomment"> * {@inheritDoc}</em> +<a class="jxr_linenumber" name="L127" href="#L127">127</a> <em class="jxr_javadoccomment"> *</em> +<a class="jxr_linenumber" name="L128" href="#L128">128</a> <em class="jxr_javadoccomment"> * <p>Returns {@code 0} when {@code p == 1} and</em> +<a class="jxr_linenumber" name="L129" href="#L129">129</a> <em class="jxr_javadoccomment"> * {@link Double#POSITIVE_INFINITY} when {@code p == 0}.</em> +<a class="jxr_linenumber" name="L130" href="#L130">130</a> <em class="jxr_javadoccomment"> */</em> +<a class="jxr_linenumber" name="L131" href="#L131">131</a> @Override +<a class="jxr_linenumber" name="L132" href="#L132">132</a> <strong class="jxr_keyword">public</strong> <strong class="jxr_keyword">double</strong> inverseSurvivalProbability(<strong class="jxr_keyword">double</strong> p) { +<a class="jxr_linenumber" name="L133" href="#L133">133</a> ArgumentUtils.checkProbability(p); +<a class="jxr_linenumber" name="L134" href="#L134">134</a> <strong class="jxr_keyword">if</strong> (p == 0) { +<a class="jxr_linenumber" name="L135" href="#L135">135</a> <strong class="jxr_keyword">return</strong> Double.POSITIVE_INFINITY; +<a class="jxr_linenumber" name="L136" href="#L136">136</a> } +<a class="jxr_linenumber" name="L137" href="#L137">137</a> <em class="jxr_comment">// Subtract from zero to prevent returning -0.0 for p=1</em> +<a class="jxr_linenumber" name="L138" href="#L138">138</a> <strong class="jxr_keyword">return</strong> 0 - mean * Math.log(p); +<a class="jxr_linenumber" name="L139" href="#L139">139</a> } +<a class="jxr_linenumber" name="L140" href="#L140">140</a> +<a class="jxr_linenumber" name="L141" href="#L141">141</a> <em class="jxr_javadoccomment">/** {@inheritDoc} */</em> +<a class="jxr_linenumber" name="L142" href="#L142">142</a> @Override +<a class="jxr_linenumber" name="L143" href="#L143">143</a> <strong class="jxr_keyword">public</strong> <strong class="jxr_keyword">double</strong> getMean() { +<a class="jxr_linenumber" name="L144" href="#L144">144</a> <strong class="jxr_keyword">return</strong> mean; +<a class="jxr_linenumber" name="L145" href="#L145">145</a> } +<a class="jxr_linenumber" name="L146" href="#L146">146</a> +<a class="jxr_linenumber" name="L147" href="#L147">147</a> <em class="jxr_javadoccomment">/**</em> +<a class="jxr_linenumber" name="L148" href="#L148">148</a> <em class="jxr_javadoccomment"> * {@inheritDoc}</em> +<a class="jxr_linenumber" name="L149" href="#L149">149</a> <em class="jxr_javadoccomment"> *</em> +<a class="jxr_linenumber" name="L150" href="#L150">150</a> <em class="jxr_javadoccomment"> * <p>For mean \( \mu \), the variance is \( \mu^2 \).</em> +<a class="jxr_linenumber" name="L151" href="#L151">151</a> <em class="jxr_javadoccomment"> */</em> +<a class="jxr_linenumber" name="L152" href="#L152">152</a> @Override +<a class="jxr_linenumber" name="L153" href="#L153">153</a> <strong class="jxr_keyword">public</strong> <strong class="jxr_keyword">double</strong> getVariance() { +<a class="jxr_linenumber" name="L154" href="#L154">154</a> <strong class="jxr_keyword">return</strong> mean * mean; +<a class="jxr_linenumber" name="L155" href="#L155">155</a> } +<a class="jxr_linenumber" name="L156" href="#L156">156</a> +<a class="jxr_linenumber" name="L157" href="#L157">157</a> <em class="jxr_javadoccomment">/**</em> +<a class="jxr_linenumber" name="L158" href="#L158">158</a> <em class="jxr_javadoccomment"> * {@inheritDoc}</em> +<a class="jxr_linenumber" name="L159" href="#L159">159</a> <em class="jxr_javadoccomment"> *</em> +<a class="jxr_linenumber" name="L160" href="#L160">160</a> <em class="jxr_javadoccomment"> * <p>The lower bound of the support is always 0.</em> +<a class="jxr_linenumber" name="L161" href="#L161">161</a> <em class="jxr_javadoccomment"> *</em> +<a class="jxr_linenumber" name="L162" href="#L162">162</a> <em class="jxr_javadoccomment"> * @return 0.</em> +<a class="jxr_linenumber" name="L163" href="#L163">163</a> <em class="jxr_javadoccomment"> */</em> +<a class="jxr_linenumber" name="L164" href="#L164">164</a> @Override +<a class="jxr_linenumber" name="L165" href="#L165">165</a> <strong class="jxr_keyword">public</strong> <strong class="jxr_keyword">double</strong> getSupportLowerBound() { +<a class="jxr_linenumber" name="L166" href="#L166">166</a> <strong class="jxr_keyword">return</strong> SUPPORT_LO; +<a class="jxr_linenumber" name="L167" href="#L167">167</a> } +<a class="jxr_linenumber" name="L168" href="#L168">168</a> +<a class="jxr_linenumber" name="L169" href="#L169">169</a> <em class="jxr_javadoccomment">/**</em> +<a class="jxr_linenumber" name="L170" href="#L170">170</a> <em class="jxr_javadoccomment"> * {@inheritDoc}</em> +<a class="jxr_linenumber" name="L171" href="#L171">171</a> <em class="jxr_javadoccomment"> *</em> +<a class="jxr_linenumber" name="L172" href="#L172">172</a> <em class="jxr_javadoccomment"> * <p>The upper bound of the support is always positive infinity.</em> +<a class="jxr_linenumber" name="L173" href="#L173">173</a> <em class="jxr_javadoccomment"> *</em> +<a class="jxr_linenumber" name="L174" href="#L174">174</a> <em class="jxr_javadoccomment"> * @return {@link Double#POSITIVE_INFINITY positive infinity}.</em> +<a class="jxr_linenumber" name="L175" href="#L175">175</a> <em class="jxr_javadoccomment"> */</em> +<a class="jxr_linenumber" name="L176" href="#L176">176</a> @Override +<a class="jxr_linenumber" name="L177" href="#L177">177</a> <strong class="jxr_keyword">public</strong> <strong class="jxr_keyword">double</strong> getSupportUpperBound() { +<a class="jxr_linenumber" name="L178" href="#L178">178</a> <strong class="jxr_keyword">return</strong> SUPPORT_HI; +<a class="jxr_linenumber" name="L179" href="#L179">179</a> } +<a class="jxr_linenumber" name="L180" href="#L180">180</a> +<a class="jxr_linenumber" name="L181" href="#L181">181</a> <em class="jxr_javadoccomment">/** {@inheritDoc} */</em> +<a class="jxr_linenumber" name="L182" href="#L182">182</a> @Override +<a class="jxr_linenumber" name="L183" href="#L183">183</a> <strong class="jxr_keyword">double</strong> getMedian() { +<a class="jxr_linenumber" name="L184" href="#L184">184</a> <em class="jxr_comment">// Overridden for the probability(double, double) method.</em> +<a class="jxr_linenumber" name="L185" href="#L185">185</a> <em class="jxr_comment">// This is intentionally not a public method.</em> +<a class="jxr_linenumber" name="L186" href="#L186">186</a> <em class="jxr_comment">// ln(2) / rate = mean * ln(2)</em> +<a class="jxr_linenumber" name="L187" href="#L187">187</a> <strong class="jxr_keyword">return</strong> mean * LN_2; +<a class="jxr_linenumber" name="L188" href="#L188">188</a> } +<a class="jxr_linenumber" name="L189" href="#L189">189</a> +<a class="jxr_linenumber" name="L190" href="#L190">190</a> <em class="jxr_javadoccomment">/** {@inheritDoc} */</em> +<a class="jxr_linenumber" name="L191" href="#L191">191</a> @Override +<a class="jxr_linenumber" name="L192" href="#L192">192</a> <strong class="jxr_keyword">public</strong> ContinuousDistribution.Sampler createSampler(<strong class="jxr_keyword">final</strong> UniformRandomProvider rng) { +<a class="jxr_linenumber" name="L193" href="#L193">193</a> <em class="jxr_comment">// Exponential distribution sampler.</em> +<a class="jxr_linenumber" name="L194" href="#L194">194</a> <strong class="jxr_keyword">return</strong> ZigguratSampler.Exponential.of(rng, getMean())::sample; +<a class="jxr_linenumber" name="L195" href="#L195">195</a> } +<a class="jxr_linenumber" name="L196" href="#L196">196</a> } +</pre> +<hr/> +<div id="footer">Copyright © 2018–2022 <a href="https://www.apache.org/";>The Apache Software Foundation</a>. All rights reserved.</div> +</body> +</html>
Added: dev/commons/statistics/1.0-RC1/site/xref/org/apache/commons/statistics/distribution/ExtendedPrecision.html ============================================================================== --- dev/commons/statistics/1.0-RC1/site/xref/org/apache/commons/statistics/distribution/ExtendedPrecision.html (added) +++ dev/commons/statistics/1.0-RC1/site/xref/org/apache/commons/statistics/distribution/ExtendedPrecision.html Thu Dec 1 16:47:12 2022 @@ -0,0 +1,336 @@ +<!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"; xml:lang="en" lang="en"> +<head><meta http-equiv="content-type" content="text/html; charset=UTF-8" /> +<title>ExtendedPrecision xref</title> +<link type="text/css" rel="stylesheet" href="../../../../../stylesheet.css" /> +</head> +<body> +<div id="overview"><a href="../../../../../../apidocs/org/apache/commons/statistics/distribution/ExtendedPrecision.html">View Javadoc</a></div><pre> +<a class="jxr_linenumber" name="L1" href="#L1">1</a> <em class="jxr_comment">/*</em> +<a class="jxr_linenumber" name="L2" href="#L2">2</a> <em class="jxr_comment"> * Licensed to the Apache Software Foundation (ASF) under one or more</em> +<a class="jxr_linenumber" name="L3" href="#L3">3</a> <em class="jxr_comment"> * contributor license agreements. See the NOTICE file distributed with</em> +<a class="jxr_linenumber" name="L4" href="#L4">4</a> <em class="jxr_comment"> * this work for additional information regarding copyright ownership.</em> +<a class="jxr_linenumber" name="L5" href="#L5">5</a> <em class="jxr_comment"> * The ASF licenses this file to You under the Apache License, Version 2.0</em> +<a class="jxr_linenumber" name="L6" href="#L6">6</a> <em class="jxr_comment"> * (the "License"); you may not use this file except in compliance with</em> +<a class="jxr_linenumber" name="L7" href="#L7">7</a> <em class="jxr_comment"> * the License. You may obtain a copy of the License at</em> +<a class="jxr_linenumber" name="L8" href="#L8">8</a> <em class="jxr_comment"> *</em> +<a class="jxr_linenumber" name="L9" href="#L9">9</a> <em class="jxr_comment"> * <a href="http://www.apache.org/licenses/LICENSE-2.0"; target="alexandria_uri">http://www.apache.org/licenses/LICENSE-2.0</a></em> +<a class="jxr_linenumber" name="L10" href="#L10">10</a> <em class="jxr_comment"> *</em> +<a class="jxr_linenumber" name="L11" href="#L11">11</a> <em class="jxr_comment"> * Unless required by applicable law or agreed to in writing, software</em> +<a class="jxr_linenumber" name="L12" href="#L12">12</a> <em class="jxr_comment"> * distributed under the License is distributed on an "AS IS" BASIS,</em> +<a class="jxr_linenumber" name="L13" href="#L13">13</a> <em class="jxr_comment"> * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.</em> +<a class="jxr_linenumber" name="L14" href="#L14">14</a> <em class="jxr_comment"> * See the License for the specific language governing permissions and</em> +<a class="jxr_linenumber" name="L15" href="#L15">15</a> <em class="jxr_comment"> * limitations under the License.</em> +<a class="jxr_linenumber" name="L16" href="#L16">16</a> <em class="jxr_comment"> */</em> +<a class="jxr_linenumber" name="L17" href="#L17">17</a> <strong class="jxr_keyword">package</strong> org.apache.commons.statistics.distribution; +<a class="jxr_linenumber" name="L18" href="#L18">18</a> +<a class="jxr_linenumber" name="L19" href="#L19">19</a> <em class="jxr_javadoccomment">/**</em> +<a class="jxr_linenumber" name="L20" href="#L20">20</a> <em class="jxr_javadoccomment"> * Computes extended precision floating-point operations.</em> +<a class="jxr_linenumber" name="L21" href="#L21">21</a> <em class="jxr_javadoccomment"> *</em> +<a class="jxr_linenumber" name="L22" href="#L22">22</a> <em class="jxr_javadoccomment"> * <p>It is based on the 1971 paper</em> +<a class="jxr_linenumber" name="L23" href="#L23">23</a> <em class="jxr_javadoccomment"> * <a href="<a href="https://doi.org/10.1007/BF01397083"; target="alexandria_uri">https://doi.org/10.1007/BF01397083</a>"></em> +<a class="jxr_linenumber" name="L24" href="#L24">24</a> <em class="jxr_javadoccomment"> * Dekker (1971) A floating-point technique for extending the available precision</a></em> +<a class="jxr_linenumber" name="L25" href="#L25">25</a> <em class="jxr_javadoccomment"> * Numer. Math. 18, 224-242.</em> +<a class="jxr_linenumber" name="L26" href="#L26">26</a> <em class="jxr_javadoccomment"> *</em> +<a class="jxr_linenumber" name="L27" href="#L27">27</a> <em class="jxr_javadoccomment"> * <p>Adapted from {@code org.apache.commons.numbers.core.ExtendedPrecision}.</em> +<a class="jxr_linenumber" name="L28" href="#L28">28</a> <em class="jxr_javadoccomment"> */</em> +<a class="jxr_linenumber" name="L29" href="#L29">29</a> <strong class="jxr_keyword">final</strong> <strong class="jxr_keyword">class</strong> <a name="ExtendedPrecision" href="../../../../../org/apache/commons/statistics/distribution/ExtendedPrecision.html#ExtendedPrecision">ExtendedPrecision</a> { +<a class="jxr_linenumber" name="L30" href="#L30">30</a> <em class="jxr_javadoccomment">/** sqrt(2 pi) as a double. Computed to 64-digits precision and converted to double. */</em> +<a class="jxr_linenumber" name="L31" href="#L31">31</a> <strong class="jxr_keyword">static</strong> <strong class="jxr_keyword">final</strong> <strong class="jxr_keyword">double</strong> SQRT2PI = 2.5066282746310007; +<a class="jxr_linenumber" name="L32" href="#L32">32</a> <em class="jxr_javadoccomment">/** Round-off from sqrt(2 pi) as a double.</em> +<a class="jxr_linenumber" name="L33" href="#L33">33</a> <em class="jxr_javadoccomment"> * Computed from the value sqrt(2 pi) to 64-digits precision minus {@link #SQRT2PI}. */</em> +<a class="jxr_linenumber" name="L34" href="#L34">34</a> <strong class="jxr_keyword">static</strong> <strong class="jxr_keyword">final</strong> <strong class="jxr_keyword">double</strong> SQRT2PI_R = -1.8328579980459167e-16; +<a class="jxr_linenumber" name="L35" href="#L35">35</a> +<a class="jxr_linenumber" name="L36" href="#L36">36</a> <em class="jxr_javadoccomment">/**</em> +<a class="jxr_linenumber" name="L37" href="#L37">37</a> <em class="jxr_javadoccomment"> * The multiplier used to split the double value into high and low parts. From</em> +<a class="jxr_linenumber" name="L38" href="#L38">38</a> <em class="jxr_javadoccomment"> * Dekker (1971): "The constant should be chosen equal to 2^(p - p/2) + 1,</em> +<a class="jxr_linenumber" name="L39" href="#L39">39</a> <em class="jxr_javadoccomment"> * where p is the number of binary digits in the mantissa". Here p is 53</em> +<a class="jxr_linenumber" name="L40" href="#L40">40</a> <em class="jxr_javadoccomment"> * and the multiplier is {@code 2^27 + 1}.</em> +<a class="jxr_linenumber" name="L41" href="#L41">41</a> <em class="jxr_javadoccomment"> */</em> +<a class="jxr_linenumber" name="L42" href="#L42">42</a> <strong class="jxr_keyword">private</strong> <strong class="jxr_keyword">static</strong> <strong class="jxr_keyword">final</strong> <strong class="jxr_keyword">double</strong> MULTIPLIER = 1.0 + 0x1.0p27; +<a class="jxr_linenumber" name="L43" href="#L43">43</a> <em class="jxr_javadoccomment">/** Threshold for a big number that may overflow when squared. 2^500. */</em> +<a class="jxr_linenumber" name="L44" href="#L44">44</a> <strong class="jxr_keyword">private</strong> <strong class="jxr_keyword">static</strong> <strong class="jxr_keyword">final</strong> <strong class="jxr_keyword">double</strong> BIG = 0x1.0p500; +<a class="jxr_linenumber" name="L45" href="#L45">45</a> <em class="jxr_javadoccomment">/** Threshold for a small number that may underflow when squared. 2^-500. */</em> +<a class="jxr_linenumber" name="L46" href="#L46">46</a> <strong class="jxr_keyword">private</strong> <strong class="jxr_keyword">static</strong> <strong class="jxr_keyword">final</strong> <strong class="jxr_keyword">double</strong> SMALL = 0x1.0p-500; +<a class="jxr_linenumber" name="L47" href="#L47">47</a> <em class="jxr_javadoccomment">/** Scale up by 2^600. */</em> +<a class="jxr_linenumber" name="L48" href="#L48">48</a> <strong class="jxr_keyword">private</strong> <strong class="jxr_keyword">static</strong> <strong class="jxr_keyword">final</strong> <strong class="jxr_keyword">double</strong> SCALE_UP = 0x1.0p600; +<a class="jxr_linenumber" name="L49" href="#L49">49</a> <em class="jxr_javadoccomment">/** Scale down by 2^600. */</em> +<a class="jxr_linenumber" name="L50" href="#L50">50</a> <strong class="jxr_keyword">private</strong> <strong class="jxr_keyword">static</strong> <strong class="jxr_keyword">final</strong> <strong class="jxr_keyword">double</strong> SCALE_DOWN = 0x1.0p-600; +<a class="jxr_linenumber" name="L51" href="#L51">51</a> <em class="jxr_javadoccomment">/** Upper bits of sqrt(2 pi). */</em> +<a class="jxr_linenumber" name="L52" href="#L52">52</a> <strong class="jxr_keyword">private</strong> <strong class="jxr_keyword">static</strong> <strong class="jxr_keyword">final</strong> <strong class="jxr_keyword">double</strong> SQRT2PI_H; +<a class="jxr_linenumber" name="L53" href="#L53">53</a> <em class="jxr_javadoccomment">/** Lower bits of sqrt(2 pi). */</em> +<a class="jxr_linenumber" name="L54" href="#L54">54</a> <strong class="jxr_keyword">private</strong> <strong class="jxr_keyword">static</strong> <strong class="jxr_keyword">final</strong> <strong class="jxr_keyword">double</strong> SQRT2PI_L; +<a class="jxr_linenumber" name="L55" href="#L55">55</a> <em class="jxr_javadoccomment">/** X squared value where {@code exp(-0.5*x*x)} cannot increase accuracy using the round-off</em> +<a class="jxr_linenumber" name="L56" href="#L56">56</a> <em class="jxr_javadoccomment"> * from x squared. */</em> +<a class="jxr_linenumber" name="L57" href="#L57">57</a> <strong class="jxr_keyword">private</strong> <strong class="jxr_keyword">static</strong> <strong class="jxr_keyword">final</strong> <strong class="jxr_keyword">int</strong> EXP_M_HALF_XX_MIN_VALUE = 2; +<a class="jxr_linenumber" name="L58" href="#L58">58</a> <em class="jxr_javadoccomment">/** Approximate x squared value where {@code exp(-0.5*x*x) == 0}. This is above</em> +<a class="jxr_linenumber" name="L59" href="#L59">59</a> <em class="jxr_javadoccomment"> * {@code -2 * ln(2^-1074)} due to rounding performed within the exp function. */</em> +<a class="jxr_linenumber" name="L60" href="#L60">60</a> <strong class="jxr_keyword">private</strong> <strong class="jxr_keyword">static</strong> <strong class="jxr_keyword">final</strong> <strong class="jxr_keyword">int</strong> EXP_M_HALF_XX_MAX_VALUE = 1491; +<a class="jxr_linenumber" name="L61" href="#L61">61</a> +<a class="jxr_linenumber" name="L62" href="#L62">62</a> <strong class="jxr_keyword">static</strong> { +<a class="jxr_linenumber" name="L63" href="#L63">63</a> <em class="jxr_comment">// Initialise constants</em> +<a class="jxr_linenumber" name="L64" href="#L64">64</a> +<a class="jxr_linenumber" name="L65" href="#L65">65</a> <em class="jxr_comment">// Split the upper 53-bits for extended precision multiplication</em> +<a class="jxr_linenumber" name="L66" href="#L66">66</a> SQRT2PI_H = highPartUnscaled(SQRT2PI); +<a class="jxr_linenumber" name="L67" href="#L67">67</a> SQRT2PI_L = SQRT2PI - SQRT2PI_H; +<a class="jxr_linenumber" name="L68" href="#L68">68</a> } +<a class="jxr_linenumber" name="L69" href="#L69">69</a> +<a class="jxr_linenumber" name="L70" href="#L70">70</a> <em class="jxr_javadoccomment">/** No instances. */</em> +<a class="jxr_linenumber" name="L71" href="#L71">71</a> <strong class="jxr_keyword">private</strong> <a name="ExtendedPrecision" href="../../../../../org/apache/commons/statistics/distribution/ExtendedPrecision.html#ExtendedPrecision">ExtendedPrecision</a>() {} +<a class="jxr_linenumber" name="L72" href="#L72">72</a> +<a class="jxr_linenumber" name="L73" href="#L73">73</a> <em class="jxr_javadoccomment">/**</em> +<a class="jxr_linenumber" name="L74" href="#L74">74</a> <em class="jxr_javadoccomment"> * Multiply the term by sqrt(2 pi).</em> +<a class="jxr_linenumber" name="L75" href="#L75">75</a> <em class="jxr_javadoccomment"> *</em> +<a class="jxr_linenumber" name="L76" href="#L76">76</a> <em class="jxr_javadoccomment"> * @param x Value (assumed to be positive)</em> +<a class="jxr_linenumber" name="L77" href="#L77">77</a> <em class="jxr_javadoccomment"> * @return x * sqrt(2 pi)</em> +<a class="jxr_linenumber" name="L78" href="#L78">78</a> <em class="jxr_javadoccomment"> */</em> +<a class="jxr_linenumber" name="L79" href="#L79">79</a> <strong class="jxr_keyword">static</strong> <strong class="jxr_keyword">double</strong> xsqrt2pi(<strong class="jxr_keyword">double</strong> x) { +<a class="jxr_linenumber" name="L80" href="#L80">80</a> <em class="jxr_comment">// Note: Do not convert x to absolute for this use case</em> +<a class="jxr_linenumber" name="L81" href="#L81">81</a> <strong class="jxr_keyword">if</strong> (x > BIG) { +<a class="jxr_linenumber" name="L82" href="#L82">82</a> <strong class="jxr_keyword">if</strong> (x == Double.POSITIVE_INFINITY) { +<a class="jxr_linenumber" name="L83" href="#L83">83</a> <strong class="jxr_keyword">return</strong> Double.POSITIVE_INFINITY; +<a class="jxr_linenumber" name="L84" href="#L84">84</a> } +<a class="jxr_linenumber" name="L85" href="#L85">85</a> <strong class="jxr_keyword">return</strong> computeXsqrt2pi(x * SCALE_DOWN) * SCALE_UP; +<a class="jxr_linenumber" name="L86" href="#L86">86</a> } <strong class="jxr_keyword">else</strong> <strong class="jxr_keyword">if</strong> (x < SMALL) { +<a class="jxr_linenumber" name="L87" href="#L87">87</a> <em class="jxr_comment">// Note: Ignore possible zero for this use case</em> +<a class="jxr_linenumber" name="L88" href="#L88">88</a> <strong class="jxr_keyword">return</strong> computeXsqrt2pi(x * SCALE_UP) * SCALE_DOWN; +<a class="jxr_linenumber" name="L89" href="#L89">89</a> } <strong class="jxr_keyword">else</strong> { +<a class="jxr_linenumber" name="L90" href="#L90">90</a> <strong class="jxr_keyword">return</strong> computeXsqrt2pi(x); +<a class="jxr_linenumber" name="L91" href="#L91">91</a> } +<a class="jxr_linenumber" name="L92" href="#L92">92</a> } +<a class="jxr_linenumber" name="L93" href="#L93">93</a> +<a class="jxr_linenumber" name="L94" href="#L94">94</a> <em class="jxr_javadoccomment">/**</em> +<a class="jxr_linenumber" name="L95" href="#L95">95</a> <em class="jxr_javadoccomment"> * Compute {@code a * sqrt(2 * pi)}.</em> +<a class="jxr_linenumber" name="L96" href="#L96">96</a> <em class="jxr_javadoccomment"> *</em> +<a class="jxr_linenumber" name="L97" href="#L97">97</a> <em class="jxr_javadoccomment"> * @param a Value</em> +<a class="jxr_linenumber" name="L98" href="#L98">98</a> <em class="jxr_javadoccomment"> * @return the result</em> +<a class="jxr_linenumber" name="L99" href="#L99">99</a> <em class="jxr_javadoccomment"> */</em> +<a class="jxr_linenumber" name="L100" href="#L100">100</a> <strong class="jxr_keyword">private</strong> <strong class="jxr_keyword">static</strong> <strong class="jxr_keyword">double</strong> computeXsqrt2pi(<strong class="jxr_keyword">double</strong> a) { +<a class="jxr_linenumber" name="L101" href="#L101">101</a> <em class="jxr_comment">// Split the number</em> +<a class="jxr_linenumber" name="L102" href="#L102">102</a> <strong class="jxr_keyword">final</strong> <strong class="jxr_keyword">double</strong> ha = highPartUnscaled(a); +<a class="jxr_linenumber" name="L103" href="#L103">103</a> <strong class="jxr_keyword">final</strong> <strong class="jxr_keyword">double</strong> la = a - ha; +<a class="jxr_linenumber" name="L104" href="#L104">104</a> +<a class="jxr_linenumber" name="L105" href="#L105">105</a> <em class="jxr_comment">// Extended precision product with sqrt(2 * pi)</em> +<a class="jxr_linenumber" name="L106" href="#L106">106</a> <strong class="jxr_keyword">final</strong> <strong class="jxr_keyword">double</strong> x = a * SQRT2PI; +<a class="jxr_linenumber" name="L107" href="#L107">107</a> <strong class="jxr_keyword">final</strong> <strong class="jxr_keyword">double</strong> xx = productLow(ha, la, SQRT2PI_H, SQRT2PI_L, x); +<a class="jxr_linenumber" name="L108" href="#L108">108</a> +<a class="jxr_linenumber" name="L109" href="#L109">109</a> <em class="jxr_comment">// Add the term a multiplied by the round-off from sqrt(2 * pi)</em> +<a class="jxr_linenumber" name="L110" href="#L110">110</a> <em class="jxr_comment">// result = a * (SQRT2PI + SQRT2PI_R)</em> +<a class="jxr_linenumber" name="L111" href="#L111">111</a> <em class="jxr_comment">// Sum from small to high</em> +<a class="jxr_linenumber" name="L112" href="#L112">112</a> <strong class="jxr_keyword">return</strong> a * SQRT2PI_R + xx + x; +<a class="jxr_linenumber" name="L113" href="#L113">113</a> } +<a class="jxr_linenumber" name="L114" href="#L114">114</a> +<a class="jxr_linenumber" name="L115" href="#L115">115</a> <em class="jxr_javadoccomment">/**</em> +<a class="jxr_linenumber" name="L116" href="#L116">116</a> <em class="jxr_javadoccomment"> * Compute {@code sqrt(2 * x * x)}.</em> +<a class="jxr_linenumber" name="L117" href="#L117">117</a> <em class="jxr_javadoccomment"> *</em> +<a class="jxr_linenumber" name="L118" href="#L118">118</a> <em class="jxr_javadoccomment"> * <p>The result is computed using a high precision computation of</em> +<a class="jxr_linenumber" name="L119" href="#L119">119</a> <em class="jxr_javadoccomment"> * {@code sqrt(2 * x * x)} avoiding underflow or overflow of {@code x}</em> +<a class="jxr_linenumber" name="L120" href="#L120">120</a> <em class="jxr_javadoccomment"> * squared.</em> +<a class="jxr_linenumber" name="L121" href="#L121">121</a> <em class="jxr_javadoccomment"> *</em> +<a class="jxr_linenumber" name="L122" href="#L122">122</a> <em class="jxr_javadoccomment"> * @param x Value (assumed to be positive)</em> +<a class="jxr_linenumber" name="L123" href="#L123">123</a> <em class="jxr_javadoccomment"> * @return {@code sqrt(2 * x * x)}</em> +<a class="jxr_linenumber" name="L124" href="#L124">124</a> <em class="jxr_javadoccomment"> */</em> +<a class="jxr_linenumber" name="L125" href="#L125">125</a> <strong class="jxr_keyword">static</strong> <strong class="jxr_keyword">double</strong> sqrt2xx(<strong class="jxr_keyword">double</strong> x) { +<a class="jxr_linenumber" name="L126" href="#L126">126</a> <em class="jxr_comment">// Note: Do not convert x to absolute for this use case</em> +<a class="jxr_linenumber" name="L127" href="#L127">127</a> <strong class="jxr_keyword">if</strong> (x > BIG) { +<a class="jxr_linenumber" name="L128" href="#L128">128</a> <strong class="jxr_keyword">if</strong> (x == Double.POSITIVE_INFINITY) { +<a class="jxr_linenumber" name="L129" href="#L129">129</a> <strong class="jxr_keyword">return</strong> Double.POSITIVE_INFINITY; +<a class="jxr_linenumber" name="L130" href="#L130">130</a> } +<a class="jxr_linenumber" name="L131" href="#L131">131</a> <strong class="jxr_keyword">return</strong> computeSqrt2aa(x * SCALE_DOWN) * SCALE_UP; +<a class="jxr_linenumber" name="L132" href="#L132">132</a> } <strong class="jxr_keyword">else</strong> <strong class="jxr_keyword">if</strong> (x < SMALL) { +<a class="jxr_linenumber" name="L133" href="#L133">133</a> <em class="jxr_comment">// Note: Ignore possible zero for this use case</em> +<a class="jxr_linenumber" name="L134" href="#L134">134</a> <strong class="jxr_keyword">return</strong> computeSqrt2aa(x * SCALE_UP) * SCALE_DOWN; +<a class="jxr_linenumber" name="L135" href="#L135">135</a> } <strong class="jxr_keyword">else</strong> { +<a class="jxr_linenumber" name="L136" href="#L136">136</a> <strong class="jxr_keyword">return</strong> computeSqrt2aa(x); +<a class="jxr_linenumber" name="L137" href="#L137">137</a> } +<a class="jxr_linenumber" name="L138" href="#L138">138</a> } +<a class="jxr_linenumber" name="L139" href="#L139">139</a> +<a class="jxr_linenumber" name="L140" href="#L140">140</a> <em class="jxr_javadoccomment">/**</em> +<a class="jxr_linenumber" name="L141" href="#L141">141</a> <em class="jxr_javadoccomment"> * Compute {@code sqrt(2 * a * a)}.</em> +<a class="jxr_linenumber" name="L142" href="#L142">142</a> <em class="jxr_javadoccomment"> *</em> +<a class="jxr_linenumber" name="L143" href="#L143">143</a> <em class="jxr_javadoccomment"> * @param a Value</em> +<a class="jxr_linenumber" name="L144" href="#L144">144</a> <em class="jxr_javadoccomment"> * @return the result</em> +<a class="jxr_linenumber" name="L145" href="#L145">145</a> <em class="jxr_javadoccomment"> */</em> +<a class="jxr_linenumber" name="L146" href="#L146">146</a> <strong class="jxr_keyword">private</strong> <strong class="jxr_keyword">static</strong> <strong class="jxr_keyword">double</strong> computeSqrt2aa(<strong class="jxr_keyword">double</strong> a) { +<a class="jxr_linenumber" name="L147" href="#L147">147</a> <em class="jxr_comment">// Split the number</em> +<a class="jxr_linenumber" name="L148" href="#L148">148</a> <strong class="jxr_keyword">final</strong> <strong class="jxr_keyword">double</strong> ha = highPartUnscaled(a); +<a class="jxr_linenumber" name="L149" href="#L149">149</a> <strong class="jxr_keyword">final</strong> <strong class="jxr_keyword">double</strong> la = a - ha; +<a class="jxr_linenumber" name="L150" href="#L150">150</a> +<a class="jxr_linenumber" name="L151" href="#L151">151</a> <em class="jxr_comment">// Extended precision product</em> +<a class="jxr_linenumber" name="L152" href="#L152">152</a> <strong class="jxr_keyword">final</strong> <strong class="jxr_keyword">double</strong> x = 2 * a * a; +<a class="jxr_linenumber" name="L153" href="#L153">153</a> <strong class="jxr_keyword">final</strong> <strong class="jxr_keyword">double</strong> xx = productLow(ha, la, 2 * ha, 2 * la, x); +<a class="jxr_linenumber" name="L154" href="#L154">154</a> +<a class="jxr_linenumber" name="L155" href="#L155">155</a> <em class="jxr_comment">// Standard sqrt</em> +<a class="jxr_linenumber" name="L156" href="#L156">156</a> <strong class="jxr_keyword">final</strong> <strong class="jxr_keyword">double</strong> c = Math.sqrt(x); +<a class="jxr_linenumber" name="L157" href="#L157">157</a> +<a class="jxr_linenumber" name="L158" href="#L158">158</a> <em class="jxr_comment">// Edge case.</em> +<a class="jxr_linenumber" name="L159" href="#L159">159</a> <em class="jxr_comment">// Occurs if a has limited precision in the mantissa including</em> +<a class="jxr_linenumber" name="L160" href="#L160">160</a> <em class="jxr_comment">// the special cases of 0 and 1.</em> +<a class="jxr_linenumber" name="L161" href="#L161">161</a> <strong class="jxr_keyword">if</strong> (xx == 0) { +<a class="jxr_linenumber" name="L162" href="#L162">162</a> <strong class="jxr_keyword">return</strong> c; +<a class="jxr_linenumber" name="L163" href="#L163">163</a> } +<a class="jxr_linenumber" name="L164" href="#L164">164</a> +<a class="jxr_linenumber" name="L165" href="#L165">165</a> <em class="jxr_comment">// Dekker's double precision sqrt2 algorithm.</em> +<a class="jxr_linenumber" name="L166" href="#L166">166</a> <em class="jxr_comment">// See Dekker, 1971, pp 242.</em> +<a class="jxr_linenumber" name="L167" href="#L167">167</a> <strong class="jxr_keyword">final</strong> <strong class="jxr_keyword">double</strong> hc = highPartUnscaled(c); +<a class="jxr_linenumber" name="L168" href="#L168">168</a> <strong class="jxr_keyword">final</strong> <strong class="jxr_keyword">double</strong> lc = c - hc; +<a class="jxr_linenumber" name="L169" href="#L169">169</a> <strong class="jxr_keyword">final</strong> <strong class="jxr_keyword">double</strong> u = c * c; +<a class="jxr_linenumber" name="L170" href="#L170">170</a> <strong class="jxr_keyword">final</strong> <strong class="jxr_keyword">double</strong> uu = productLow(hc, lc, hc, lc, u); +<a class="jxr_linenumber" name="L171" href="#L171">171</a> <strong class="jxr_keyword">final</strong> <strong class="jxr_keyword">double</strong> cc = (x - u - uu + xx) * 0.5 / c; +<a class="jxr_linenumber" name="L172" href="#L172">172</a> +<a class="jxr_linenumber" name="L173" href="#L173">173</a> <em class="jxr_comment">// Extended precision result:</em> +<a class="jxr_linenumber" name="L174" href="#L174">174</a> <em class="jxr_comment">// y = c + cc</em> +<a class="jxr_linenumber" name="L175" href="#L175">175</a> <em class="jxr_comment">// yy = c - y + cc</em> +<a class="jxr_linenumber" name="L176" href="#L176">176</a> <em class="jxr_comment">// Return only y</em> +<a class="jxr_linenumber" name="L177" href="#L177">177</a> <strong class="jxr_keyword">return</strong> c + cc; +<a class="jxr_linenumber" name="L178" href="#L178">178</a> } +<a class="jxr_linenumber" name="L179" href="#L179">179</a> +<a class="jxr_linenumber" name="L180" href="#L180">180</a> <em class="jxr_javadoccomment">/**</em> +<a class="jxr_linenumber" name="L181" href="#L181">181</a> <em class="jxr_javadoccomment"> * Compute {@code exp(-0.5*x*x)} with high accuracy. This is performed using information in the</em> +<a class="jxr_linenumber" name="L182" href="#L182">182</a> <em class="jxr_javadoccomment"> * round-off from {@code x*x}.</em> +<a class="jxr_linenumber" name="L183" href="#L183">183</a> <em class="jxr_javadoccomment"> *</em> +<a class="jxr_linenumber" name="L184" href="#L184">184</a> <em class="jxr_javadoccomment"> * <p>This is accurate at large x to 1 ulp until exp(-0.5*x*x) is close to sub-normal. For very</em> +<a class="jxr_linenumber" name="L185" href="#L185">185</a> <em class="jxr_javadoccomment"> * small exp(-0.5*x*x) the adjustment is sub-normal and bits can be lost in the adjustment for a</em> +<a class="jxr_linenumber" name="L186" href="#L186">186</a> <em class="jxr_javadoccomment"> * max observed error of {@code < 2} ulp.</em> +<a class="jxr_linenumber" name="L187" href="#L187">187</a> <em class="jxr_javadoccomment"> *</em> +<a class="jxr_linenumber" name="L188" href="#L188">188</a> <em class="jxr_javadoccomment"> * <p>At small x the accuracy cannot be improved over using exp(-0.5*x*x). This occurs at</em> +<a class="jxr_linenumber" name="L189" href="#L189">189</a> <em class="jxr_javadoccomment"> * {@code x <= sqrt(2)}.</em> +<a class="jxr_linenumber" name="L190" href="#L190">190</a> <em class="jxr_javadoccomment"> *</em> +<a class="jxr_linenumber" name="L191" href="#L191">191</a> <em class="jxr_javadoccomment"> * @param x Value</em> +<a class="jxr_linenumber" name="L192" href="#L192">192</a> <em class="jxr_javadoccomment"> * @return exp(-0.5*x*x)</em> +<a class="jxr_linenumber" name="L193" href="#L193">193</a> <em class="jxr_javadoccomment"> * @see <a href="<a href="https://issues.apache.org/jira/browse/STATISTICS-52"; target="alexandria_uri">https://issues.apache.org/jira/browse/STATISTICS-52</a>">STATISTICS-52</a></em> +<a class="jxr_linenumber" name="L194" href="#L194">194</a> <em class="jxr_javadoccomment"> */</em> +<a class="jxr_linenumber" name="L195" href="#L195">195</a> <strong class="jxr_keyword">static</strong> <strong class="jxr_keyword">double</strong> expmhxx(<strong class="jxr_keyword">double</strong> x) { +<a class="jxr_linenumber" name="L196" href="#L196">196</a> <strong class="jxr_keyword">final</strong> <strong class="jxr_keyword">double</strong> z = x * x; +<a class="jxr_linenumber" name="L197" href="#L197">197</a> <strong class="jxr_keyword">if</strong> (z <= EXP_M_HALF_XX_MIN_VALUE) { +<a class="jxr_linenumber" name="L198" href="#L198">198</a> <strong class="jxr_keyword">return</strong> Math.exp(-0.5 * z); +<a class="jxr_linenumber" name="L199" href="#L199">199</a> } <strong class="jxr_keyword">else</strong> <strong class="jxr_keyword">if</strong> (z >= EXP_M_HALF_XX_MAX_VALUE) { +<a class="jxr_linenumber" name="L200" href="#L200">200</a> <em class="jxr_comment">// exp(-745.5) == 0</em> +<a class="jxr_linenumber" name="L201" href="#L201">201</a> <strong class="jxr_keyword">return</strong> 0; +<a class="jxr_linenumber" name="L202" href="#L202">202</a> } +<a class="jxr_linenumber" name="L203" href="#L203">203</a> <em class="jxr_comment">// Split the number</em> +<a class="jxr_linenumber" name="L204" href="#L204">204</a> <strong class="jxr_keyword">final</strong> <strong class="jxr_keyword">double</strong> hx = highPartUnscaled(x); +<a class="jxr_linenumber" name="L205" href="#L205">205</a> <strong class="jxr_keyword">final</strong> <strong class="jxr_keyword">double</strong> lx = x - hx; +<a class="jxr_linenumber" name="L206" href="#L206">206</a> <em class="jxr_comment">// Compute the round-off</em> +<a class="jxr_linenumber" name="L207" href="#L207">207</a> <strong class="jxr_keyword">final</strong> <strong class="jxr_keyword">double</strong> zz = squareLow(hx, lx, z); +<a class="jxr_linenumber" name="L208" href="#L208">208</a> <strong class="jxr_keyword">return</strong> expxx(-0.5 * z, -0.5 * zz); +<a class="jxr_linenumber" name="L209" href="#L209">209</a> } +<a class="jxr_linenumber" name="L210" href="#L210">210</a> +<a class="jxr_linenumber" name="L211" href="#L211">211</a> <em class="jxr_javadoccomment">/**</em> +<a class="jxr_linenumber" name="L212" href="#L212">212</a> <em class="jxr_javadoccomment"> * Compute {@code exp(a+b)} with high accuracy assuming {@code a+b = a}.</em> +<a class="jxr_linenumber" name="L213" href="#L213">213</a> <em class="jxr_javadoccomment"> *</em> +<a class="jxr_linenumber" name="L214" href="#L214">214</a> <em class="jxr_javadoccomment"> * <p>This is accurate at large positive a to 1 ulp. If a is negative and exp(a) is close to</em> +<a class="jxr_linenumber" name="L215" href="#L215">215</a> <em class="jxr_javadoccomment"> * sub-normal a bit of precision may be lost when adjusting result as the adjustment is sub-normal</em> +<a class="jxr_linenumber" name="L216" href="#L216">216</a> <em class="jxr_javadoccomment"> * (max observed error {@code < 2} ulp). For the use case of multiplication of a number less than</em> +<a class="jxr_linenumber" name="L217" href="#L217">217</a> <em class="jxr_javadoccomment"> * 1 by exp(-x*x), a = -x*x, the result will be sub-normal and the rounding error is lost.</em> +<a class="jxr_linenumber" name="L218" href="#L218">218</a> <em class="jxr_javadoccomment"> *</em> +<a class="jxr_linenumber" name="L219" href="#L219">219</a> <em class="jxr_javadoccomment"> * <p>At small |a| the accuracy cannot be improved over using exp(a) as the round-off is too small</em> +<a class="jxr_linenumber" name="L220" href="#L220">220</a> <em class="jxr_javadoccomment"> * to create terms that can adjust the standard result by more than 0.5 ulp. This occurs at</em> +<a class="jxr_linenumber" name="L221" href="#L221">221</a> <em class="jxr_javadoccomment"> * {@code |a| <= 1}.</em> +<a class="jxr_linenumber" name="L222" href="#L222">222</a> <em class="jxr_javadoccomment"> *</em> +<a class="jxr_linenumber" name="L223" href="#L223">223</a> <em class="jxr_javadoccomment"> * @param a High bits of a split number</em> +<a class="jxr_linenumber" name="L224" href="#L224">224</a> <em class="jxr_javadoccomment"> * @param b Low bits of a split number</em> +<a class="jxr_linenumber" name="L225" href="#L225">225</a> <em class="jxr_javadoccomment"> * @return exp(a+b)</em> +<a class="jxr_linenumber" name="L226" href="#L226">226</a> <em class="jxr_javadoccomment"> * @see <a href="<a href="https://issues.apache.org/jira/projects/NUMBERS/issues/NUMBERS-177"; target="alexandria_uri">https://issues.apache.org/jira/projects/NUMBERS/issues/NUMBERS-177</a>"></em> +<a class="jxr_linenumber" name="L227" href="#L227">227</a> <em class="jxr_javadoccomment"> * Numbers-177: Accurate scaling by exp(z*z)</a></em> +<a class="jxr_linenumber" name="L228" href="#L228">228</a> <em class="jxr_javadoccomment"> */</em> +<a class="jxr_linenumber" name="L229" href="#L229">229</a> <strong class="jxr_keyword">private</strong> <strong class="jxr_keyword">static</strong> <strong class="jxr_keyword">double</strong> expxx(<strong class="jxr_keyword">double</strong> a, <strong class="jxr_keyword">double</strong> b) { +<a class="jxr_linenumber" name="L230" href="#L230">230</a> <em class="jxr_comment">// exp(a+b) = exp(a) * exp(b)</em> +<a class="jxr_linenumber" name="L231" href="#L231">231</a> <em class="jxr_comment">// = exp(a) * (exp(b) - 1) + exp(a)</em> +<a class="jxr_linenumber" name="L232" href="#L232">232</a> <em class="jxr_comment">// Assuming:</em> +<a class="jxr_linenumber" name="L233" href="#L233">233</a> <em class="jxr_comment">// 1. -746 < a < 710 for no under/overflow of exp(a)</em> +<a class="jxr_linenumber" name="L234" href="#L234">234</a> <em class="jxr_comment">// 2. a+b = a</em> +<a class="jxr_linenumber" name="L235" href="#L235">235</a> <em class="jxr_comment">// As b -> 0 then exp(b) -> 1; expm1(b) -> b</em> +<a class="jxr_linenumber" name="L236" href="#L236">236</a> <em class="jxr_comment">// The round-off b is limited to ~ 0.5 * ulp(746) ~ 5.68e-14</em> +<a class="jxr_linenumber" name="L237" href="#L237">237</a> <em class="jxr_comment">// and we can use an approximation for expm1 (x/1! + x^2/2! + ...)</em> +<a class="jxr_linenumber" name="L238" href="#L238">238</a> <em class="jxr_comment">// The second term is required for the expm1 result but the</em> +<a class="jxr_linenumber" name="L239" href="#L239">239</a> <em class="jxr_comment">// bits are not significant to change the product with exp(a)</em> +<a class="jxr_linenumber" name="L240" href="#L240">240</a> +<a class="jxr_linenumber" name="L241" href="#L241">241</a> <strong class="jxr_keyword">final</strong> <strong class="jxr_keyword">double</strong> ea = Math.exp(a); +<a class="jxr_linenumber" name="L242" href="#L242">242</a> <em class="jxr_comment">// b ~ expm1(b)</em> +<a class="jxr_linenumber" name="L243" href="#L243">243</a> <strong class="jxr_keyword">return</strong> ea * b + ea; +<a class="jxr_linenumber" name="L244" href="#L244">244</a> } +<a class="jxr_linenumber" name="L245" href="#L245">245</a> +<a class="jxr_linenumber" name="L246" href="#L246">246</a> <em class="jxr_javadoccomment">/**</em> +<a class="jxr_linenumber" name="L247" href="#L247">247</a> <em class="jxr_javadoccomment"> * Implement Dekker's method to split a value into two parts. Multiplying by (2^s + 1) creates</em> +<a class="jxr_linenumber" name="L248" href="#L248">248</a> <em class="jxr_javadoccomment"> * a big value from which to derive the two split parts.</em> +<a class="jxr_linenumber" name="L249" href="#L249">249</a> <em class="jxr_javadoccomment"> * <pre></em> +<a class="jxr_linenumber" name="L250" href="#L250">250</a> <em class="jxr_javadoccomment"> * c = (2^s + 1) * a</em> +<a class="jxr_linenumber" name="L251" href="#L251">251</a> <em class="jxr_javadoccomment"> * a_big = c - a</em> +<a class="jxr_linenumber" name="L252" href="#L252">252</a> <em class="jxr_javadoccomment"> * a_hi = c - a_big</em> +<a class="jxr_linenumber" name="L253" href="#L253">253</a> <em class="jxr_javadoccomment"> * a_lo = a - a_hi</em> +<a class="jxr_linenumber" name="L254" href="#L254">254</a> <em class="jxr_javadoccomment"> * a = a_hi + a_lo</em> +<a class="jxr_linenumber" name="L255" href="#L255">255</a> <em class="jxr_javadoccomment"> * </pre></em> +<a class="jxr_linenumber" name="L256" href="#L256">256</a> <em class="jxr_javadoccomment"> *</em> +<a class="jxr_linenumber" name="L257" href="#L257">257</a> <em class="jxr_javadoccomment"> * <p>The multiplicand allows a p-bit value to be split into</em> +<a class="jxr_linenumber" name="L258" href="#L258">258</a> <em class="jxr_javadoccomment"> * (p-s)-bit value {@code a_hi} and a non-overlapping (s-1)-bit value {@code a_lo}.</em> +<a class="jxr_linenumber" name="L259" href="#L259">259</a> <em class="jxr_javadoccomment"> * Combined they have (p-1) bits of significand but the sign bit of {@code a_lo}</em> +<a class="jxr_linenumber" name="L260" href="#L260">260</a> <em class="jxr_javadoccomment"> * contains a bit of information. The constant is chosen so that s is ceil(p/2) where</em> +<a class="jxr_linenumber" name="L261" href="#L261">261</a> <em class="jxr_javadoccomment"> * the precision p for a double is 53-bits (1-bit of the mantissa is assumed to be</em> +<a class="jxr_linenumber" name="L262" href="#L262">262</a> <em class="jxr_javadoccomment"> * 1 for a non sub-normal number) and s is 27.</em> +<a class="jxr_linenumber" name="L263" href="#L263">263</a> <em class="jxr_javadoccomment"> *</em> +<a class="jxr_linenumber" name="L264" href="#L264">264</a> <em class="jxr_javadoccomment"> * <p>This conversion does not use scaling and the result of overflow is NaN. Overflow</em> +<a class="jxr_linenumber" name="L265" href="#L265">265</a> <em class="jxr_javadoccomment"> * may occur when the exponent of the input value is above 996.</em> +<a class="jxr_linenumber" name="L266" href="#L266">266</a> <em class="jxr_javadoccomment"> *</em> +<a class="jxr_linenumber" name="L267" href="#L267">267</a> <em class="jxr_javadoccomment"> * <p>Splitting a NaN or infinite value will return NaN.</em> +<a class="jxr_linenumber" name="L268" href="#L268">268</a> <em class="jxr_javadoccomment"> *</em> +<a class="jxr_linenumber" name="L269" href="#L269">269</a> <em class="jxr_javadoccomment"> * @param value Value.</em> +<a class="jxr_linenumber" name="L270" href="#L270">270</a> <em class="jxr_javadoccomment"> * @return the high part of the value.</em> +<a class="jxr_linenumber" name="L271" href="#L271">271</a> <em class="jxr_javadoccomment"> * @see Math#getExponent(double)</em> +<a class="jxr_linenumber" name="L272" href="#L272">272</a> <em class="jxr_javadoccomment"> */</em> +<a class="jxr_linenumber" name="L273" href="#L273">273</a> <strong class="jxr_keyword">private</strong> <strong class="jxr_keyword">static</strong> <strong class="jxr_keyword">double</strong> highPartUnscaled(<strong class="jxr_keyword">double</strong> value) { +<a class="jxr_linenumber" name="L274" href="#L274">274</a> <strong class="jxr_keyword">final</strong> <strong class="jxr_keyword">double</strong> c = MULTIPLIER * value; +<a class="jxr_linenumber" name="L275" href="#L275">275</a> <strong class="jxr_keyword">return</strong> c - (c - value); +<a class="jxr_linenumber" name="L276" href="#L276">276</a> } +<a class="jxr_linenumber" name="L277" href="#L277">277</a> +<a class="jxr_linenumber" name="L278" href="#L278">278</a> <em class="jxr_javadoccomment">/**</em> +<a class="jxr_linenumber" name="L279" href="#L279">279</a> <em class="jxr_javadoccomment"> * Compute the low part of the double length number {@code (z,zz)} for the exact</em> +<a class="jxr_linenumber" name="L280" href="#L280">280</a> <em class="jxr_javadoccomment"> * product of {@code x} and {@code y} using Dekker's mult12 algorithm. The standard</em> +<a class="jxr_linenumber" name="L281" href="#L281">281</a> <em class="jxr_javadoccomment"> * precision product {@code x*y} must be provided. The numbers {@code x} and {@code y}</em> +<a class="jxr_linenumber" name="L282" href="#L282">282</a> <em class="jxr_javadoccomment"> * should already be split into low and high parts.</em> +<a class="jxr_linenumber" name="L283" href="#L283">283</a> <em class="jxr_javadoccomment"> *</em> +<a class="jxr_linenumber" name="L284" href="#L284">284</a> <em class="jxr_javadoccomment"> * <p>Note: This uses the high part of the result {@code (z,zz)} as {@code x * y} and not</em> +<a class="jxr_linenumber" name="L285" href="#L285">285</a> <em class="jxr_javadoccomment"> * {@code hx * hy + hx * ty + tx * hy} as specified in Dekker's original paper.</em> +<a class="jxr_linenumber" name="L286" href="#L286">286</a> <em class="jxr_javadoccomment"> * See Shewchuk (1997) for working examples.</em> +<a class="jxr_linenumber" name="L287" href="#L287">287</a> <em class="jxr_javadoccomment"> *</em> +<a class="jxr_linenumber" name="L288" href="#L288">288</a> <em class="jxr_javadoccomment"> * @param hx High part of first factor.</em> +<a class="jxr_linenumber" name="L289" href="#L289">289</a> <em class="jxr_javadoccomment"> * @param lx Low part of first factor.</em> +<a class="jxr_linenumber" name="L290" href="#L290">290</a> <em class="jxr_javadoccomment"> * @param hy High part of second factor.</em> +<a class="jxr_linenumber" name="L291" href="#L291">291</a> <em class="jxr_javadoccomment"> * @param ly Low part of second factor.</em> +<a class="jxr_linenumber" name="L292" href="#L292">292</a> <em class="jxr_javadoccomment"> * @param xy Product of the factors.</em> +<a class="jxr_linenumber" name="L293" href="#L293">293</a> <em class="jxr_javadoccomment"> * @return <code>lx * ly - (((xy - hx * hy) - lx * hy) - hx * ly)</code></em> +<a class="jxr_linenumber" name="L294" href="#L294">294</a> <em class="jxr_javadoccomment"> * @see <a href="<a href="https://www-2.cs.cmu.edu/afs/cs/project/quake/public/papers/robust-arithmetic.ps"; target="alexandria_uri">https://www-2.cs.cmu.edu/afs/cs/project/quake/public/papers/robust-arithmetic.ps</a>"></em> +<a class="jxr_linenumber" name="L295" href="#L295">295</a> <em class="jxr_javadoccomment"> * Shewchuk (1997) Theorum 18</a></em> +<a class="jxr_linenumber" name="L296" href="#L296">296</a> <em class="jxr_javadoccomment"> */</em> +<a class="jxr_linenumber" name="L297" href="#L297">297</a> <strong class="jxr_keyword">private</strong> <strong class="jxr_keyword">static</strong> <strong class="jxr_keyword">double</strong> productLow(<strong class="jxr_keyword">double</strong> hx, <strong class="jxr_keyword">double</strong> lx, <strong class="jxr_keyword">double</strong> hy, <strong class="jxr_keyword">double</strong> ly, <strong class="jxr_keyword">double</strong> xy) { +<a class="jxr_linenumber" name="L298" href="#L298">298</a> <em class="jxr_comment">// Compute the multiply low part:</em> +<a class="jxr_linenumber" name="L299" href="#L299">299</a> <em class="jxr_comment">// err1 = xy - hx * hy</em> +<a class="jxr_linenumber" name="L300" href="#L300">300</a> <em class="jxr_comment">// err2 = err1 - lx * hy</em> +<a class="jxr_linenumber" name="L301" href="#L301">301</a> <em class="jxr_comment">// err3 = err2 - hx * ly</em> +<a class="jxr_linenumber" name="L302" href="#L302">302</a> <em class="jxr_comment">// low = lx * ly - err3</em> +<a class="jxr_linenumber" name="L303" href="#L303">303</a> <strong class="jxr_keyword">return</strong> lx * ly - (((xy - hx * hy) - lx * hy) - hx * ly); +<a class="jxr_linenumber" name="L304" href="#L304">304</a> } +<a class="jxr_linenumber" name="L305" href="#L305">305</a> +<a class="jxr_linenumber" name="L306" href="#L306">306</a> <em class="jxr_javadoccomment">/**</em> +<a class="jxr_linenumber" name="L307" href="#L307">307</a> <em class="jxr_javadoccomment"> * Compute the low part of the double length number {@code (z,zz)} for the exact</em> +<a class="jxr_linenumber" name="L308" href="#L308">308</a> <em class="jxr_javadoccomment"> * square of {@code x} using Dekker's mult12 algorithm. The standard precision product</em> +<a class="jxr_linenumber" name="L309" href="#L309">309</a> <em class="jxr_javadoccomment"> * {@code x*x} must be provided. The number {@code x} should already be split into low</em> +<a class="jxr_linenumber" name="L310" href="#L310">310</a> <em class="jxr_javadoccomment"> * and high parts.</em> +<a class="jxr_linenumber" name="L311" href="#L311">311</a> <em class="jxr_javadoccomment"> *</em> +<a class="jxr_linenumber" name="L312" href="#L312">312</a> <em class="jxr_javadoccomment"> * <p>Note: This is a specialisation of</em> +<a class="jxr_linenumber" name="L313" href="#L313">313</a> <em class="jxr_javadoccomment"> * {@link #productLow(double, double, double, double, double)}.</em> +<a class="jxr_linenumber" name="L314" href="#L314">314</a> <em class="jxr_javadoccomment"> *</em> +<a class="jxr_linenumber" name="L315" href="#L315">315</a> <em class="jxr_javadoccomment"> * @param hx High part of factor.</em> +<a class="jxr_linenumber" name="L316" href="#L316">316</a> <em class="jxr_javadoccomment"> * @param lx Low part of factor.</em> +<a class="jxr_linenumber" name="L317" href="#L317">317</a> <em class="jxr_javadoccomment"> * @param xx Square of the factor.</em> +<a class="jxr_linenumber" name="L318" href="#L318">318</a> <em class="jxr_javadoccomment"> * @return <code>lx * lx - (((xx - hx * hx) - lx * hx) - hx * lx)</code></em> +<a class="jxr_linenumber" name="L319" href="#L319">319</a> <em class="jxr_javadoccomment"> */</em> +<a class="jxr_linenumber" name="L320" href="#L320">320</a> <strong class="jxr_keyword">private</strong> <strong class="jxr_keyword">static</strong> <strong class="jxr_keyword">double</strong> squareLow(<strong class="jxr_keyword">double</strong> hx, <strong class="jxr_keyword">double</strong> lx, <strong class="jxr_keyword">double</strong> xx) { +<a class="jxr_linenumber" name="L321" href="#L321">321</a> <strong class="jxr_keyword">return</strong> lx * lx - ((xx - hx * hx) - 2 * lx * hx); +<a class="jxr_linenumber" name="L322" href="#L322">322</a> } +<a class="jxr_linenumber" name="L323" href="#L323">323</a> } +</pre> +<hr/> +<div id="footer">Copyright © 2018–2022 <a href="https://www.apache.org/";>The Apache Software Foundation</a>. All rights reserved.</div> +</body> +</html>
