Main Page   Data Structures   File List   Data Fields   Globals  

scs2double.c

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

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