00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028 #include "scs.h"
00029 #include "scs_private.h"
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
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
00056
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
00062 i = X_IND;
00063 nb.d = 0;
00064 if (X_IND > 0){
00065
00066 *res *=SCS_RADIX_RNG_DOUBLE;
00067 i -= SCS_MAX_RANGE;
00068 while((i-->0)&&(*res <= MAX_DOUBLE)) {
00069
00070 *res *= SCS_RADIX_ONE_DOUBLE;
00071 }
00072 }else {
00073
00074 *res *=SCS_RADIX_MRNG_DOUBLE;
00075 i += SCS_MAX_RANGE;
00076 while((i++<0)&&(*res != 0)) {
00077 *res *=SCS_RADIX_MONE_DOUBLE;
00078 }
00079 }
00080 }
00081
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
00094 nb.d = (double)X_HW[0];
00095
00096
00097 t1 = X_HW[1];
00098 lowpart = (t1 << SCS_NB_BITS) + X_HW[2];
00099
00100
00101 if (X_EXP != 1){*result = X_EXP; return;}
00102
00103
00104 exp = ((nb.i[HI_ENDIAN] & 0x7ff00000)>>20) - 1023;
00105
00106 lowpart = lowpart >> (exp+(2*SCS_NB_BITS)-53);
00107
00108
00109 if (lowpart & 0x0000000000000001){
00110
00111 rndcorr.i[LO_ENDIAN] = 0;
00112 rndcorr.i[HI_ENDIAN] = (exp+971) << 20;
00113 }else{
00114
00115 rndcorr.l = 0;
00116 }
00117
00118 lowpart = lowpart >> 1;
00119 nb.l = nb.l | lowpart;
00120 *result = nb.d + rndcorr.d;
00121
00122
00123
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
00143 if (X_EXP != 1){
00144 *res = X_EXP;
00145 return;
00146 }
00147
00148
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
00154 lowpart = lowpart >> (exp + 2*SCS_NB_BITS - 52);
00155
00156
00157 nb.l = nb.l | lowpart;
00158
00159
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
00165 rndcorr.i[LO_ENDIAN] = 0;
00166 rndcorr.i[HI_ENDIAN] = (exp+971)<<20;
00167 } else {
00168 rndcorr.d = 0.0;
00169 }
00170
00171 *res = nb.d + rndcorr.d;
00172 computeExponent(x, res);
00173 }
00174
00175
00176
00177
00178
00179 void scs_get_d_minf(double *result, scs_ptr x){
00180
00181
00182 scs_get_d_directed(x, (X_SGN<0), result);
00183
00184 return;
00185 }
00186
00187
00188
00189
00190
00191
00192 void scs_get_d_pinf(double *result, scs_ptr x){
00193
00194
00195 scs_get_d_directed(x, (X_SGN>=0), result);
00196
00197 return;
00198 }
00199
00200
00201
00202
00203
00204
00205 void scs_get_d_zero(double *result, scs_ptr x){
00206
00207 scs_get_d_directed(x, 0, result);
00208 }