1 #include <math.h> 2 #include <errno.h> 3 #define _RESEARCH_SOURCE 4 #include <float.h> 5 6 #define MASK 0x7ffL 7 #define SHIFT 20 8 #define BIAS 1022L 9 #define SIG 52 10 11 typedef union 12 { 13 double d; 14 struct 15 { 16 #ifdef IEEE_8087 17 long ls; 18 long ms; 19 #else 20 long ms; 21 long ls; 22 #endif 23 }; 24 } Cheat; 25 26 double frexp(double d,int * ep)27frexp(double d, int *ep) 28 { 29 Cheat x, a; 30 31 *ep = 0; 32 /* order matters: only isNaN can operate on NaN */ 33 if(isNaN(d) || isInf(d, 0) || d == 0) 34 return d; 35 x.d = d; 36 a.d = fabs(d); 37 if((a.ms >> SHIFT) == 0){ /* normalize subnormal numbers */ 38 x.d = (double)(1ULL<<SIG) * d; 39 *ep = -SIG; 40 } 41 *ep = ((x.ms >> SHIFT) & MASK) - BIAS; 42 x.ms &= ~(MASK << SHIFT); 43 x.ms |= BIAS << SHIFT; 44 return x.d; 45 } 46 47 double ldexp(double d,int e)48ldexp(double d, int e) 49 { 50 Cheat x; 51 52 if(d == 0) 53 return 0; 54 x.d = d; 55 e += (x.ms >> SHIFT) & MASK; 56 if(e <= 0) 57 return 0; 58 if(e >= MASK){ 59 errno = ERANGE; 60 if(d < 0) 61 return -HUGE_VAL; 62 return HUGE_VAL; 63 } 64 x.ms &= ~(MASK << SHIFT); 65 x.ms |= (long)e << SHIFT; 66 return x.d; 67 } 68 69 double modf(double d,double * ip)70modf(double d, double *ip) 71 { 72 double f; 73 Cheat x; 74 int e; 75 76 x.d = d; 77 e = (x.ms >> SHIFT) & MASK; 78 if(e == MASK){ 79 *ip = d; 80 if(x.ls != 0 || (x.ms & 0xfffffL) != 0) /* NaN */ 81 return d; 82 /* ±Inf */ 83 x.ms &= 0x80000000L; 84 return x.d; 85 } 86 if(d < 1) { 87 if(d < 0) { 88 f = modf(-d, ip); 89 *ip = -*ip; 90 return -f; 91 } 92 *ip = 0; 93 return d; 94 } 95 e -= BIAS; 96 if(e <= SHIFT+1) { 97 x.ms &= ~(0x1fffffL >> e); 98 x.ls = 0; 99 } else 100 if(e <= SHIFT+33) 101 x.ls &= ~(0x7fffffffL >> (e-SHIFT-2)); 102 *ip = x.d; 103 return d - x.d; 104 } 105