<!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>scs2double.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> <a class="qindex" href="annotated.html">Data Structures</a> <a class="qindex" href="files.html">File List</a> <a class="qindex" href="functions.html">Data Fields</a> <a class="qindex" href="globals.html">Globals</a> </center>
<hr><h1>scs2double.c</h1><a href="scs2double_8c.html">Go to the documentation of this file.</a><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 <font class="comment">/* computes the exponent from the index */</font>
00042 <font class="comment">/* in principle an inline function would be cleaner, but </font>
00043 <font class="comment"> this leads to faster and smaller code </font>
00044 <font class="comment">*/</font>
00045
00046
00047
00048
00049
<a name="l00050"></a><a class="code" href="scs2double_8c.html#a0">00050</a> <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){
00051 scs_db_number nb, rndcorr;
00052 <font class="keywordtype">unsigned</font> <font class="keywordtype">long</font> <font class="keywordtype">long</font> <font class="keywordtype">int</font> lowpart, t1;
00053 <font class="keywordtype">int</font> exp,expfinal;
00054 <font class="keywordtype">double</font> res;
00055
00056 <font class="comment">/* convert the MSB digit into a double, and store it in nb.d */</font>
00057 nb.d = (double)X_HW[0];
00058
00059 <font class="comment">/* place the two next digits in lowpart */</font>
00060 t1 = X_HW[1];
00061 lowpart = (t1 << SCS_NB_BITS) + X_HW[2];
00062 <font class="comment">/* there is at least one significant bit in nb, </font>
00063 <font class="comment"> and at least 2*SCS_NB_BITS in lowpart, </font>
00064 <font class="comment"> so provided SCS_NB_BITS >= 27</font>
00065 <font class="comment"> together they include the 53+ guard bits to decide rounding </font>
00066 <font class="comment"> */</font>
00067
00068 <font class="comment">/* test for s/qNan, +/- Inf, +/- 0, placed here for obscure performance reasons */</font>
00069 <font class="keywordflow">if</font> (X_EXP != 1){
00070 *result = X_EXP;
00071 <font class="keywordflow">return</font>;
00072 }
00073
00074 <font class="comment">/* take the exponent of nb.d (will be in [0:SCS_NB_BITS])*/</font>
00075 exp = ((nb.i[HI_ENDIAN] & 0x7ff00000)>>20) - 1023;
00076
00077 <font class="comment">/* compute the exponent of the result */</font>
00078 expfinal = exp + SCS_NB_BITS*X_IND;
00079
00080 <font class="comment">/* Is the SCS number not too large for the IEEE exponent range ? */</font>
00081 <font class="keywordflow">if</font> (expfinal > 1023) {
00082 <font class="comment">/* return an infinity */</font>
00083 res = SCS_RADIX_RNG_DOUBLE*SCS_RADIX_RNG_DOUBLE;
00084 }
00085
00086 <font class="comment">/* Is our SCS number a denormal ? */</font>
00087 <font class="keywordflow">else</font> <font class="keywordflow">if</font> (expfinal >= -1022){
00088 <font class="comment">/* x is in the normal range */</font>
00089
00090 <font class="comment">/* align the rest of the mantissa to nb : shift by (2*SCS_NB_BITS)-53-exp */</font>
00091 lowpart = lowpart >> (exp+(2*SCS_NB_BITS)-53);
00092 <font class="comment">/* Look at the last bit to decide rounding */</font>
00093 <font class="keywordflow">if</font> (lowpart & 0x0000000000000001){
00094 <font class="comment">/* need to add an half-ulp */</font>
00095 rndcorr.i[LO_ENDIAN] = 0;
00096 rndcorr.i[HI_ENDIAN] = (exp-52+1023)<<20; <font class="comment">/* 2^(exp-52) */</font>
00097 }<font class="keywordflow">else</font>{
00098 <font class="comment">/* need to add nothing*/</font>
00099 rndcorr.d = 0.0;
00100 }
00101
00102 lowpart = lowpart >> 1;
00103 nb.l = nb.l | lowpart; <font class="comment">/* Finish to fill the mantissa */</font>
00104 res = nb.d + rndcorr.d; <font class="comment">/* rounded to nearest */</font>
00105
00106 <font class="comment">/* now compute the exponent from the index :</font>
00107 <font class="comment"> we need to multiply res by 2^(X_IND*SCS_NB_BITS)</font>
00108 <font class="comment"> First check this number won't be a denormal itself */</font>
00109 <font class="keywordflow">if</font>((X_IND)*SCS_NB_BITS +1023 > 0) {
00110 <font class="comment">/* build the double 2^(X_IND*SCS_NB_BITS) */</font>
00111 nb.i[HI_ENDIAN] = ((X_IND)*SCS_NB_BITS +1023) << 20;
00112 nb.i[LO_ENDIAN] = 0;
00113 res *= nb.d; <font class="comment">/* exact multiplication */</font>
00114 }
00115 <font class="keywordflow">else</font> { <font class="comment">/*offset the previous computation by 2^(2*SCS_NB_BITS) */</font>
00116 <font class="comment">/* build the double 2^(X_IND*SCS_NB_BITS) */</font>
00117 nb.i[HI_ENDIAN] = ((X_IND)*SCS_NB_BITS +1023 + 2*SCS_NB_BITS) << 20;
00118 nb.i[LO_ENDIAN] = 0;
00119 res *= SCS_RADIX_MTWO_DOUBLE; <font class="comment">/* exact multiplication */</font>
00120 res *= nb.d; <font class="comment">/* exact multiplication */</font>
00121 }
00122 }
00123
00124
00125 <font class="keywordflow">else</font> {
00126 <font class="comment">/* the final number is a denormal with 52-(expfinal+1022)</font>
00127 <font class="comment"> significant bits. */</font>
00128
00129 <font class="keywordflow">if</font> (expfinal < -1022 - 53 ) {
00130 res = 0.0;
00131 }
00132 <font class="keywordflow">else</font> {
00133
00134 <font class="comment">/* align the rest of the mantissa to nb */</font>
00135 lowpart = lowpart >> (exp+(2*SCS_NB_BITS)-52);
00136 <font class="comment">/* Finish to fill the mantissa */</font>
00137 nb.l = nb.l | lowpart;
00138
00139 <font class="comment">/* this is still a normal number. </font>
00140 <font class="comment"> Now remove its exponent and add back the implicit one */</font>
00141 nb.l = (nb.l & 0x000FFFFFFFFFFFFF) | 0x0010000000000000;
00142
00143 <font class="comment">/* keep only the significant bits */</font>
00144 nb.l = nb.l >> (-1023 - expfinal);
00145 <font class="comment">/* Look at the last bit to decide rounding */</font>
00146 <font class="keywordflow">if</font> (nb.i[LO_ENDIAN] & 0x00000001){
00147 <font class="comment">/* need to add an half-ulp */</font>
00148 rndcorr.l = 1; <font class="comment">/* this is a full ulp but we multiply by 0.5 in the end */</font>
00149 }<font class="keywordflow">else</font>{
00150 <font class="comment">/* need to add nothing*/</font>
00151 rndcorr.d = 0.0;
00152
00153 }
00154 res = 0.5*(nb.d + rndcorr.d); <font class="comment">/* rounded to nearest */</font>
00155
00156 <font class="comment">/* the exponent field is already set to zero so that's all */</font>
00157 }
00158 }
00159
00160 <font class="comment">/* sign management */</font>
00161 <font class="keywordflow">if</font> (X_SGN < 0)
00162 *result = - res;
00163 <font class="keywordflow">else</font>
00164 *result = res;
00165 }
00166
00167
00168
00169
00170
00171 <font class="comment">/* All the directed roundings boil down to the same computation, which</font>
00172 <font class="comment">is: first build the truncated mantissa. if the SCS number is exactly</font>
00173 <font class="comment">a double precision number, return that. Otherwise, either return the</font>
00174 <font class="comment">truncated mantissa, or return this mantissa plus an ulp, rounded to</font>
00175 <font class="comment">the nearest. Plus handle the infinities and denormals.</font>
00176 <font class="comment">*/</font>
00177
00178 <font class="keywordtype">void</font> get_d_directed(<font class="keywordtype">double</font> *result, <a class="code" href="structscs.html">scs_ptr</a> x, <font class="keywordtype">int</font> rndMantissaUp){
00179 scs_db_number nb, rndcorr;
00180 <font class="keywordtype">unsigned</font> <font class="keywordtype">long</font> <font class="keywordtype">long</font> <font class="keywordtype">int</font> lowpart, t1;
00181 <font class="keywordtype">int</font> exp,expfinal,i, not_null;
00182 <font class="keywordtype">double</font> res;
00183
00184 <font class="comment">/* convert the MSB digit into a double, and store it in nb.d */</font>
00185 nb.d = (double)X_HW[0];
00186
00187 <font class="comment">/* place the two next digits in lowpart */</font>
00188 t1 = X_HW[1];
00189 lowpart = (t1 << SCS_NB_BITS) + X_HW[2];
00190
00191 <font class="comment">/* test for s/qNan, +/- Inf, +/- 0, placed here for obscure performance reasons */</font>
00192 <font class="keywordflow">if</font> (X_EXP != 1){
00193 *result = X_EXP;
00194 <font class="keywordflow">return</font>;
00195 }
00196
00197 <font class="comment">/* take the exponent of nb.d (will be in [0:SCS_NB_BITS])*/</font>
00198 exp = ((nb.i[HI_ENDIAN] & 0x7ff00000)>>20) - 1023;
00199 not_null = ((lowpart << (64+52 - 2*SCS_NB_BITS - exp)) != 0 );
00200 <font class="comment">/* Test if we are not on an exact double precision number */</font>
00201 <font class="keywordflow">for</font> (i=3; i<SCS_NB_WORDS; i++)
00202 <font class="keywordflow">if</font> (X_HW[i]!=0) not_null = 1;
00203
00204 <font class="comment">/* compute the exponent of the result */</font>
00205 expfinal = exp + SCS_NB_BITS*X_IND;
00206
00207 <font class="comment">/* Is the SCS number not too large for the IEEE exponent range ? */</font>
00208 <font class="keywordflow">if</font> (expfinal > 1023) {
00209 <font class="keywordflow">if</font> (rndMantissaUp)
00210 <font class="comment">/* return an infinity */</font>
00211 res = SCS_RADIX_RNG_DOUBLE*SCS_RADIX_RNG_DOUBLE;
00212 <font class="keywordflow">else</font>
00213 <font class="comment">/* infinity, rounded down, is SCS_MAX_DOUBLE */</font>
00214 res = SCS_MAX_DOUBLE;
00215 }
00216
00217 <font class="comment">/* Is our SCS number a denormal ? */</font>
00218 <font class="keywordflow">else</font> <font class="keywordflow">if</font> (expfinal >= -1022){
00219 <font class="comment">/* x is in the normal range */</font>
00220
00221 <font class="comment">/* align the rest of the mantissa to nb : shift by (2*SCS_NB_BITS)-53-exp */</font>
00222 lowpart = lowpart >> (exp+(2*SCS_NB_BITS)-52);
00223 <font class="comment">/* Finish to fill the mantissa */</font>
00224 nb.l = nb.l | lowpart;
00225 <font class="keywordflow">if</font> (rndMantissaUp && (not_null)){
00226 rndcorr.i[LO_ENDIAN] = 0;
00227 rndcorr.i[HI_ENDIAN] = (exp-52+1023)<<20; <font class="comment">/* 2^(exp-52) */</font>
00228 } <font class="keywordflow">else</font> {
00229 rndcorr.d = 0.0;
00230 }
00231 res = nb.d + rndcorr.d; <font class="comment">/* rounded to nearest */</font>
00232
00233 <font class="comment">/* now compute the exponent from the index :</font>
00234 <font class="comment"> we need to multiply res by 2^(X_IND*SCS_NB_BITS)</font>
00235 <font class="comment"> First check this number won't be a denormal itself */</font>
00236 <font class="keywordflow">if</font>((X_IND)*SCS_NB_BITS +1023 > 0) {
00237 <font class="comment">/* build the double 2^(X_IND*SCS_NB_BITS) */</font>
00238 nb.i[HI_ENDIAN] = ((X_IND)*SCS_NB_BITS +1023) << 20;
00239 nb.i[LO_ENDIAN] = 0;
00240 res *= nb.d; <font class="comment">/* exact multiplication */</font>
00241 }
00242 <font class="keywordflow">else</font> { <font class="comment">/*offset the previous computation by 2^(2*SCS_NB_BITS) */</font>
00243 <font class="comment">/* build the double 2^(X_IND*SCS_NB_BITS) */</font>
00244 nb.i[HI_ENDIAN] = ((X_IND)*SCS_NB_BITS +1023 + 2*SCS_NB_BITS) << 20;
00245 nb.i[LO_ENDIAN] = 0;
00246 res *= SCS_RADIX_MTWO_DOUBLE; <font class="comment">/* exact multiplication */</font>
00247 res *= nb.d; <font class="comment">/* exact multiplication */</font>
00248 }
00249 }
00250
00251
00252 <font class="keywordflow">else</font> {
00253 <font class="comment">/* the final number is a denormal with 52-(expfinal+1022)</font>
00254 <font class="comment"> significant bits. */</font>
00255
00256 <font class="keywordflow">if</font> (expfinal < -1022 - 53 ) {
00257 <font class="keywordflow">if</font>(rndMantissaUp)
00258 res = SCS_MIN_DOUBLE;
00259 <font class="keywordflow">else</font>
00260 res = 0.0;
00261 }
00262 <font class="keywordflow">else</font> {
00263
00264 <font class="comment">/* align the rest of the mantissa to nb */</font>
00265 lowpart = lowpart >> (exp+(2*SCS_NB_BITS)-52);
00266 <font class="comment">/* Finish to fill the mantissa */</font>
00267 nb.l = nb.l | lowpart;
00268
00269 <font class="comment">/* this is still a normal number. </font>
00270 <font class="comment"> Now remove its exponent and add back the implicit one */</font>
00271 nb.l = (nb.l & 0x000FFFFFFFFFFFFF) | 0x0010000000000000;
00272
00273 <font class="keywordflow">if</font> (rndMantissaUp && (not_null)){
00274 nb.l = nb.l >> (-1022 - expfinal);
00275 nb.l = nb.l +1; <font class="comment">/* works even if we move back into the normals*/</font>
00276 }
00277 <font class="keywordflow">else</font>
00278 <font class="comment">/* keep only the significant bits */</font>
00279 nb.l = nb.l >> (-1022 - expfinal);
00280
00281 res = nb.d;
00282
00283 <font class="comment">/* the exponent field is already set to zero so that's all */</font>
00284 }
00285 }
00286
00287 <font class="comment">/* sign management */</font>
00288 <font class="keywordflow">if</font> (X_SGN < 0)
00289 *result = - res;
00290 <font class="keywordflow">else</font>
00291 *result = res;
00292 }
00293
00294 <font class="preprocessor">#if 0</font>
00295 <font class="preprocessor"></font><font class="keywordtype">void</font> get_d_directed0(<font class="keywordtype">double</font> *result, <a class="code" href="structscs.html">scs_ptr</a> x,<font class="keywordtype">int</font> rndMantissaUp)
00296 {
00297 <font class="keywordtype">unsigned</font> <font class="keywordtype">long</font> <font class="keywordtype">long</font> <font class="keywordtype">int</font> lowpart, t1;
00298 scs_db_number nb, rndcorr;
00299 <font class="keywordtype">int</font> i, exp, not_null;
00300 <font class="keywordtype">double</font> res;
00301 <font class="comment">/* convert the MSB digit into a double, and store it in nb.d */</font>
00302 nb.d = (double)X_HW[0];
00303 <font class="comment">/* place the two next digits in lowpart */</font>
00304 t1 = X_HW[1];
00305 lowpart = (t1 << SCS_NB_BITS) + X_HW[2];
00306 <font class="comment">/* s/qNan, +/- Inf, +/- 0 */</font>
00307 <font class="keywordflow">if</font> (X_EXP != 1){
00308 *result = X_EXP;
00309 <font class="keywordflow">return</font>;
00310 }
00311 <font class="comment">/* take the exponent */</font>
00312 exp = ((nb.i[HI_ENDIAN] & 0x7ff00000)>>20) - 1023;
00313 not_null = ((lowpart << (64+52 - 2*SCS_NB_BITS - exp)) != 0 );
00314 <font class="comment">/* align the rest of the mantissa */</font>
00315 lowpart = lowpart >> (exp + 2*SCS_NB_BITS - 52);
00316 <font class="comment">/* Finish to fill the mantissa */</font>
00317 nb.l = nb.l | lowpart;
00318 <font class="comment">/* Test if we are not on an exact double precision number */</font>
00319 <font class="keywordflow">for</font> (i=3; i<SCS_NB_WORDS; i++)
00320 <font class="keywordflow">if</font> (X_HW[i]!=0) not_null = 1;
00321 <font class="keywordflow">if</font> (rndMantissaUp && (not_null)){
00322 rndcorr.i[LO_ENDIAN] = 0;
00323 rndcorr.i[HI_ENDIAN] = (exp-52+1023)<<20; <font class="comment">/* 2^(exp-52) */</font>
00324 } <font class="keywordflow">else</font> {
00325 rndcorr.d = 0.0;
00326 }
00327 res = nb.d + rndcorr.d; <font class="comment">/* make a rounded to nearest */</font>
00328 <font class="keywordflow">if</font> ((X_IND < SCS_MAX_RANGE) && (X_IND > -SCS_MAX_RANGE)){
00329 <font class="comment">/* x is comfortably in the double-precision range */</font>
00330 <font class="comment">/* build the double 2^(X_IND*SCS_NB_BITS) */</font>
00331 nb.i[HI_ENDIAN] = ((X_IND)*SCS_NB_BITS +1023) << 20;
00332 nb.i[LO_ENDIAN] = 0;
00333 res *= nb.d;
00334 }<font class="keywordflow">else</font> {
00335 <font class="comment">/* x may end up being a denormal or overflow */</font>
00336 i = X_IND;
00337 nb.d = 0;
00338 <font class="keywordflow">if</font> (X_IND > 0){
00339 <font class="comment">/* one of the following computations may lead to an overflow */</font>
00340 res *=SCS_RADIX_RNG_DOUBLE; <font class="comment">/* 2^(SCS_NB_BITS.SCS_MAX_RANGE) */</font>
00341 i -= SCS_MAX_RANGE;
00342 <font class="keywordflow">while</font>((i-->0)&&(res <= SCS_MAX_DOUBLE)) {
00343 <font class="comment">/* second test means: This loop stops on overflow */</font>
00344 res *= SCS_RADIX_ONE_DOUBLE;
00345 }
00346 }<font class="keywordflow">else</font> {
00347 <font class="comment">/* One of the computations may lead to denormal/underflow */</font>
00348 res *=SCS_RADIX_MRNG_DOUBLE; <font class="comment">/* 2^-(SCS_NB_BITS.SCS_MAX_RANGE)*/</font>
00349 i += SCS_MAX_RANGE;
00350 <font class="keywordflow">while</font>((i++<0)&&(res != 0)) {
00351 res *=SCS_RADIX_MONE_DOUBLE;
00352 }
00353 }
00354 }
00355 <font class="comment">/* sign management */</font>
00356 <font class="keywordflow">if</font> (X_SGN < 0)
00357 *result = - res;
00358 <font class="keywordflow">else</font>
00359 *result = res;
00360 }
00361
00362 <font class="preprocessor">#endif</font>
00363 <font class="preprocessor"></font><font class="comment">/*</font>
00364 <font class="comment"> * Rounded toward -Inf</font>
00365 <font class="comment"> */</font>
<a name="l00366"></a><a class="code" href="scs2double_8c.html#a5">00366</a> <font class="keywordtype">void</font> <a class="code" href="scs2double_8c.html#a5">scs_get_d_minf</a>(<font class="keywordtype">double</font> *result, <a class="code" href="structscs.html">scs_ptr</a> x){
00367
00368 <font class="comment">/* round up the mantissa if negative */</font>
00369 get_d_directed(result, x, (X_SGN<0));
00370 }
00371
00372
00373
00374 <font class="comment">/*</font>
00375 <font class="comment"> * Rounded toward +Inf</font>
00376 <font class="comment"> */</font>
<a name="l00377"></a><a class="code" href="scs2double_8c.html#a6">00377</a> <font class="keywordtype">void</font> <a class="code" href="scs2double_8c.html#a6">scs_get_d_pinf</a>(<font class="keywordtype">double</font> *result, <a class="code" href="structscs.html">scs_ptr</a> x){
00378
00379 <font class="comment">/* round up the mantissa if positive */</font>
00380 get_d_directed(result, x, (X_SGN>=0));
00381 }
00382
00383
00384
00385 <font class="comment">/*</font>
00386 <font class="comment"> * Rounded toward zero</font>
00387 <font class="comment"> */</font>
<a name="l00388"></a><a class="code" href="scs2double_8c.html#a7">00388</a> <font class="keywordtype">void</font> <a class="code" href="scs2double_8c.html#a7">scs_get_d_zero</a>(<font class="keywordtype">double</font> *result, <a class="code" href="structscs.html">scs_ptr</a> x){
00389 <font class="comment">/* never round up the mantissa */</font>
00390 get_d_directed(result, x, 0);
00391 }
</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>
