Advertisement
Guest User

Untitled

a guest
Apr 19th, 2014
63
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
text 2.54 KB | None | 0 0
  1. typedef union {
  2. double f;
  3. struct {
  4. unsigned long mantisa : 52;
  5. unsigned long exponent : 11;
  6. unsigned long sign : 1;
  7. } parts;
  8. } double_cast;
  9.  
  10. double myfrexp(double number, int *exp)
  11. {
  12. double_cast d1;
  13. d1.f = number;
  14. unsigned long dd;
  15. printf("n %x n", d1.parts.exponent);
  16. *exp = d1.parts.exponent - 1022;
  17. printf("n%dnn", *exp);
  18. printf("n %lf n", (double)d1.parts.mantisa);
  19. return d1.parts.mantisa;
  20. }
  21.  
  22. double myfrexp(double number, int *exp) {
  23. static const uint64_t mantissa_mask = 0x000FFFFFFFFFFFFFllu;
  24. static const uint64_t mantissa_impliedBit = 0x0010000000000000llu;
  25. static const uint64_t expo_mask = 0x7FF0000000000000llu;
  26. static const uint64_t expo_norm = 0x3FE0000000000000llu;
  27. static const uint64_t sign_mask = 0x8000000000000000llu;
  28. static const int expo_NaN = 0x07FF;
  29. static const int expo_Bias = 1023;
  30.  
  31. union {
  32. double d;
  33. uint64_t u;
  34. } x = { number };
  35. uint64_t mantissa = x.u & mantissa_mask;
  36. int expo = (x.u & expo_mask) >> 52;
  37.  
  38. if (expo == expo_NaN) { // Behavior for Infinity and NaN is unspecified.
  39. *exp = 0;
  40. return number;
  41. }
  42. if (expo > 0) {
  43. mantissa |= mantissa_impliedBit; // This line is illustrative, not needed.
  44. expo -= expo_Bias;
  45. }
  46. else if (mantissa == 0) {
  47. *exp = 0;
  48. return number; // Do not return 0.0 as that does not preserve -0.0
  49. }
  50. else {
  51. // de-normal or sub-normal numbers
  52. expo = 1 - expo_Bias; // Bias different when biased exponent is 0
  53. while (mantissa < mantissa_impliedBit) {
  54. mantissa <<= 1;
  55. expo--;
  56. }
  57. }
  58. *exp = expo + 1;
  59. mantissa &= ~mantissa_impliedBit;
  60. x.u = (x.u & sign_mask) | expo_norm | mantissa;
  61. return x.d;
  62. }
  63.  
  64. #include <limits.h>
  65. #include <math.h>
  66. #include <memory.h>
  67. #include <stdio.h>
  68. #include <float.h>
  69.  
  70. void frexp_test(double d) {
  71. int i1,i2;
  72. double d1,d2;
  73. d1 = frexp(d, &i1);
  74. d2 = myfrexp(d, &i2);
  75. if (memcmp(&d1,&d2,sizeof(d1)) != 0 || (i1 != i2)) {
  76. printf("%a (%a %x) (%a %x)n", d, d1, i1, d2, i2);
  77. }
  78. }
  79.  
  80. int main() {
  81. frexp_test(1.0);
  82. frexp_test(0.0);
  83. frexp_test(-0.0);
  84. frexp_test(DBL_MAX);
  85. frexp_test(-DBL_MAX);
  86. frexp_test(DBL_EPSILON);
  87. frexp_test(DBL_MIN);
  88. frexp_test(DBL_MIN/1024);
  89. frexp_test(DBL_MIN/1024/1024);
  90. frexp_test(INFINITY);
  91. //frexp_test(DBL_TRUE_MIN);
  92. return 0;
  93. }
  94.  
  95. uint64_t bits = *((uint64_t*)&number);
  96. uint64_t mantissa = bits & ((1<<52)-1);
  97. unsigned exponent = (bits>>52) & 0x7FF;
  98. unsigned sign = bits>>63;
  99.  
  100. printf("%c1.%"PRIu64"e%u", sign?'-':'+', mantissa, exponent);
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement