Sophie

Sophie

distrib > Fedora > 15 > i386 > by-pkgid > ff8e7344076b5fbaa54d805766a057bd > files > 26

libscs-devel-1.4.1-4.fc15.i686.rpm

<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.0 Transitional//EN">
<html><head><meta http-equiv="Content-Type" content="text/html;charset=iso-8859-1">
<title>scs2doubleCLEAN_BUT_SLOW.c Source File</title>
<link href="doxygen.css" rel="stylesheet" type="text/css">
</head><body>
<!-- Generated by Doxygen 1.2.15 -->
<center>
<a class="qindex" href="index.html">Main Page</a> &nbsp; <a class="qindex" href="annotated.html">Data Structures</a> &nbsp; <a class="qindex" href="files.html">File List</a> &nbsp; <a class="qindex" href="functions.html">Data Fields</a> &nbsp; <a class="qindex" href="globals.html">Globals</a> &nbsp; </center>
<hr><h1>scs2doubleCLEAN_BUT_SLOW.c</h1><div class="fragment"><pre>00001 <font class="comment">/** Conversion of SCS to floating-point double </font>
00002 <font class="comment">@file scs2double.c</font>
00003 <font class="comment"></font>
00004 <font class="comment">@author Defour David David.Defour@ens-lyon.fr</font>
00005 <font class="comment">@author Florent de Dinechin Florent.de.Dinechin@ens-lyon.fr </font>
00006 <font class="comment"></font>
00007 <font class="comment">This file is part of the SCS library.</font>
00008 <font class="comment">*/</font>
00009 
00010 <font class="comment">/*</font>
00011 <font class="comment">Copyright (C) 2002  David Defour and Florent de Dinechin</font>
00012 <font class="comment"></font>
00013 <font class="comment">    This library is free software; you can redistribute it and/or</font>
00014 <font class="comment">    modify it under the terms of the GNU Lesser General Public</font>
00015 <font class="comment">    License as published by the Free Software Foundation; either</font>
00016 <font class="comment">    version 2.1 of the License, or (at your option) any later version.</font>
00017 <font class="comment"></font>
00018 <font class="comment">    This library is distributed in the hope that it will be useful,</font>
00019 <font class="comment">    but WITHOUT ANY WARRANTY; without even the implied warranty of</font>
00020 <font class="comment">    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU</font>
00021 <font class="comment">    Lesser General Public License for more details.</font>
00022 <font class="comment"></font>
00023 <font class="comment">    You should have received a copy of the GNU Lesser General Public</font>
00024 <font class="comment">    License along with this library; if not, write to the Free Software</font>
00025 <font class="comment">    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA</font>
00026 <font class="comment"></font>
00027 <font class="comment"> */</font>
00028 <font class="preprocessor">#include "<a class="code" href="scs_8h.html">scs.h</a>"</font>
00029 <font class="preprocessor">#include "<a class="code" href="scs__private_8h.html">scs_private.h</a>"</font>
00030 
00031 <font class="comment"></font>
00032 <font class="comment">/** Convert a multiple precision number in scs format into a double</font>
00033 <font class="comment"> precision number.</font>
00034 <font class="comment"></font>
00035 <font class="comment">@warning  "x" need to be normalized</font>
00036 <font class="comment">  */</font> 
00037 
00038 
00039 
00040 
00041 
00042 
00043 
00044 
00045 
00046 
00047 
00048 
00049 
00050 <font class="comment">/*  computes the exponent from the index */</font>
00051 <font class="keyword">inline</font>  <font class="keywordtype">void</font> computeExponent(<a class="code" href="structscs.html">scs_ptr</a> x, <font class="keywordtype">double</font> *res) {
00052   scs_db_number nb;
00053   <font class="keywordtype">int</font> i;
00054   <font class="keywordflow">if</font> ((X_IND &lt; SCS_MAX_RANGE) &amp;&amp; (X_IND &gt; -SCS_MAX_RANGE)){
00055     <font class="comment">/* x is comfortably in the double-precision range */</font>
00056     <font class="comment">/* build the double 2^(X_IND*SCS_NB_BITS)   */</font>
00057     nb.i[HI_ENDIAN] = ((X_IND)*SCS_NB_BITS +1023)  &lt;&lt; 20;
00058     nb.i[LO_ENDIAN] = 0;
00059     *res *= nb.d;
00060   }<font class="keywordflow">else</font> {
00061     <font class="comment">/* x may end up being a denormal or overflow */</font>
00062     i    = X_IND;
00063     nb.d = 0;
00064     <font class="keywordflow">if</font> (X_IND &gt; 0){
00065       <font class="comment">/* one of the following computations may lead to an overflow */</font>
00066       *res *=SCS_RADIX_RNG_DOUBLE; <font class="comment">/* 2^(SCS_NB_BITS.SCS_MAX_RANGE) */</font>
00067       i   -= SCS_MAX_RANGE;
00068       <font class="keywordflow">while</font>((i--&gt;0)&amp;&amp;(*res &lt;= MAX_DOUBLE)) {
00069         <font class="comment">/* Thanks to the second test, this loop stops when an overflow is encountered */</font>
00070         *res *= SCS_RADIX_ONE_DOUBLE;
00071       }
00072     }<font class="keywordflow">else</font> {
00073       <font class="comment">/* One of the following computations may need to a denormal/underflow */</font>
00074       *res *=SCS_RADIX_MRNG_DOUBLE; <font class="comment">/* 2^-(SCS_NB_BITS.SCS_MAX_RANGE)*/</font>
00075       i   += SCS_MAX_RANGE;
00076       <font class="keywordflow">while</font>((i++&lt;0)&amp;&amp;(*res != 0)) {
00077         *res *=SCS_RADIX_MONE_DOUBLE;
00078       }
00079     }                    
00080   }
00081   <font class="comment">/* sign management */</font>
00082   <font class="keywordflow">if</font> (X_SGN &lt; 0) 
00083     *res = - *res;
00084 }
00085 
00086 
00087 
00088 <font class="keywordtype">void</font> <a class="code" href="scs2double_8c.html#a0">scs_get_d</a>(<font class="keywordtype">double</font> *result, <a class="code" href="structscs.html">scs_ptr</a> x){ 
00089   scs_db_number nb, rndcorr;
00090   <font class="keywordtype">unsigned</font> <font class="keywordtype">long</font> <font class="keywordtype">long</font> <font class="keywordtype">int</font> lowpart, t1;
00091   <font class="keywordtype">int</font> exp;
00092   
00093   <font class="comment">/* convert the MSB digit into a double, and store it in nb.d */</font>
00094   nb.d = (double)X_HW[0]; 
00095 
00096   <font class="comment">/* place the two next digits in lowpart */</font>
00097   t1   = X_HW[1];
00098   lowpart  = (t1 &lt;&lt; SCS_NB_BITS) + X_HW[2];
00099 
00100   <font class="comment">/* test for  s/qNan, +/- Inf, +/- 0, placed here for obscure performance reasons */</font>
00101   <font class="keywordflow">if</font> (X_EXP != 1){*result = X_EXP; <font class="keywordflow">return</font>;}
00102   
00103   <font class="comment">/* take the exponent of nb.d (will be in [0:SCS_NB_BITS])*/</font>
00104   exp = ((nb.i[HI_ENDIAN] &amp; 0x7ff00000)&gt;&gt;20) - 1023; 
00105   <font class="comment">/* align the rest of the mantissa to nb : shift by (2*SCS_NB_BITS)-53-exp */</font>
00106   lowpart = lowpart &gt;&gt; (exp+(2*SCS_NB_BITS)-53);     
00107   
00108   <font class="comment">/* Now look at the last bit to decide rounding */</font>
00109   <font class="keywordflow">if</font> (lowpart &amp; 0x0000000000000001){
00110    <font class="comment">/* need to add an half-ulp */</font>
00111     rndcorr.i[LO_ENDIAN] = 0; 
00112     rndcorr.i[HI_ENDIAN] = (exp+971) &lt;&lt; 20;    <font class="comment">/*((exp-52)+1023) &lt;&lt; 20;*/</font>
00113   }<font class="keywordflow">else</font>{
00114     <font class="comment">/* need to add nothing*/</font>
00115     rndcorr.l = 0;
00116   }
00117   
00118   lowpart = lowpart &gt;&gt; 1;
00119   nb.l = nb.l | lowpart;    <font class="comment">/* Finish to fill the mantissa */</font>
00120   *result  = nb.d + rndcorr.d;  <font class="comment">/* rounded to nearest   */</font>
00121  
00122 
00123   <font class="comment">/* now compute the exponent from the index */</font>
00124   computeExponent(x, result);
00125 
00126   <font class="keywordflow">return</font>;
00127 }
00128 
00129 
00130 
00131 <font class="keywordtype">void</font> scs_get_d_directed(<a class="code" href="structscs.html">scs_ptr</a> x, <font class="keywordtype">int</font> rndUp, <font class="keywordtype">double</font>* res)
00132  {
00133   <font class="keywordtype">unsigned</font> <font class="keywordtype">long</font> <font class="keywordtype">long</font> <font class="keywordtype">int</font> lowpart, t1;
00134   scs_db_number nb, rndcorr;
00135   <font class="keywordtype">int</font> i, exp, not_null;
00136 
00137   nb.d = (double)X_HW[0];
00138 
00139   t1   = X_HW[1];
00140   lowpart  = (t1 &lt;&lt; SCS_NB_BITS) + X_HW[2];
00141 
00142   <font class="comment">/* s/qNan, +/- Inf, +/- 0 */</font>
00143   <font class="keywordflow">if</font> (X_EXP != 1){
00144     *res = X_EXP; 
00145     <font class="keywordflow">return</font>;
00146   }
00147   
00148   <font class="comment">/* take the exponent */</font>
00149   exp = ((nb.i[HI_ENDIAN] &amp; 0x7ff00000)&gt;&gt;20) - 1023; 
00150   
00151   not_null = ((lowpart &lt;&lt; (64+52 - 2*SCS_NB_BITS - exp))==0) ? 0 : 1;
00152 
00153   <font class="comment">/* align the rest of the mantissa  */</font>
00154   lowpart = lowpart &gt;&gt; (exp + 2*SCS_NB_BITS - 52);     
00155 
00156  <font class="comment">/* Finish to fill the mantissa */</font>
00157   nb.l = nb.l | lowpart;   
00158 
00159   <font class="comment">/* Test if we are not on an exact double precision number */</font>
00160   <font class="keywordflow">for</font> (i=3; i&lt;SCS_NB_WORDS; i++)
00161       <font class="keywordflow">if</font> (X_HW[i]!=0)  not_null = 1;
00162 
00163   <font class="keywordflow">if</font> (rndUp &amp;&amp; (not_null)){
00164    <font class="comment">/* we have to perform a rounded toward +inf (ulp &gt;= 0.5) */</font>
00165     rndcorr.i[LO_ENDIAN] = 0; 
00166     rndcorr.i[HI_ENDIAN] = (exp+971)&lt;&lt;20;    <font class="comment">/*((exp-52)+1023) &lt;&lt; 20;*/</font>
00167   } <font class="keywordflow">else</font> {
00168       rndcorr.d = 0.0;
00169   }
00170 
00171   *res  = nb.d + rndcorr.d;  <font class="comment">/* make a rounded to nearest   */</font>
00172   computeExponent(x, res);
00173 }
00174 
00175 
00176 <font class="comment">/*</font>
00177 <font class="comment"> * Rounded toward -Inf</font>
00178 <font class="comment"> */</font>
00179 <font class="keywordtype">void</font> scs_get_d_minf(<font class="keywordtype">double</font> *result, <a class="code" href="structscs.html">scs_ptr</a> x){ 
00180 
00181   <font class="comment">/* round up the mantissa if negative  */</font>
00182   scs_get_d_directed(x, (X_SGN&lt;0), result);
00183  
00184   <font class="keywordflow">return</font>;
00185 }
00186 
00187 
00188 
00189 <font class="comment">/*</font>
00190 <font class="comment"> * Rounded toward +Inf</font>
00191 <font class="comment"> */</font>
00192 <font class="keywordtype">void</font> scs_get_d_pinf(<font class="keywordtype">double</font> *result, <a class="code" href="structscs.html">scs_ptr</a> x){ 
00193 
00194   <font class="comment">/* round up the mantissa if positive  */</font>
00195   scs_get_d_directed(x, (X_SGN&gt;=0), result);
00196  
00197   <font class="keywordflow">return</font>;
00198 }
00199 
00200 
00201 
00202 <font class="comment">/*</font>
00203 <font class="comment"> * Rounded toward zero</font>
00204 <font class="comment"> */</font>
00205 <font class="keywordtype">void</font> scs_get_d_zero(<font class="keywordtype">double</font> *result, <a class="code" href="structscs.html">scs_ptr</a> x){ 
00206   <font class="comment">/* never round up the mantissa  */</font>
00207   scs_get_d_directed(x, 0, result);
00208 }
</pre></div><hr><address align="right"><small>Generated on Tue Jun 17 10:15:51 2003 for SCSLib by
<a href="http://www.doxygen.org/index.html">
<img src="doxygen.png" alt="doxygen" align="middle" border=0 
width=110 height=53></a>1.2.15 </small></address>
</body>
</html>