#include #include /* * this is big/little endian non-portable * it gets the endian from the FPdbleword * union in u.h. */ #define MASK 0x7ffL #define SHIFT 20 #define BIAS 1022L #define SIG 52 double frexp(double d, int *ep) { FPdbleword x, a; *ep = 0; /* order matters: only isNaN can operate on NaN */ if(isNaN(d) || isInf(d, 0) || d == 0) return d; x.x = d; a.x = fabs(d); if((a.hi >> SHIFT) == 0){ /* normalize subnormal numbers */ x.x = (double)(1ULL<> SHIFT) & MASK) - BIAS; x.hi &= ~(MASK << SHIFT); x.hi |= BIAS << SHIFT; return x.x; } double ldexp(double d, int deltae) { int e, bits; FPdbleword x; ulong z; if(d == 0) return 0; x.x = d; e = (x.hi >> SHIFT) & MASK; if(deltae >= 0 || e+deltae >= 1){ /* no underflow */ e += deltae; if(e >= MASK){ /* overflow */ if(d < 0) return Inf(-1); return Inf(1); } }else{ /* underflow gracefully */ deltae = -deltae; /* need to shift d right deltae */ if(e > 1){ /* shift e-1 by exponent manipulation */ deltae -= e-1; e = 1; } if(deltae > 0 && e==1){ /* shift 1 by switch from 1.fff to 0.1ff */ deltae--; e = 0; x.lo >>= 1; x.lo |= (x.hi&1)<<31; z = x.hi & ((1<>1); } while(deltae > 0){ /* shift bits down */ bits = deltae; if(bits > SHIFT) bits = SHIFT; x.lo >>= bits; x.lo |= (x.hi&((1<>bits; deltae -= bits; } } x.hi &= ~(MASK << SHIFT); x.hi |= (long)e << SHIFT; return x.x; } double modf(double d, double *ip) { FPdbleword x; int e; x.x = d; e = (x.hi >> SHIFT) & MASK; if(e == MASK){ *ip = d; if(x.lo != 0 || (x.hi & 0xfffffL) != 0) /* NaN */ return d; /* ±Inf */ x.hi &= 0x80000000L; return x.x; } if(d < 1) { if(d < 0) { x.x = modf(-d, ip); *ip = -*ip; return -x.x; } *ip = 0; return d; } e -= BIAS; if(e <= SHIFT+1) { x.hi &= ~(0x1fffffL >> e); x.lo = 0; } else if(e <= SHIFT+33) x.lo &= ~(0x7fffffffL >> (e-SHIFT-2)); *ip = x.x; return d - x.x; }