mirror of https://gitlab.com/rnger/amath
448 lines
50 KiB
HTML
448 lines
50 KiB
HTML
<!doctype html>
|
|
<html>
|
|
<head>
|
|
<title>lib/real/pow.c</title>
|
|
<style type="text/css">
|
|
body { color:#000000; background-color:#ffffff }
|
|
body { font-family:Helvetica, sans-serif; font-size:10pt }
|
|
h1 { font-size:14pt }
|
|
.code { border-collapse:collapse; width:100%; }
|
|
.code { font-family: "Monospace", monospace; font-size:10pt }
|
|
.code { line-height: 1.2em }
|
|
.comment { color: green; font-style: oblique }
|
|
.keyword { color: blue }
|
|
.string_literal { color: red }
|
|
.directive { color: darkmagenta }
|
|
.expansion { display: none; }
|
|
.macro:hover .expansion { display: block; border: 2px solid #FF0000; padding: 2px; background-color:#FFF0F0; font-weight: normal; -webkit-border-radius:5px; -webkit-box-shadow:1px 1px 7px #000; position: absolute; top: -1em; left:10em; z-index: 1 }
|
|
.macro { color: darkmagenta; background-color:LemonChiffon; position: relative }
|
|
.num { width:2.5em; padding-right:2ex; background-color:#eeeeee }
|
|
.num { text-align:right; font-size:8pt }
|
|
.num { color:#444444 }
|
|
.line { padding-left: 1ex; border-left: 3px solid #ccc }
|
|
.line { white-space: pre }
|
|
.msg { -webkit-box-shadow:1px 1px 7px #000 }
|
|
.msg { -webkit-border-radius:5px }
|
|
.msg { font-family:Helvetica, sans-serif; font-size:8pt }
|
|
.msg { float:left }
|
|
.msg { padding:0.25em 1ex 0.25em 1ex }
|
|
.msg { margin-top:10px; margin-bottom:10px }
|
|
.msg { font-weight:bold }
|
|
.msg { max-width:60em; word-wrap: break-word; white-space: pre-wrap }
|
|
.msgT { padding:0x; spacing:0x }
|
|
.msgEvent { background-color:#fff8b4; color:#000000 }
|
|
.msgControl { background-color:#bbbbbb; color:#000000 }
|
|
.mrange { background-color:#dfddf3 }
|
|
.mrange { border-bottom:1px solid #6F9DBE }
|
|
.PathIndex { font-weight: bold; padding:0px 5px; margin-right:5px; }
|
|
.PathIndex { -webkit-border-radius:8px }
|
|
.PathIndexEvent { background-color:#bfba87 }
|
|
.PathIndexControl { background-color:#8c8c8c }
|
|
.PathNav a { text-decoration:none; font-size: larger }
|
|
.CodeInsertionHint { font-weight: bold; background-color: #10dd10 }
|
|
.CodeRemovalHint { background-color:#de1010 }
|
|
.CodeRemovalHint { border-bottom:1px solid #6F9DBE }
|
|
table.simpletable {
|
|
padding: 5px;
|
|
font-size:12pt;
|
|
margin:20px;
|
|
border-collapse: collapse; border-spacing: 0px;
|
|
}
|
|
td.rowname {
|
|
text-align:right; font-weight:bold; color:#444444;
|
|
padding-right:2ex; }
|
|
</style>
|
|
</head>
|
|
<body>
|
|
<!-- BUGDESC Value stored to 'i1' is never read -->
|
|
|
|
<!-- BUGTYPE Dead assignment -->
|
|
|
|
<!-- BUGCATEGORY Dead store -->
|
|
|
|
<!-- BUGFILE /home/carsten/amath/lib/real/pow.c -->
|
|
|
|
<!-- FILENAME pow.c -->
|
|
|
|
<!-- FUNCTIONNAME pow -->
|
|
|
|
<!-- ISSUEHASHCONTENTOFLINEINCONTEXT e1b4b10697e719e7dd943d54d2e75047 -->
|
|
|
|
<!-- BUGLINE 148 -->
|
|
|
|
<!-- BUGCOLUMN 5 -->
|
|
|
|
<!-- BUGPATHLENGTH 1 -->
|
|
|
|
<!-- BUGMETAEND -->
|
|
<!-- REPORTHEADER -->
|
|
<h3>Bug Summary</h3>
|
|
<table class="simpletable">
|
|
<tr><td class="rowname">File:</td><td>real/pow.c</td></tr>
|
|
<tr><td class="rowname">Location:</td><td><a href="#EndPath">line 148, column 5</a></td></tr>
|
|
<tr><td class="rowname">Description:</td><td>Value stored to 'i1' is never read</td></tr>
|
|
</table>
|
|
<!-- REPORTSUMMARYEXTRA -->
|
|
<h3>Annotated Source Code</h3>
|
|
<table class="code">
|
|
<tr><td class="num" id="LN1">1</td><td class="line"><span class='comment'>/* @(#)e_pow.c 1.5 04/04/22 SMI */</span></td></tr>
|
|
<tr><td class="num" id="LN2">2</td><td class="line"> </td></tr>
|
|
<tr><td class="num" id="LN3">3</td><td class="line"><span class='comment'>/*</span></td></tr>
|
|
<tr><td class="num" id="LN4">4</td><td class="line"> <span class='comment'>* Copyright (c) 2015-2017 Carsten Sonne Larsen <cs@innolan.dk></span></td></tr>
|
|
<tr><td class="num" id="LN5">5</td><td class="line"> <span class='comment'>* All rights reserved.</span></td></tr>
|
|
<tr><td class="num" id="LN6">6</td><td class="line"> <span class='comment'>*</span></td></tr>
|
|
<tr><td class="num" id="LN7">7</td><td class="line"> <span class='comment'>* Redistribution and use in source and binary forms, with or without</span></td></tr>
|
|
<tr><td class="num" id="LN8">8</td><td class="line"> <span class='comment'>* modification, are permitted provided that the following conditions</span></td></tr>
|
|
<tr><td class="num" id="LN9">9</td><td class="line"> <span class='comment'>* are met:</span></td></tr>
|
|
<tr><td class="num" id="LN10">10</td><td class="line"> <span class='comment'>* 1. Redistributions of source code must retain the above copyright</span></td></tr>
|
|
<tr><td class="num" id="LN11">11</td><td class="line"> <span class='comment'>* notice, this list of conditions and the following disclaimer.</span></td></tr>
|
|
<tr><td class="num" id="LN12">12</td><td class="line"> <span class='comment'>* 2. Redistributions in binary form must reproduce the above copyright</span></td></tr>
|
|
<tr><td class="num" id="LN13">13</td><td class="line"> <span class='comment'>* notice, this list of conditions and the following disclaimer in the</span></td></tr>
|
|
<tr><td class="num" id="LN14">14</td><td class="line"> <span class='comment'>* documentation and/or other materials provided with the distribution.</span></td></tr>
|
|
<tr><td class="num" id="LN15">15</td><td class="line"> <span class='comment'>*</span></td></tr>
|
|
<tr><td class="num" id="LN16">16</td><td class="line"> <span class='comment'>* THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR</span></td></tr>
|
|
<tr><td class="num" id="LN17">17</td><td class="line"> <span class='comment'>* IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES</span></td></tr>
|
|
<tr><td class="num" id="LN18">18</td><td class="line"> <span class='comment'>* OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.</span></td></tr>
|
|
<tr><td class="num" id="LN19">19</td><td class="line"> <span class='comment'>* IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,</span></td></tr>
|
|
<tr><td class="num" id="LN20">20</td><td class="line"> <span class='comment'>* INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT</span></td></tr>
|
|
<tr><td class="num" id="LN21">21</td><td class="line"> <span class='comment'>* NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,</span></td></tr>
|
|
<tr><td class="num" id="LN22">22</td><td class="line"> <span class='comment'>* DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY</span></td></tr>
|
|
<tr><td class="num" id="LN23">23</td><td class="line"> <span class='comment'>* THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT</span></td></tr>
|
|
<tr><td class="num" id="LN24">24</td><td class="line"> <span class='comment'>* (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF</span></td></tr>
|
|
<tr><td class="num" id="LN25">25</td><td class="line"> <span class='comment'>* THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.</span></td></tr>
|
|
<tr><td class="num" id="LN26">26</td><td class="line"> <span class='comment'>*</span></td></tr>
|
|
<tr><td class="num" id="LN27">27</td><td class="line"> <span class='comment'>* The origin source code can be obtained from:</span></td></tr>
|
|
<tr><td class="num" id="LN28">28</td><td class="line"> <span class='comment'>* http://www.netlib.org/fdlibm/e_pow.c</span></td></tr>
|
|
<tr><td class="num" id="LN29">29</td><td class="line"> <span class='comment'>*</span> </td></tr>
|
|
<tr><td class="num" id="LN30">30</td><td class="line"> <span class='comment'>*/</span></td></tr>
|
|
<tr><td class="num" id="LN31">31</td><td class="line"> </td></tr>
|
|
<tr><td class="num" id="LN32">32</td><td class="line"><span class='comment'>/*</span></td></tr>
|
|
<tr><td class="num" id="LN33">33</td><td class="line"> <span class='comment'>* ====================================================</span></td></tr>
|
|
<tr><td class="num" id="LN34">34</td><td class="line"> <span class='comment'>* Copyright (C) 2004 by Sun Microsystems, Inc. All rights reserved.</span></td></tr>
|
|
<tr><td class="num" id="LN35">35</td><td class="line"> <span class='comment'>*</span></td></tr>
|
|
<tr><td class="num" id="LN36">36</td><td class="line"> <span class='comment'>* Permission to use, copy, modify, and distribute this</span></td></tr>
|
|
<tr><td class="num" id="LN37">37</td><td class="line"> <span class='comment'>* software is freely granted, provided that this notice</span></td></tr>
|
|
<tr><td class="num" id="LN38">38</td><td class="line"> <span class='comment'>* is preserved.</span></td></tr>
|
|
<tr><td class="num" id="LN39">39</td><td class="line"> <span class='comment'>* ====================================================</span></td></tr>
|
|
<tr><td class="num" id="LN40">40</td><td class="line"> <span class='comment'>*/</span></td></tr>
|
|
<tr><td class="num" id="LN41">41</td><td class="line"> </td></tr>
|
|
<tr><td class="num" id="LN42">42</td><td class="line"><span class='directive'>#include "prim.h"</span></td></tr>
|
|
<tr><td class="num" id="LN43">43</td><td class="line"><span class='directive'>#include "math.h"</span></td></tr>
|
|
<tr><td class="num" id="LN44">44</td><td class="line"> </td></tr>
|
|
<tr><td class="num" id="LN45">45</td><td class="line"><span class='directive'>#ifdef <span class='macro'>__clang__<span class='expansion'>1</span></span></span></td></tr>
|
|
<tr><td class="num" id="LN46">46</td><td class="line"><span class='directive'># pragma clang diagnostic ignored "-Wunused-variable"</span></td></tr>
|
|
<tr><td class="num" id="LN47">47</td><td class="line"><span class='directive'># pragma clang diagnostic ignored "-Wstrict-aliasing"</span></td></tr>
|
|
<tr><td class="num" id="LN48">48</td><td class="line"><span class='directive'>#elif defined(<span class='macro'>__GNUC__<span class='expansion'>4</span></span>) && <span class='macro'>__GNUC__<span class='expansion'>4</span></span> > 2</span></td></tr>
|
|
<tr><td class="num" id="LN49">49</td><td class="line"><span class='directive'># pragma GCC diagnostic ignored "-Wunused-but-set-variable"</span></td></tr>
|
|
<tr><td class="num" id="LN50">50</td><td class="line"><span class='directive'># pragma GCC diagnostic ignored "-Wstrict-aliasing"</span></td></tr>
|
|
<tr><td class="num" id="LN51">51</td><td class="line"><span class='directive'>#endif</span></td></tr>
|
|
<tr><td class="num" id="LN52">52</td><td class="line"> </td></tr>
|
|
<tr><td class="num" id="LN53">53</td><td class="line"><span class='keyword'>static</span> <span class='keyword'>const</span> <span class='keyword'>double</span></td></tr>
|
|
<tr><td class="num" id="LN54">54</td><td class="line">bp[] = {1.0, 1.5,},</td></tr>
|
|
<tr><td class="num" id="LN55">55</td><td class="line">dp_h[] = { 0.0, 5.84962487220764160156e-01,}, <span class='comment'>/* 0x3FE2B803, 0x40000000 */</span></td></tr>
|
|
<tr><td class="num" id="LN56">56</td><td class="line">dp_l[] = { 0.0, 1.35003920212974897128e-08,}, <span class='comment'>/* 0x3E4CFDEB, 0x43CFD006 */</span></td></tr>
|
|
<tr><td class="num" id="LN57">57</td><td class="line">zero = 0.0,</td></tr>
|
|
<tr><td class="num" id="LN58">58</td><td class="line">one = 1.0,</td></tr>
|
|
<tr><td class="num" id="LN59">59</td><td class="line">two = 2.0,</td></tr>
|
|
<tr><td class="num" id="LN60">60</td><td class="line">two53 = 9007199254740992.0, <span class='comment'>/* 0x43400000, 0x00000000 */</span></td></tr>
|
|
<tr><td class="num" id="LN61">61</td><td class="line">huge = 1.0e300,</td></tr>
|
|
<tr><td class="num" id="LN62">62</td><td class="line">tiny = 1.0e-300,</td></tr>
|
|
<tr><td class="num" id="LN63">63</td><td class="line"> <span class='comment'>/* poly coefs for (3/2)*(log(x)-2s-2/3*s**3 */</span></td></tr>
|
|
<tr><td class="num" id="LN64">64</td><td class="line">L1 = 5.99999999999994648725e-01, <span class='comment'>/* 0x3FE33333, 0x33333303 */</span></td></tr>
|
|
<tr><td class="num" id="LN65">65</td><td class="line">L2 = 4.28571428578550184252e-01, <span class='comment'>/* 0x3FDB6DB6, 0xDB6FABFF */</span></td></tr>
|
|
<tr><td class="num" id="LN66">66</td><td class="line">L3 = 3.33333329818377432918e-01, <span class='comment'>/* 0x3FD55555, 0x518F264D */</span></td></tr>
|
|
<tr><td class="num" id="LN67">67</td><td class="line">L4 = 2.72728123808534006489e-01, <span class='comment'>/* 0x3FD17460, 0xA91D4101 */</span></td></tr>
|
|
<tr><td class="num" id="LN68">68</td><td class="line">L5 = 2.30660745775561754067e-01, <span class='comment'>/* 0x3FCD864A, 0x93C9DB65 */</span></td></tr>
|
|
<tr><td class="num" id="LN69">69</td><td class="line">L6 = 2.06975017800338417784e-01, <span class='comment'>/* 0x3FCA7E28, 0x4A454EEF */</span></td></tr>
|
|
<tr><td class="num" id="LN70">70</td><td class="line">P1 = 1.66666666666666019037e-01, <span class='comment'>/* 0x3FC55555, 0x5555553E */</span></td></tr>
|
|
<tr><td class="num" id="LN71">71</td><td class="line">P2 = -2.77777777770155933842e-03, <span class='comment'>/* 0xBF66C16C, 0x16BEBD93 */</span></td></tr>
|
|
<tr><td class="num" id="LN72">72</td><td class="line">P3 = 6.61375632143793436117e-05, <span class='comment'>/* 0x3F11566A, 0xAF25DE2C */</span></td></tr>
|
|
<tr><td class="num" id="LN73">73</td><td class="line">P4 = -1.65339022054652515390e-06, <span class='comment'>/* 0xBEBBBD41, 0xC5D26BF1 */</span></td></tr>
|
|
<tr><td class="num" id="LN74">74</td><td class="line">P5 = 4.13813679705723846039e-08, <span class='comment'>/* 0x3E663769, 0x72BEA4D0 */</span></td></tr>
|
|
<tr><td class="num" id="LN75">75</td><td class="line">lg2 = 6.93147180559945286227e-01, <span class='comment'>/* 0x3FE62E42, 0xFEFA39EF */</span></td></tr>
|
|
<tr><td class="num" id="LN76">76</td><td class="line">lg2_h = 6.93147182464599609375e-01, <span class='comment'>/* 0x3FE62E43, 0x00000000 */</span></td></tr>
|
|
<tr><td class="num" id="LN77">77</td><td class="line">lg2_l = -1.90465429995776804525e-09, <span class='comment'>/* 0xBE205C61, 0x0CA86C39 */</span></td></tr>
|
|
<tr><td class="num" id="LN78">78</td><td class="line">ovt = 8.0085662595372944372e-0017, <span class='comment'>/* -(1024-log2(ovfl+.5ulp)) */</span></td></tr>
|
|
<tr><td class="num" id="LN79">79</td><td class="line">cp = 9.61796693925975554329e-01, <span class='comment'>/* 0x3FEEC709, 0xDC3A03FD =2/(3ln2) */</span></td></tr>
|
|
<tr><td class="num" id="LN80">80</td><td class="line">cp_h = 9.61796700954437255859e-01, <span class='comment'>/* 0x3FEEC709, 0xE0000000 =(float)cp */</span></td></tr>
|
|
<tr><td class="num" id="LN81">81</td><td class="line">cp_l = -7.02846165095275826516e-09, <span class='comment'>/* 0xBE3E2FE0, 0x145B01F5 =tail of cp_h*/</span></td></tr>
|
|
<tr><td class="num" id="LN82">82</td><td class="line">ivln2 = 1.44269504088896338700e+00, <span class='comment'>/* 0x3FF71547, 0x652B82FE =1/ln2 */</span></td></tr>
|
|
<tr><td class="num" id="LN83">83</td><td class="line">ivln2_h = 1.44269502162933349609e+00, <span class='comment'>/* 0x3FF71547, 0x60000000 =24b 1/ln2*/</span></td></tr>
|
|
<tr><td class="num" id="LN84">84</td><td class="line">ivln2_l = 1.92596299112661746887e-08; <span class='comment'>/* 0x3E54AE0B, 0xF85DDF44 =1/ln2 tail*/</span></td></tr>
|
|
<tr><td class="num" id="LN85">85</td><td class="line"> </td></tr>
|
|
<tr><td class="num" id="LN86">86</td><td class="line"><span class='comment'>/**</span></td></tr>
|
|
<tr><td class="num" id="LN87">87</td><td class="line"> <span class='comment'>* @brief Expontation function.</span></td></tr>
|
|
<tr><td class="num" id="LN88">88</td><td class="line"> <span class='comment'>* @version 1.3</span></td></tr>
|
|
<tr><td class="num" id="LN89">89</td><td class="line"> <span class='comment'>* @date 95/01/18</span></td></tr>
|
|
<tr><td class="num" id="LN90">90</td><td class="line"> <span class='comment'>* @details</span></td></tr>
|
|
<tr><td class="num" id="LN91">91</td><td class="line"> <span class='comment'>* <pre></span></td></tr>
|
|
<tr><td class="num" id="LN92">92</td><td class="line"> <span class='comment'>* Method: Let x = 2 * (1+f)</span></td></tr>
|
|
<tr><td class="num" id="LN93">93</td><td class="line"> <span class='comment'>* 1. Compute and return log2(x) in two pieces:</span></td></tr>
|
|
<tr><td class="num" id="LN94">94</td><td class="line"> <span class='comment'>* log2(x) = w1 + w2,</span></td></tr>
|
|
<tr><td class="num" id="LN95">95</td><td class="line"> <span class='comment'>* where w1 has 53-24 = 29 bit trailing zeros.</span></td></tr>
|
|
<tr><td class="num" id="LN96">96</td><td class="line"> <span class='comment'>* 2. Perform y*log2(x) = n+y' by simulating muti-precision</span></td></tr>
|
|
<tr><td class="num" id="LN97">97</td><td class="line"> <span class='comment'>* arithmetic, where |y'|<=0.5.</span></td></tr>
|
|
<tr><td class="num" id="LN98">98</td><td class="line"> <span class='comment'>* 3. Return x**y = 2**n*exp(y'*log2)</span></td></tr>
|
|
<tr><td class="num" id="LN99">99</td><td class="line"> <span class='comment'>*</span></td></tr>
|
|
<tr><td class="num" id="LN100">100</td><td class="line"> <span class='comment'>* Special cases:</span></td></tr>
|
|
<tr><td class="num" id="LN101">101</td><td class="line"> <span class='comment'>* 1. (anything) ** 0 is 1</span></td></tr>
|
|
<tr><td class="num" id="LN102">102</td><td class="line"> <span class='comment'>* 2. (anything) ** 1 is itself</span></td></tr>
|
|
<tr><td class="num" id="LN103">103</td><td class="line"> <span class='comment'>* 3. (anything) ** NAN is NAN</span></td></tr>
|
|
<tr><td class="num" id="LN104">104</td><td class="line"> <span class='comment'>* 4. NAN ** (anything except 0) is NAN</span></td></tr>
|
|
<tr><td class="num" id="LN105">105</td><td class="line"> <span class='comment'>* 5. +-(|x| > 1) ** +INF is +INF</span></td></tr>
|
|
<tr><td class="num" id="LN106">106</td><td class="line"> <span class='comment'>* 6. +-(|x| > 1) ** -INF is +0</span></td></tr>
|
|
<tr><td class="num" id="LN107">107</td><td class="line"> <span class='comment'>* 7. +-(|x| < 1) ** +INF is +0</span></td></tr>
|
|
<tr><td class="num" id="LN108">108</td><td class="line"> <span class='comment'>* 8. +-(|x| < 1) ** -INF is +INF</span></td></tr>
|
|
<tr><td class="num" id="LN109">109</td><td class="line"> <span class='comment'>* 9. +-1 ** +-INF is NAN</span></td></tr>
|
|
<tr><td class="num" id="LN110">110</td><td class="line"> <span class='comment'>* 10. +0 ** (+anything except 0, NAN) is +0</span></td></tr>
|
|
<tr><td class="num" id="LN111">111</td><td class="line"> <span class='comment'>* 11. -0 ** (+anything except 0, NAN, odd integer) is +0</span></td></tr>
|
|
<tr><td class="num" id="LN112">112</td><td class="line"> <span class='comment'>* 12. +0 ** (-anything except 0, NAN) is +INF</span></td></tr>
|
|
<tr><td class="num" id="LN113">113</td><td class="line"> <span class='comment'>* 13. -0 ** (-anything except 0, NAN, odd integer) is +INF</span></td></tr>
|
|
<tr><td class="num" id="LN114">114</td><td class="line"> <span class='comment'>* 14. -0 ** (odd integer) = -( +0 ** (odd integer) )</span></td></tr>
|
|
<tr><td class="num" id="LN115">115</td><td class="line"> <span class='comment'>* 15. +INF ** (+anything except 0,NAN) is +INF</span></td></tr>
|
|
<tr><td class="num" id="LN116">116</td><td class="line"> <span class='comment'>* 16. +INF ** (-anything except 0,NAN) is +0</span></td></tr>
|
|
<tr><td class="num" id="LN117">117</td><td class="line"> <span class='comment'>* 17. -INF ** (anything) = -0 ** (-anything)</span></td></tr>
|
|
<tr><td class="num" id="LN118">118</td><td class="line"> <span class='comment'>* 18. (-anything) ** (integer) is (-1)**(integer)*(+anything**integer)</span></td></tr>
|
|
<tr><td class="num" id="LN119">119</td><td class="line"> <span class='comment'>* 19. (-anything except 0 and inf) ** (non-integer) is NAN</span></td></tr>
|
|
<tr><td class="num" id="LN120">120</td><td class="line"> <span class='comment'>*</span></td></tr>
|
|
<tr><td class="num" id="LN121">121</td><td class="line"> <span class='comment'>* Accuracy:</span></td></tr>
|
|
<tr><td class="num" id="LN122">122</td><td class="line"> <span class='comment'>* pow(x,y) returns x**y nearly rounded. In particular</span></td></tr>
|
|
<tr><td class="num" id="LN123">123</td><td class="line"> <span class='comment'>* pow(integer,integer)</span></td></tr>
|
|
<tr><td class="num" id="LN124">124</td><td class="line"> <span class='comment'>* always returns the correct integer provided it is</span></td></tr>
|
|
<tr><td class="num" id="LN125">125</td><td class="line"> <span class='comment'>* representable.</span></td></tr>
|
|
<tr><td class="num" id="LN126">126</td><td class="line"> <span class='comment'>*</span></td></tr>
|
|
<tr><td class="num" id="LN127">127</td><td class="line"> <span class='comment'>* Constants :</span></td></tr>
|
|
<tr><td class="num" id="LN128">128</td><td class="line"> <span class='comment'>* The hexadecimal values are the intended ones for the following</span></td></tr>
|
|
<tr><td class="num" id="LN129">129</td><td class="line"> <span class='comment'>* constants. The decimal values may be used, provided that the</span></td></tr>
|
|
<tr><td class="num" id="LN130">130</td><td class="line"> <span class='comment'>* compiler will convert from decimal to binary accurately enough</span></td></tr>
|
|
<tr><td class="num" id="LN131">131</td><td class="line"> <span class='comment'>* to produce the hexadecimal values shown.</span></td></tr>
|
|
<tr><td class="num" id="LN132">132</td><td class="line"> <span class='comment'>* </pre></span></td></tr>
|
|
<tr><td class="num" id="LN133">133</td><td class="line"> <span class='comment'>* @copyright Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.</span></td></tr>
|
|
<tr><td class="num" id="LN134">134</td><td class="line"> <span class='comment'>* @license Developed at SunSoft, a Sun Microsystems, Inc. business. Permission</span></td></tr>
|
|
<tr><td class="num" id="LN135">135</td><td class="line"> <span class='comment'>* to use, copy, modify, and distribute this software is freely granted,</span></td></tr>
|
|
<tr><td class="num" id="LN136">136</td><td class="line"> <span class='comment'>* provided that this notice is preserved.</span></td></tr>
|
|
<tr><td class="num" id="LN137">137</td><td class="line"> <span class='comment'>*/</span></td></tr>
|
|
<tr><td class="num" id="LN138">138</td><td class="line"> </td></tr>
|
|
<tr><td class="num" id="LN139">139</td><td class="line"><span class='keyword'>double</span> pow(<span class='keyword'>double</span> x, <span class='keyword'>double</span> y)</td></tr>
|
|
<tr><td class="num" id="LN140">140</td><td class="line">{</td></tr>
|
|
<tr><td class="num" id="LN141">141</td><td class="line"> <span class='keyword'>double</span> z,ax,z_h,z_l,p_h,p_l;</td></tr>
|
|
<tr><td class="num" id="LN142">142</td><td class="line"> <span class='keyword'>double</span> y1,t1,t2,r,s,t,u,v,w;</td></tr>
|
|
<tr><td class="num" id="LN143">143</td><td class="line"> sword i0,i1,i,j,k,yisint,n;</td></tr>
|
|
<tr><td class="num" id="LN144">144</td><td class="line"> sword hx,hy,ix,iy;</td></tr>
|
|
<tr><td class="num" id="LN145">145</td><td class="line"> uword lx,ly;</td></tr>
|
|
<tr><td class="num" id="LN146">146</td><td class="line"> </td></tr>
|
|
<tr><td class="num" id="LN147">147</td><td class="line"> i0 = ((*(<span class='keyword'>int</span>*)&one)>>29)^1;</td></tr>
|
|
<tr><td class="num" id="LN148">148</td><td class="line"> i1=<span class="mrange">1-i0</span>;</td></tr>
|
|
<tr><td class="num"></td><td class="line"><div id="EndPath" class="msg msgEvent" style="margin-left:5ex">Value stored to 'i1' is never read</div></td></tr>
|
|
<tr><td class="num" id="LN149">149</td><td class="line"> <span class='macro'>EXTRACT_WORDS(hx,lx,x)<span class='expansion'>do { ieee_double_shape_type ew_u; ew_u.value = (x); (hx) = ew_u<br>.parts.msw; (lx) = ew_u.parts.lsw; } while (0)</span></span>;</td></tr>
|
|
<tr><td class="num" id="LN150">150</td><td class="line"> <span class='macro'>EXTRACT_WORDS(hy,ly,y)<span class='expansion'>do { ieee_double_shape_type ew_u; ew_u.value = (y); (hy) = ew_u<br>.parts.msw; (ly) = ew_u.parts.lsw; } while (0)</span></span>;</td></tr>
|
|
<tr><td class="num" id="LN151">151</td><td class="line"> ix = hx&0x7fffffff;</td></tr>
|
|
<tr><td class="num" id="LN152">152</td><td class="line"> iy = hy&0x7fffffff;</td></tr>
|
|
<tr><td class="num" id="LN153">153</td><td class="line"> </td></tr>
|
|
<tr><td class="num" id="LN154">154</td><td class="line"> <span class='comment'>/* y==zero: x**0 = 1 */</span></td></tr>
|
|
<tr><td class="num" id="LN155">155</td><td class="line"> <span class='keyword'>if</span>((iy|ly)==0) <span class='keyword'>return</span> one;</td></tr>
|
|
<tr><td class="num" id="LN156">156</td><td class="line"> </td></tr>
|
|
<tr><td class="num" id="LN157">157</td><td class="line"> <span class='comment'>/* +-NaN return x+y */</span></td></tr>
|
|
<tr><td class="num" id="LN158">158</td><td class="line"> <span class='keyword'>if</span>(ix > 0x7ff00000 || ((ix==0x7ff00000)&&(lx!=0)) ||</td></tr>
|
|
<tr><td class="num" id="LN159">159</td><td class="line"> iy > 0x7ff00000 || ((iy==0x7ff00000)&&(ly!=0)))</td></tr>
|
|
<tr><td class="num" id="LN160">160</td><td class="line"> <span class='keyword'>return</span> x+y;</td></tr>
|
|
<tr><td class="num" id="LN161">161</td><td class="line"> </td></tr>
|
|
<tr><td class="num" id="LN162">162</td><td class="line"> <span class='comment'>/* determine if y is an odd int when x < 0</span></td></tr>
|
|
<tr><td class="num" id="LN163">163</td><td class="line"> <span class='comment'>* yisint = 0 ... y is not an integer</span></td></tr>
|
|
<tr><td class="num" id="LN164">164</td><td class="line"> <span class='comment'>* yisint = 1 ... y is an odd int</span></td></tr>
|
|
<tr><td class="num" id="LN165">165</td><td class="line"> <span class='comment'>* yisint = 2 ... y is an even int</span></td></tr>
|
|
<tr><td class="num" id="LN166">166</td><td class="line"> <span class='comment'>*/</span></td></tr>
|
|
<tr><td class="num" id="LN167">167</td><td class="line"> yisint = 0;</td></tr>
|
|
<tr><td class="num" id="LN168">168</td><td class="line"> <span class='keyword'>if</span>(hx<0) {</td></tr>
|
|
<tr><td class="num" id="LN169">169</td><td class="line"> <span class='keyword'>if</span>(iy>=0x43400000) yisint = 2; <span class='comment'>/* even integer y */</span></td></tr>
|
|
<tr><td class="num" id="LN170">170</td><td class="line"> <span class='keyword'>else</span> <span class='keyword'>if</span>(iy>=0x3ff00000) {</td></tr>
|
|
<tr><td class="num" id="LN171">171</td><td class="line"> k = (iy>>20)-0x3ff; <span class='comment'>/* exponent */</span></td></tr>
|
|
<tr><td class="num" id="LN172">172</td><td class="line"> <span class='keyword'>if</span>(k>20) {</td></tr>
|
|
<tr><td class="num" id="LN173">173</td><td class="line"> j = ly>>(52-k);</td></tr>
|
|
<tr><td class="num" id="LN174">174</td><td class="line"> <span class='keyword'>if</span>((uword)(j<<(52-k))==ly) yisint = 2-(j&1);</td></tr>
|
|
<tr><td class="num" id="LN175">175</td><td class="line"> } <span class='keyword'>else</span> <span class='keyword'>if</span>(ly==0) {</td></tr>
|
|
<tr><td class="num" id="LN176">176</td><td class="line"> j = iy>>(20-k);</td></tr>
|
|
<tr><td class="num" id="LN177">177</td><td class="line"> <span class='keyword'>if</span>((j<<(20-k))==iy) yisint = 2-(j&1);</td></tr>
|
|
<tr><td class="num" id="LN178">178</td><td class="line"> }</td></tr>
|
|
<tr><td class="num" id="LN179">179</td><td class="line"> }</td></tr>
|
|
<tr><td class="num" id="LN180">180</td><td class="line"> }</td></tr>
|
|
<tr><td class="num" id="LN181">181</td><td class="line"> </td></tr>
|
|
<tr><td class="num" id="LN182">182</td><td class="line"> <span class='comment'>/* special value of y */</span></td></tr>
|
|
<tr><td class="num" id="LN183">183</td><td class="line"> <span class='keyword'>if</span>(ly==0) {</td></tr>
|
|
<tr><td class="num" id="LN184">184</td><td class="line"> <span class='keyword'>if</span> (iy==0x7ff00000) { <span class='comment'>/* y is +-inf */</span></td></tr>
|
|
<tr><td class="num" id="LN185">185</td><td class="line"> <span class='keyword'>if</span>(((ix-0x3ff00000)|lx)==0)</td></tr>
|
|
<tr><td class="num" id="LN186">186</td><td class="line"> <span class='keyword'>return</span> y - y; <span class='comment'>/* inf**+-1 is NaN */</span></td></tr>
|
|
<tr><td class="num" id="LN187">187</td><td class="line"> <span class='keyword'>else</span> <span class='keyword'>if</span> (ix >= 0x3ff00000)<span class='comment'>/* (|x|>1)**+-inf = inf,0 */</span></td></tr>
|
|
<tr><td class="num" id="LN188">188</td><td class="line"> <span class='keyword'>return</span> (hy>=0)? y: zero;</td></tr>
|
|
<tr><td class="num" id="LN189">189</td><td class="line"> <span class='keyword'>else</span> <span class='comment'>/* (|x|<1)**-,+inf = inf,0 */</span></td></tr>
|
|
<tr><td class="num" id="LN190">190</td><td class="line"> <span class='keyword'>return</span> (hy<0)?-y: zero;</td></tr>
|
|
<tr><td class="num" id="LN191">191</td><td class="line"> }</td></tr>
|
|
<tr><td class="num" id="LN192">192</td><td class="line"> <span class='keyword'>if</span>(iy==0x3ff00000) { <span class='comment'>/* y is +-1 */</span></td></tr>
|
|
<tr><td class="num" id="LN193">193</td><td class="line"> <span class='keyword'>if</span>(hy<0) <span class='keyword'>return</span> one/x;</td></tr>
|
|
<tr><td class="num" id="LN194">194</td><td class="line"> <span class='keyword'>else</span> <span class='keyword'>return</span> x;</td></tr>
|
|
<tr><td class="num" id="LN195">195</td><td class="line"> }</td></tr>
|
|
<tr><td class="num" id="LN196">196</td><td class="line"> <span class='keyword'>if</span>(hy==0x40000000) <span class='keyword'>return</span> x*x; <span class='comment'>/* y is 2 */</span></td></tr>
|
|
<tr><td class="num" id="LN197">197</td><td class="line"> <span class='keyword'>if</span>(hy==0x3fe00000) { <span class='comment'>/* y is 0.5 */</span></td></tr>
|
|
<tr><td class="num" id="LN198">198</td><td class="line"> <span class='keyword'>if</span>(hx>=0) <span class='comment'>/* x >= +0 */</span></td></tr>
|
|
<tr><td class="num" id="LN199">199</td><td class="line"> <span class='keyword'>return</span> sqrt(x);</td></tr>
|
|
<tr><td class="num" id="LN200">200</td><td class="line"> }</td></tr>
|
|
<tr><td class="num" id="LN201">201</td><td class="line"> }</td></tr>
|
|
<tr><td class="num" id="LN202">202</td><td class="line"> </td></tr>
|
|
<tr><td class="num" id="LN203">203</td><td class="line"> ax = fabs(x);</td></tr>
|
|
<tr><td class="num" id="LN204">204</td><td class="line"> <span class='comment'>/* special value of x */</span></td></tr>
|
|
<tr><td class="num" id="LN205">205</td><td class="line"> <span class='keyword'>if</span>(lx==0) {</td></tr>
|
|
<tr><td class="num" id="LN206">206</td><td class="line"> <span class='keyword'>if</span>(ix==0x7ff00000||ix==0||ix==0x3ff00000) {</td></tr>
|
|
<tr><td class="num" id="LN207">207</td><td class="line"> z = ax; <span class='comment'>/*x is +-0,+-inf,+-1*/</span></td></tr>
|
|
<tr><td class="num" id="LN208">208</td><td class="line"> <span class='keyword'>if</span>(hy<0) z = one/z; <span class='comment'>/* z = (1/|x|) */</span></td></tr>
|
|
<tr><td class="num" id="LN209">209</td><td class="line"> <span class='keyword'>if</span>(hx<0) {</td></tr>
|
|
<tr><td class="num" id="LN210">210</td><td class="line"> <span class='keyword'>if</span>(((ix-0x3ff00000)|yisint)==0) {</td></tr>
|
|
<tr><td class="num" id="LN211">211</td><td class="line"> z = (z-z)/(z-z); <span class='comment'>/* (-1)**non-int is NaN */</span></td></tr>
|
|
<tr><td class="num" id="LN212">212</td><td class="line"> } <span class='keyword'>else</span> <span class='keyword'>if</span>(yisint==1)</td></tr>
|
|
<tr><td class="num" id="LN213">213</td><td class="line"> z = -z; <span class='comment'>/* (x<0)**odd = -(|x|**odd) */</span></td></tr>
|
|
<tr><td class="num" id="LN214">214</td><td class="line"> }</td></tr>
|
|
<tr><td class="num" id="LN215">215</td><td class="line"> <span class='keyword'>return</span> z;</td></tr>
|
|
<tr><td class="num" id="LN216">216</td><td class="line"> }</td></tr>
|
|
<tr><td class="num" id="LN217">217</td><td class="line"> }</td></tr>
|
|
<tr><td class="num" id="LN218">218</td><td class="line"> </td></tr>
|
|
<tr><td class="num" id="LN219">219</td><td class="line"> n = (hx>>31)+1;</td></tr>
|
|
<tr><td class="num" id="LN220">220</td><td class="line"> </td></tr>
|
|
<tr><td class="num" id="LN221">221</td><td class="line"> <span class='comment'>/* (x<0)**(non-int) is NaN */</span></td></tr>
|
|
<tr><td class="num" id="LN222">222</td><td class="line"> <span class='keyword'>if</span>((n|yisint)==0) <span class='keyword'>return</span> (x-x)/(x-x);</td></tr>
|
|
<tr><td class="num" id="LN223">223</td><td class="line"> </td></tr>
|
|
<tr><td class="num" id="LN224">224</td><td class="line"> s = one; <span class='comment'>/* s (sign of result -ve**odd) = -1 else = 1 */</span></td></tr>
|
|
<tr><td class="num" id="LN225">225</td><td class="line"> <span class='keyword'>if</span>((n|(yisint-1))==0) s = -one;<span class='comment'>/* (-ve)**(odd int) */</span></td></tr>
|
|
<tr><td class="num" id="LN226">226</td><td class="line"> </td></tr>
|
|
<tr><td class="num" id="LN227">227</td><td class="line"> <span class='comment'>/* |y| is huge */</span></td></tr>
|
|
<tr><td class="num" id="LN228">228</td><td class="line"> <span class='keyword'>if</span>(iy>0x41e00000) { <span class='comment'>/* if |y| > 2**31 */</span></td></tr>
|
|
<tr><td class="num" id="LN229">229</td><td class="line"> <span class='keyword'>if</span>(iy>0x43f00000) { <span class='comment'>/* if |y| > 2**64, must o/uflow */</span></td></tr>
|
|
<tr><td class="num" id="LN230">230</td><td class="line"> <span class='keyword'>if</span>(ix<=0x3fefffff) <span class='keyword'>return</span> (hy<0)? huge*huge:tiny*tiny;</td></tr>
|
|
<tr><td class="num" id="LN231">231</td><td class="line"> <span class='keyword'>if</span>(ix>=0x3ff00000) <span class='keyword'>return</span> (hy>0)? huge*huge:tiny*tiny;</td></tr>
|
|
<tr><td class="num" id="LN232">232</td><td class="line"> }</td></tr>
|
|
<tr><td class="num" id="LN233">233</td><td class="line"> <span class='comment'>/* over/underflow if x is not close to one */</span></td></tr>
|
|
<tr><td class="num" id="LN234">234</td><td class="line"> <span class='keyword'>if</span>(ix<0x3fefffff) <span class='keyword'>return</span> (hy<0)? s*huge*huge:s*tiny*tiny;</td></tr>
|
|
<tr><td class="num" id="LN235">235</td><td class="line"> <span class='keyword'>if</span>(ix>0x3ff00000) <span class='keyword'>return</span> (hy>0)? s*huge*huge:s*tiny*tiny;</td></tr>
|
|
<tr><td class="num" id="LN236">236</td><td class="line"> <span class='comment'>/* now |1-x| is tiny <= 2**-20, suffice to compute</span></td></tr>
|
|
<tr><td class="num" id="LN237">237</td><td class="line"> <span class='comment'>log(x) by x-x^2/2+x^3/3-x^4/4 */</span></td></tr>
|
|
<tr><td class="num" id="LN238">238</td><td class="line"> t = ax-one; <span class='comment'>/* t has 20 trailing zeros */</span></td></tr>
|
|
<tr><td class="num" id="LN239">239</td><td class="line"> w = (t*t)*(0.5-t*(0.3333333333333333333333-t*0.25));</td></tr>
|
|
<tr><td class="num" id="LN240">240</td><td class="line"> u = ivln2_h*t; <span class='comment'>/* ivln2_h has 21 sig. bits */</span></td></tr>
|
|
<tr><td class="num" id="LN241">241</td><td class="line"> v = t*ivln2_l-w*ivln2;</td></tr>
|
|
<tr><td class="num" id="LN242">242</td><td class="line"> t1 = u+v;</td></tr>
|
|
<tr><td class="num" id="LN243">243</td><td class="line"> <span class='macro'>SET_LOW_WORD(t1,0)<span class='expansion'>do { ieee_double_shape_type sl_u; sl_u.value = (t1); sl_u.parts<br>.lsw = (0); (t1) = sl_u.value; } while (0)</span></span>;</td></tr>
|
|
<tr><td class="num" id="LN244">244</td><td class="line"> t2 = v-(t1-u);</td></tr>
|
|
<tr><td class="num" id="LN245">245</td><td class="line"> } <span class='keyword'>else</span> {</td></tr>
|
|
<tr><td class="num" id="LN246">246</td><td class="line"> <span class='keyword'>double</span> ss,s2,s_h,s_l,t_h,t_l;</td></tr>
|
|
<tr><td class="num" id="LN247">247</td><td class="line"> n = 0;</td></tr>
|
|
<tr><td class="num" id="LN248">248</td><td class="line"> <span class='comment'>/* take care subnormal number */</span></td></tr>
|
|
<tr><td class="num" id="LN249">249</td><td class="line"> <span class='keyword'>if</span>(ix<0x00100000)</td></tr>
|
|
<tr><td class="num" id="LN250">250</td><td class="line"> {</td></tr>
|
|
<tr><td class="num" id="LN251">251</td><td class="line"> ax *= two53;</td></tr>
|
|
<tr><td class="num" id="LN252">252</td><td class="line"> n -= 53;</td></tr>
|
|
<tr><td class="num" id="LN253">253</td><td class="line"> <span class='macro'>GET_HIGH_WORD(ix,ax)<span class='expansion'>do { ieee_double_shape_type gh_u; gh_u.value = (ax); (ix) = gh_u<br>.parts.msw; } while (0)</span></span>;</td></tr>
|
|
<tr><td class="num" id="LN254">254</td><td class="line"> }</td></tr>
|
|
<tr><td class="num" id="LN255">255</td><td class="line"> n += ((ix)>>20)-0x3ff;</td></tr>
|
|
<tr><td class="num" id="LN256">256</td><td class="line"> j = ix&0x000fffff;</td></tr>
|
|
<tr><td class="num" id="LN257">257</td><td class="line"> <span class='comment'>/* determine interval */</span></td></tr>
|
|
<tr><td class="num" id="LN258">258</td><td class="line"> ix = j|0x3ff00000; <span class='comment'>/* normalize ix */</span></td></tr>
|
|
<tr><td class="num" id="LN259">259</td><td class="line"> <span class='keyword'>if</span>(j<=0x3988E) k=0; <span class='comment'>/* |x|<sqrt(3/2) */</span></td></tr>
|
|
<tr><td class="num" id="LN260">260</td><td class="line"> <span class='keyword'>else</span> <span class='keyword'>if</span>(j<0xBB67A) k=1; <span class='comment'>/* |x|<sqrt(3) */</span></td></tr>
|
|
<tr><td class="num" id="LN261">261</td><td class="line"> <span class='keyword'>else</span> {</td></tr>
|
|
<tr><td class="num" id="LN262">262</td><td class="line"> k=0;</td></tr>
|
|
<tr><td class="num" id="LN263">263</td><td class="line"> n+=1;</td></tr>
|
|
<tr><td class="num" id="LN264">264</td><td class="line"> ix -= 0x00100000;</td></tr>
|
|
<tr><td class="num" id="LN265">265</td><td class="line"> }</td></tr>
|
|
<tr><td class="num" id="LN266">266</td><td class="line"> <span class='macro'>SET_HIGH_WORD(ax,ix)<span class='expansion'>do { ieee_double_shape_type sh_u; sh_u.value = (ax); sh_u.parts<br>.msw = (ix); (ax) = sh_u.value; } while (0)</span></span>;</td></tr>
|
|
<tr><td class="num" id="LN267">267</td><td class="line"> </td></tr>
|
|
<tr><td class="num" id="LN268">268</td><td class="line"> <span class='comment'>/* compute ss = s_h+s_l = (x-1)/(x+1) or (x-1.5)/(x+1.5) */</span></td></tr>
|
|
<tr><td class="num" id="LN269">269</td><td class="line"> u = ax-bp[k]; <span class='comment'>/* bp[0]=1.0, bp[1]=1.5 */</span></td></tr>
|
|
<tr><td class="num" id="LN270">270</td><td class="line"> v = one/(ax+bp[k]);</td></tr>
|
|
<tr><td class="num" id="LN271">271</td><td class="line"> ss = u*v;</td></tr>
|
|
<tr><td class="num" id="LN272">272</td><td class="line"> s_h = ss;</td></tr>
|
|
<tr><td class="num" id="LN273">273</td><td class="line"> <span class='macro'>SET_LOW_WORD(s_h,0)<span class='expansion'>do { ieee_double_shape_type sl_u; sl_u.value = (s_h); sl_u.parts<br>.lsw = (0); (s_h) = sl_u.value; } while (0)</span></span>;</td></tr>
|
|
<tr><td class="num" id="LN274">274</td><td class="line"> <span class='comment'>/* t_h=ax+bp[k] High */</span></td></tr>
|
|
<tr><td class="num" id="LN275">275</td><td class="line"> t_h = zero;</td></tr>
|
|
<tr><td class="num" id="LN276">276</td><td class="line"> <span class='macro'>SET_HIGH_WORD(t_h,((ix>>1)|0x20000000)+0x00080000+(k<<18))<span class='expansion'>do { ieee_double_shape_type sh_u; sh_u.value = (t_h); sh_u.parts<br>.msw = (((ix>>1)|0x20000000)+0x00080000+(k<<18));<br> (t_h) = sh_u.value; } while (0)</span></span>;</td></tr>
|
|
<tr><td class="num" id="LN277">277</td><td class="line"> t_l = ax - (t_h-bp[k]);</td></tr>
|
|
<tr><td class="num" id="LN278">278</td><td class="line"> s_l = v*((u-s_h*t_h)-s_h*t_l);</td></tr>
|
|
<tr><td class="num" id="LN279">279</td><td class="line"> <span class='comment'>/* compute log(ax) */</span></td></tr>
|
|
<tr><td class="num" id="LN280">280</td><td class="line"> s2 = ss*ss;</td></tr>
|
|
<tr><td class="num" id="LN281">281</td><td class="line"> r = s2*s2*(L1+s2*(L2+s2*(L3+s2*(L4+s2*(L5+s2*L6)))));</td></tr>
|
|
<tr><td class="num" id="LN282">282</td><td class="line"> r += s_l*(s_h+ss);</td></tr>
|
|
<tr><td class="num" id="LN283">283</td><td class="line"> s2 = s_h*s_h;</td></tr>
|
|
<tr><td class="num" id="LN284">284</td><td class="line"> t_h = 3.0+s2+r;</td></tr>
|
|
<tr><td class="num" id="LN285">285</td><td class="line"> <span class='macro'>SET_LOW_WORD(t_h,0)<span class='expansion'>do { ieee_double_shape_type sl_u; sl_u.value = (t_h); sl_u.parts<br>.lsw = (0); (t_h) = sl_u.value; } while (0)</span></span>;</td></tr>
|
|
<tr><td class="num" id="LN286">286</td><td class="line"> t_l = r-((t_h-3.0)-s2);</td></tr>
|
|
<tr><td class="num" id="LN287">287</td><td class="line"> <span class='comment'>/* u+v = ss*(1+...) */</span></td></tr>
|
|
<tr><td class="num" id="LN288">288</td><td class="line"> u = s_h*t_h;</td></tr>
|
|
<tr><td class="num" id="LN289">289</td><td class="line"> v = s_l*t_h+t_l*ss;</td></tr>
|
|
<tr><td class="num" id="LN290">290</td><td class="line"> <span class='comment'>/* 2/(3log2)*(ss+...) */</span></td></tr>
|
|
<tr><td class="num" id="LN291">291</td><td class="line"> p_h = u+v;</td></tr>
|
|
<tr><td class="num" id="LN292">292</td><td class="line"> <span class='macro'>SET_LOW_WORD(p_h,0)<span class='expansion'>do { ieee_double_shape_type sl_u; sl_u.value = (p_h); sl_u.parts<br>.lsw = (0); (p_h) = sl_u.value; } while (0)</span></span>;</td></tr>
|
|
<tr><td class="num" id="LN293">293</td><td class="line"> p_l = v-(p_h-u);</td></tr>
|
|
<tr><td class="num" id="LN294">294</td><td class="line"> z_h = cp_h*p_h; <span class='comment'>/* cp_h+cp_l = 2/(3*log2) */</span></td></tr>
|
|
<tr><td class="num" id="LN295">295</td><td class="line"> z_l = cp_l*p_h+p_l*cp+dp_l[k];</td></tr>
|
|
<tr><td class="num" id="LN296">296</td><td class="line"> <span class='comment'>/* log2(ax) = (ss+..)*2/(3*log2) = n + dp_h + z_h + z_l */</span></td></tr>
|
|
<tr><td class="num" id="LN297">297</td><td class="line"> t = (<span class='keyword'>double</span>)n;</td></tr>
|
|
<tr><td class="num" id="LN298">298</td><td class="line"> t1 = (((z_h+z_l)+dp_h[k])+t);</td></tr>
|
|
<tr><td class="num" id="LN299">299</td><td class="line"> <span class='macro'>SET_LOW_WORD(t1,0)<span class='expansion'>do { ieee_double_shape_type sl_u; sl_u.value = (t1); sl_u.parts<br>.lsw = (0); (t1) = sl_u.value; } while (0)</span></span>;</td></tr>
|
|
<tr><td class="num" id="LN300">300</td><td class="line"> t2 = z_l-(((t1-t)-dp_h[k])-z_h);</td></tr>
|
|
<tr><td class="num" id="LN301">301</td><td class="line"> }</td></tr>
|
|
<tr><td class="num" id="LN302">302</td><td class="line"> </td></tr>
|
|
<tr><td class="num" id="LN303">303</td><td class="line"> <span class='comment'>/* split up y into y1+y2 and compute (y1+y2)*(t1+t2) */</span></td></tr>
|
|
<tr><td class="num" id="LN304">304</td><td class="line"> y1 = y;</td></tr>
|
|
<tr><td class="num" id="LN305">305</td><td class="line"> <span class='macro'>SET_LOW_WORD(y1,0)<span class='expansion'>do { ieee_double_shape_type sl_u; sl_u.value = (y1); sl_u.parts<br>.lsw = (0); (y1) = sl_u.value; } while (0)</span></span>;</td></tr>
|
|
<tr><td class="num" id="LN306">306</td><td class="line"> p_l = (y-y1)*t1+y*t2;</td></tr>
|
|
<tr><td class="num" id="LN307">307</td><td class="line"> p_h = y1*t1;</td></tr>
|
|
<tr><td class="num" id="LN308">308</td><td class="line"> z = p_l+p_h;</td></tr>
|
|
<tr><td class="num" id="LN309">309</td><td class="line"> <span class='macro'>EXTRACT_WORDS(j,i,z)<span class='expansion'>do { ieee_double_shape_type ew_u; ew_u.value = (z); (j) = ew_u<br>.parts.msw; (i) = ew_u.parts.lsw; } while (0)</span></span>;</td></tr>
|
|
<tr><td class="num" id="LN310">310</td><td class="line"> <span class='keyword'>if</span> (j>=0x40900000) { <span class='comment'>/* z >= 1024 */</span></td></tr>
|
|
<tr><td class="num" id="LN311">311</td><td class="line"> <span class='keyword'>if</span>(((j-0x40900000)|i)!=0) <span class='comment'>/* if z > 1024 */</span></td></tr>
|
|
<tr><td class="num" id="LN312">312</td><td class="line"> <span class='keyword'>return</span> s*huge*huge; <span class='comment'>/* overflow */</span></td></tr>
|
|
<tr><td class="num" id="LN313">313</td><td class="line"> <span class='keyword'>else</span> {</td></tr>
|
|
<tr><td class="num" id="LN314">314</td><td class="line"> <span class='keyword'>if</span>(p_l+ovt>z-p_h) <span class='keyword'>return</span> s*huge*huge; <span class='comment'>/* overflow */</span></td></tr>
|
|
<tr><td class="num" id="LN315">315</td><td class="line"> }</td></tr>
|
|
<tr><td class="num" id="LN316">316</td><td class="line"> } <span class='keyword'>else</span> <span class='keyword'>if</span>((j&0x7fffffff)>=0x4090cc00 ) { <span class='comment'>/* z <= -1075 */</span></td></tr>
|
|
<tr><td class="num" id="LN317">317</td><td class="line"> <span class='keyword'>if</span>(((j-0xc090cc00)|i)!=0) <span class='comment'>/* z < -1075 */</span></td></tr>
|
|
<tr><td class="num" id="LN318">318</td><td class="line"> <span class='keyword'>return</span> s*tiny*tiny; <span class='comment'>/* underflow */</span></td></tr>
|
|
<tr><td class="num" id="LN319">319</td><td class="line"> <span class='keyword'>else</span> {</td></tr>
|
|
<tr><td class="num" id="LN320">320</td><td class="line"> <span class='keyword'>if</span>(p_l<=z-p_h) <span class='keyword'>return</span> s*tiny*tiny; <span class='comment'>/* underflow */</span></td></tr>
|
|
<tr><td class="num" id="LN321">321</td><td class="line"> }</td></tr>
|
|
<tr><td class="num" id="LN322">322</td><td class="line"> }</td></tr>
|
|
<tr><td class="num" id="LN323">323</td><td class="line"> <span class='comment'>/*</span></td></tr>
|
|
<tr><td class="num" id="LN324">324</td><td class="line"> <span class='comment'>* compute 2**(p_h+p_l)</span></td></tr>
|
|
<tr><td class="num" id="LN325">325</td><td class="line"> <span class='comment'>*/</span></td></tr>
|
|
<tr><td class="num" id="LN326">326</td><td class="line"> i = j&0x7fffffff;</td></tr>
|
|
<tr><td class="num" id="LN327">327</td><td class="line"> k = (i>>20)-0x3ff;</td></tr>
|
|
<tr><td class="num" id="LN328">328</td><td class="line"> n = 0;</td></tr>
|
|
<tr><td class="num" id="LN329">329</td><td class="line"> <span class='keyword'>if</span>(i>0x3fe00000) { <span class='comment'>/* if |z| > 0.5, set n = [z+0.5] */</span></td></tr>
|
|
<tr><td class="num" id="LN330">330</td><td class="line"> n = j+(0x00100000>>(k+1));</td></tr>
|
|
<tr><td class="num" id="LN331">331</td><td class="line"> k = ((n&0x7fffffff)>>20)-0x3ff; <span class='comment'>/* new k for n */</span></td></tr>
|
|
<tr><td class="num" id="LN332">332</td><td class="line"> t = zero;</td></tr>
|
|
<tr><td class="num" id="LN333">333</td><td class="line"> <span class='macro'>SET_HIGH_WORD(t,(n&~(0x000fffff>>k)))<span class='expansion'>do { ieee_double_shape_type sh_u; sh_u.value = (t); sh_u.parts<br>.msw = ((n&~(0x000fffff>>k))); (t) = sh_u.value; } while<br> (0)</span></span>;</td></tr>
|
|
<tr><td class="num" id="LN334">334</td><td class="line"> n = ((n&0x000fffff)|0x00100000)>>(20-k);</td></tr>
|
|
<tr><td class="num" id="LN335">335</td><td class="line"> <span class='keyword'>if</span>(j<0) n = -n;</td></tr>
|
|
<tr><td class="num" id="LN336">336</td><td class="line"> p_h -= t;</td></tr>
|
|
<tr><td class="num" id="LN337">337</td><td class="line"> }</td></tr>
|
|
<tr><td class="num" id="LN338">338</td><td class="line"> t = p_l+p_h;</td></tr>
|
|
<tr><td class="num" id="LN339">339</td><td class="line"> <span class='macro'>SET_LOW_WORD(t,0)<span class='expansion'>do { ieee_double_shape_type sl_u; sl_u.value = (t); sl_u.parts<br>.lsw = (0); (t) = sl_u.value; } while (0)</span></span>;</td></tr>
|
|
<tr><td class="num" id="LN340">340</td><td class="line"> u = t*lg2_h;</td></tr>
|
|
<tr><td class="num" id="LN341">341</td><td class="line"> v = (p_l-(t-p_h))*lg2+t*lg2_l;</td></tr>
|
|
<tr><td class="num" id="LN342">342</td><td class="line"> z = u+v;</td></tr>
|
|
<tr><td class="num" id="LN343">343</td><td class="line"> w = v-(z-u);</td></tr>
|
|
<tr><td class="num" id="LN344">344</td><td class="line"> t = z*z;</td></tr>
|
|
<tr><td class="num" id="LN345">345</td><td class="line"> t1 = z - t*(P1+t*(P2+t*(P3+t*(P4+t*P5))));</td></tr>
|
|
<tr><td class="num" id="LN346">346</td><td class="line"> r = (z*t1)/(t1-two)-(w+z*w);</td></tr>
|
|
<tr><td class="num" id="LN347">347</td><td class="line"> z = one-(r-z);</td></tr>
|
|
<tr><td class="num" id="LN348">348</td><td class="line"> <span class='macro'>GET_HIGH_WORD(j,z)<span class='expansion'>do { ieee_double_shape_type gh_u; gh_u.value = (z); (j) = gh_u<br>.parts.msw; } while (0)</span></span>;</td></tr>
|
|
<tr><td class="num" id="LN349">349</td><td class="line"> j += (n<<20);</td></tr>
|
|
<tr><td class="num" id="LN350">350</td><td class="line"> <span class='keyword'>if</span>((j>>20)<=0) z = scalbn(z,n); <span class='comment'>/* subnormal output */</span></td></tr>
|
|
<tr><td class="num" id="LN351">351</td><td class="line"> <span class='keyword'>else</span></td></tr>
|
|
<tr><td class="num" id="LN352">352</td><td class="line"> {</td></tr>
|
|
<tr><td class="num" id="LN353">353</td><td class="line"> uword hz;</td></tr>
|
|
<tr><td class="num" id="LN354">354</td><td class="line"> <span class='macro'>GET_HIGH_WORD(hz,z)<span class='expansion'>do { ieee_double_shape_type gh_u; gh_u.value = (z); (hz) = gh_u<br>.parts.msw; } while (0)</span></span>;</td></tr>
|
|
<tr><td class="num" id="LN355">355</td><td class="line"> <span class='macro'>SET_HIGH_WORD(z,hz + (n<<20))<span class='expansion'>do { ieee_double_shape_type sh_u; sh_u.value = (z); sh_u.parts<br>.msw = (hz + (n<<20)); (z) = sh_u.value; } while (0)</span></span>;</td></tr>
|
|
<tr><td class="num" id="LN356">356</td><td class="line"> }</td></tr>
|
|
<tr><td class="num" id="LN357">357</td><td class="line"> <span class='keyword'>return</span> s*z;</td></tr>
|
|
<tr><td class="num" id="LN358">358</td><td class="line">}</td></tr>
|
|
</table></body></html>
|