1 #include <u.h> 2 #include <libc.h> 3 4 /* 5 * this is big/little endian non-portable 6 * it gets the endian from the FPdbleword 7 * union in u.h. 8 */ 9 #define MASK 0x7ffL 10 #define SHIFT 20 11 #define BIAS 1022L 12 13 double 14 frexp(double d, int *ep) 15 { 16 FPdbleword x; 17 18 if(d == 0) { 19 *ep = 0; 20 return 0; 21 } 22 x.x = d; 23 *ep = ((x.hi >> SHIFT) & MASK) - BIAS; 24 x.hi &= ~(MASK << SHIFT); 25 x.hi |= BIAS << SHIFT; 26 return x.x; 27 } 28 29 double 30 ldexp(double d, int deltae) 31 { 32 int e, bits; 33 FPdbleword x; 34 ulong z; 35 36 if(d == 0) 37 return 0; 38 x.x = d; 39 e = (x.hi >> SHIFT) & MASK; 40 if(deltae >= 0 || e+deltae >= 1){ /* no underflow */ 41 e += deltae; 42 if(e >= MASK){ /* overflow */ 43 if(d < 0) 44 return Inf(-1); 45 return Inf(1); 46 } 47 }else{ /* underflow gracefully */ 48 deltae = -deltae; 49 /* need to shift d right deltae */ 50 if(e > 1){ /* shift e-1 by exponent manipulation */ 51 deltae -= e-1; 52 e = 1; 53 } 54 if(deltae > 0 && e==1){ /* shift 1 by switch from 1.fff to 0.1ff */ 55 deltae--; 56 e = 0; 57 x.lo >>= 1; 58 x.lo |= (x.hi&1)<<31; 59 z = x.hi & ((1<<SHIFT)-1); 60 x.hi &= ~((1<<SHIFT)-1); 61 x.hi |= (1<<(SHIFT-1)) | (z>>1); 62 } 63 while(deltae > 0){ /* shift bits down */ 64 bits = deltae; 65 if(bits > SHIFT) 66 bits = SHIFT; 67 x.lo >>= bits; 68 x.lo |= (x.hi&((1<<bits)-1)) << (32-bits); 69 z = x.hi & ((1<<SHIFT)-1); 70 x.hi &= ~((1<<SHIFT)-1); 71 x.hi |= z>>bits; 72 deltae -= bits; 73 } 74 } 75 x.hi &= ~(MASK << SHIFT); 76 x.hi |= (long)e << SHIFT; 77 return x.x; 78 } 79 80 double 81 modf(double d, double *ip) 82 { 83 FPdbleword x; 84 int e; 85 86 if(d < 1) { 87 if(d < 0) { 88 x.x = modf(-d, ip); 89 *ip = -*ip; 90 return -x.x; 91 } 92 *ip = 0; 93 return d; 94 } 95 x.x = d; 96 e = ((x.hi >> SHIFT) & MASK) - BIAS; 97 if(e <= SHIFT+1) { 98 x.hi &= ~(0x1fffffL >> e); 99 x.lo = 0; 100 } else 101 if(e <= SHIFT+33) 102 x.lo &= ~(0x7fffffffL >> (e-SHIFT-2)); 103 *ip = x.x; 104 return d - x.x; 105 } 106