Main Page   Data Structures   File List   Data Fields   Globals  

double2scs.c

Go to the documentation of this file.
00001 /** Conversion of floating-point double to SCS
00002 @file double2scs.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 double precision number in it SCS multiprecision
00033   representation
00034  */
00035 
00036 void scs_set_d(scs_ptr result, double x){
00037   scs_db_number nb, mantissa;
00038   int exponent, exponent_remainder;
00039   int ind, i;
00040 
00041   if(x>=0){R_SGN = 1;    nb.d = x;}
00042   else    {R_SGN = -1;   nb.d = -x;}
00043 
00044   exponent = nb.i[HI_ENDIAN] & 0x7ff00000 ;
00045 
00046   if (exponent == 0x7ff00000)  {
00047     /*
00048      * x = +/- Inf, s/qNAN
00049      */
00050     R_EXP = x;  
00051     for(i=0; i<SCS_NB_WORDS; i++)
00052       R_HW[i] = 0;
00053  
00054     R_IND = 0;
00055     R_SGN = 1;
00056   } 
00057 
00058   else {    /* Normals,  denormals, +/- 0.  */
00059 
00060     /* This number is not an exception */
00061     R_EXP = 1;
00062 
00063 #if 1
00064 
00065     if (exponent == 0){
00066       /* x is a denormal number : bring it back to the normal range */
00067       nb.d = nb.d * SCS_RADIX_TWO_DOUBLE;      /* 2^(2.SCS_NB_BITS) */
00068       exponent = nb.i[HI_ENDIAN] & 0x7ff00000 ;
00069       R_IND = -2;
00070     }else {
00071       R_IND = 0;
00072     }
00073 
00074     exponent = exponent >> 20;  /* get the actual value */
00075 
00076     ind = ((exponent +((100*SCS_NB_BITS)-1023))/SCS_NB_BITS) - 100 ;  
00077     /* means : = (exponent -1023 + 100*SCS_NB_BITS)/SCS_NB_BITS -100 
00078      The business with 100*SCS_NB_BITS is to stay within the positive
00079      range for exponent_remainder between 1 and SCS_NB_BITS */
00080 
00081     exponent_remainder = exponent - 1022 - (SCS_NB_BITS*ind);
00082 
00083     R_IND += ind;
00084 
00085     /* now get the mantissa and add the implicit 1 in fp. format*/
00086     mantissa.l = (nb.l & 0x000fffffffffffff) | 0x0010000000000000;
00087 
00088 
00089     /* and spread it over the structure
00090        Everything here is 64-bit arithmetic */
00091     R_HW[0] = mantissa.l >> (53 - exponent_remainder);
00092 
00093     /* 11 = 64-53 */
00094     mantissa.l =  (mantissa.l << (exponent_remainder+11));
00095     R_HW[1] = (mantissa.i[HI_ENDIAN] >> (32 - SCS_NB_BITS))& SCS_MASK_RADIX ;
00096     mantissa.l =  (mantissa.l << SCS_NB_BITS);
00097     R_HW[2] = (mantissa.i[HI_ENDIAN] >> (32 - SCS_NB_BITS))& SCS_MASK_RADIX ;
00098 #if SCS_NB_BITS < 27
00099     mantissa.l =  (mantissa.l << SCS_NB_BITS);
00100     R_HW[3] = (mantissa.i[HI_ENDIAN] >> (32 - SCS_NB_BITS))& SCS_MASK_RADIX ;
00101 #else
00102     R_HW[3] = 0 ;
00103 #endif
00104 
00105 #if (SCS_NB_WORDS==8)
00106       R_HW[4] = 0; R_HW[5] = 0; R_HW[6] = 0; R_HW[7] = 0;
00107 #else
00108     for(i=4; i<SCS_NB_WORDS; i++)
00109       R_HW[i] = 0;
00110 #endif
00111 
00112 
00113 
00114 #else /* Other algorithm as in the research report. Slower */
00115     R_IND = 0;
00116 
00117     while(nb.d>SCS_RADIX_ONE_DOUBLE) {
00118         R_IND++;
00119         nb.d *= SCS_RADIX_MONE_DOUBLE;
00120       }
00121 
00122     while(nb.d<1)  {
00123         R_IND--;
00124         nb.d *= SCS_RADIX_ONE_DOUBLE;
00125       }
00126 
00127     i=0;
00128     while(nb.d != 0){
00129       R_HW[i] = (unsigned int) nb.d;
00130       nb.d = (nb.d - (double)R_HW[i]) * SCS_RADIX_ONE_DOUBLE;
00131       i++;
00132     }
00133     for(; i<SCS_NB_WORDS; i++)
00134       R_HW[i] = 0;
00135 
00136 #endif
00137 
00138   } /* end if test NaN etc */
00139 
00140   return;
00141 }
00142 
00143 
00144 /**
00145   Convert an integer number in it scs multiprecision
00146   representation
00147  */
00148 void scs_set_si(scs_ptr result, signed int x){
00149   unsigned int ux;
00150   int i;
00151   
00152   if(x>=0){R_SGN = 1;   ux = x;}
00153   else    {R_SGN = -1;  ux = -x;}
00154   
00155 
00156   if (ux > SCS_RADIX){
00157     R_IND   = 1;
00158     R_HW[0] = (ux - SCS_RADIX) >> SCS_NB_BITS;
00159     R_HW[1] =  ux - (R_HW[0] << SCS_NB_BITS);
00160   }else {
00161     R_IND   = 0;
00162     R_HW[0] = ux;
00163     R_HW[1] = 0;
00164   }
00165 
00166   for(i=2; i<SCS_NB_WORDS; i++)
00167     R_HW[i] = 0;
00168 
00169   if (x != 0)  R_EXP = 1;
00170   else         R_EXP = 0;
00171   
00172   return;
00173 }
00174 

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