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 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
00057 nb.d = (double)X_HW[0];
00058
00059
00060 t1 = X_HW[1];
00061 lowpart = (t1 << SCS_NB_BITS) + X_HW[2];
00062
00063
00064
00065
00066
00067
00068
00069 if (X_EXP != 1){
00070 *result = X_EXP;
00071 return;
00072 }
00073
00074
00075 exp = ((nb.i[HI_ENDIAN] & 0x7ff00000)>>20) - 1023;
00076
00077
00078 expfinal = exp + SCS_NB_BITS*X_IND;
00079
00080
00081 if (expfinal > 1023) {
00082
00083 res = SCS_RADIX_RNG_DOUBLE*SCS_RADIX_RNG_DOUBLE;
00084 }
00085
00086
00087 else if (expfinal >= -1022){
00088
00089
00090
00091 lowpart = lowpart >> (exp+(2*SCS_NB_BITS)-53);
00092
00093 if (lowpart & 0x0000000000000001){
00094
00095 rndcorr.i[LO_ENDIAN] = 0;
00096 rndcorr.i[HI_ENDIAN] = (exp-52+1023)<<20;
00097 }else{
00098
00099 rndcorr.d = 0.0;
00100 }
00101
00102 lowpart = lowpart >> 1;
00103 nb.l = nb.l | lowpart;
00104 res = nb.d + rndcorr.d;
00105
00106
00107
00108
00109 if((X_IND)*SCS_NB_BITS +1023 > 0) {
00110
00111 nb.i[HI_ENDIAN] = ((X_IND)*SCS_NB_BITS +1023) << 20;
00112 nb.i[LO_ENDIAN] = 0;
00113 res *= nb.d;
00114 }
00115 else {
00116
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;
00120 res *= nb.d;
00121 }
00122 }
00123
00124
00125 else {
00126
00127
00128
00129 if (expfinal < -1022 - 53 ) {
00130 res = 0.0;
00131 }
00132 else {
00133
00134
00135 lowpart = lowpart >> (exp+(2*SCS_NB_BITS)-52);
00136
00137 nb.l = nb.l | lowpart;
00138
00139
00140
00141 nb.l = (nb.l & 0x000FFFFFFFFFFFFF) | 0x0010000000000000;
00142
00143
00144 nb.l = nb.l >> (-1023 - expfinal);
00145
00146 if (nb.i[LO_ENDIAN] & 0x00000001){
00147
00148 rndcorr.l = 1;
00149 }else{
00150
00151 rndcorr.d = 0.0;
00152
00153 }
00154 res = 0.5*(nb.d + rndcorr.d);
00155
00156
00157 }
00158 }
00159
00160
00161 if (X_SGN < 0)
00162 *result = - res;
00163 else
00164 *result = res;
00165 }
00166
00167
00168
00169
00170
00171
00172
00173
00174
00175
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
00185 nb.d = (double)X_HW[0];
00186
00187
00188 t1 = X_HW[1];
00189 lowpart = (t1 << SCS_NB_BITS) + X_HW[2];
00190
00191
00192 if (X_EXP != 1){
00193 *result = X_EXP;
00194 return;
00195 }
00196
00197
00198 exp = ((nb.i[HI_ENDIAN] & 0x7ff00000)>>20) - 1023;
00199 not_null = ((lowpart << (64+52 - 2*SCS_NB_BITS - exp)) != 0 );
00200
00201 for (i=3; i<SCS_NB_WORDS; i++)
00202 if (X_HW[i]!=0) not_null = 1;
00203
00204
00205 expfinal = exp + SCS_NB_BITS*X_IND;
00206
00207
00208 if (expfinal > 1023) {
00209 if (rndMantissaUp)
00210
00211 res = SCS_RADIX_RNG_DOUBLE*SCS_RADIX_RNG_DOUBLE;
00212 else
00213
00214 res = SCS_MAX_DOUBLE;
00215 }
00216
00217
00218 else if (expfinal >= -1022){
00219
00220
00221
00222 lowpart = lowpart >> (exp+(2*SCS_NB_BITS)-52);
00223
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;
00228 } else {
00229 rndcorr.d = 0.0;
00230 }
00231 res = nb.d + rndcorr.d;
00232
00233
00234
00235
00236 if((X_IND)*SCS_NB_BITS +1023 > 0) {
00237
00238 nb.i[HI_ENDIAN] = ((X_IND)*SCS_NB_BITS +1023) << 20;
00239 nb.i[LO_ENDIAN] = 0;
00240 res *= nb.d;
00241 }
00242 else {
00243
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;
00247 res *= nb.d;
00248 }
00249 }
00250
00251
00252 else {
00253
00254
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
00265 lowpart = lowpart >> (exp+(2*SCS_NB_BITS)-52);
00266
00267 nb.l = nb.l | lowpart;
00268
00269
00270
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;
00276 }
00277 else
00278
00279 nb.l = nb.l >> (-1022 - expfinal);
00280
00281 res = nb.d;
00282
00283
00284 }
00285 }
00286
00287
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
00302 nb.d = (double)X_HW[0];
00303
00304 t1 = X_HW[1];
00305 lowpart = (t1 << SCS_NB_BITS) + X_HW[2];
00306
00307 if (X_EXP != 1){
00308 *result = X_EXP;
00309 return;
00310 }
00311
00312 exp = ((nb.i[HI_ENDIAN] & 0x7ff00000)>>20) - 1023;
00313 not_null = ((lowpart << (64+52 - 2*SCS_NB_BITS - exp)) != 0 );
00314
00315 lowpart = lowpart >> (exp + 2*SCS_NB_BITS - 52);
00316
00317 nb.l = nb.l | lowpart;
00318
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;
00324 } else {
00325 rndcorr.d = 0.0;
00326 }
00327 res = nb.d + rndcorr.d;
00328 if ((X_IND < SCS_MAX_RANGE) && (X_IND > -SCS_MAX_RANGE)){
00329
00330
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
00336 i = X_IND;
00337 nb.d = 0;
00338 if (X_IND > 0){
00339
00340 res *=SCS_RADIX_RNG_DOUBLE;
00341 i -= SCS_MAX_RANGE;
00342 while((i-->0)&&(res <= SCS_MAX_DOUBLE)) {
00343
00344 res *= SCS_RADIX_ONE_DOUBLE;
00345 }
00346 }else {
00347
00348 res *=SCS_RADIX_MRNG_DOUBLE;
00349 i += SCS_MAX_RANGE;
00350 while((i++<0)&&(res != 0)) {
00351 res *=SCS_RADIX_MONE_DOUBLE;
00352 }
00353 }
00354 }
00355
00356 if (X_SGN < 0)
00357 *result = - res;
00358 else
00359 *result = res;
00360 }
00361
00362 #endif
00363
00364
00365
00366 void scs_get_d_minf(double *result, scs_ptr x){
00367
00368
00369 get_d_directed(result, x, (X_SGN<0));
00370 }
00371
00372
00373
00374
00375
00376
00377 void scs_get_d_pinf(double *result, scs_ptr x){
00378
00379
00380 get_d_directed(result, x, (X_SGN>=0));
00381 }
00382
00383
00384
00385
00386
00387
00388 void scs_get_d_zero(double *result, scs_ptr x){
00389
00390 get_d_directed(result, x, 0);
00391 }