Main Page   Data Structures   File List   Data Fields   Globals  

scs2doubleCLEAN_BUT_SLOW.c

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

Generated on Tue Jun 17 10:15:51 2003 for SCSLib by doxygen1.2.15