diff options
author | Nicholas Noll <nbnoll@eml.cc> | 2021-12-05 16:16:21 -0800 |
---|---|---|
committer | Nicholas Noll <nbnoll@eml.cc> | 2021-12-05 16:16:21 -0800 |
commit | 07e77936d535e58b0aeb4f2a11400c1050556739 (patch) | |
tree | fb50fad6436ecbc159505bee73a92cd289b99a07 /src/base/math/frexp.c | |
parent | b48327d357e0818d1a6ae2a064cfa7d1567e1242 (diff) |
Feat: added math library
Used Plan9's libc as starting point. This cleans up dangling references
due to loss of libc.
Diffstat (limited to 'src/base/math/frexp.c')
-rw-r--r-- | src/base/math/frexp.c | 120 |
1 files changed, 120 insertions, 0 deletions
diff --git a/src/base/math/frexp.c b/src/base/math/frexp.c new file mode 100644 index 0000000..7d119bb --- /dev/null +++ b/src/base/math/frexp.c @@ -0,0 +1,120 @@ +#include <u.h> +#include <base.h> + +/* + * 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 +math·frexp(double d, int *ep) +{ + sys·DoubleWord x, a; + + *ep = 0; + /* order matters: only isNaN can operate on NaN */ + if(math·isNaN(d) || math·isInf(d, 0) || d == 0) + return d; + x.x = d; + a.x = math·fabs(d); + if((a.hi >> SHIFT) == 0){ /* normalize subnormal numbers */ + x.x = (double)(1ULL<<SIG) * d; + *ep = -SIG; + } + *ep += ((x.hi >> SHIFT) & MASK) - BIAS; + x.hi &= ~(MASK << SHIFT); + x.hi |= BIAS << SHIFT; + return x.x; +} + +double +math·ldexp(double d, int deltae) +{ + int e, bits; + sys·DoubleWord 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 math·Inf(-1); + return math·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<<SHIFT)-1); + x.hi &= ~((1<<SHIFT)-1); + x.hi |= (1<<(SHIFT-1)) | (z>>1); + } + while(deltae > 0){ /* shift bits down */ + bits = deltae; + if(bits > SHIFT) + bits = SHIFT; + x.lo >>= bits; + x.lo |= (x.hi&((1<<bits)-1)) << (32-bits); + z = x.hi & ((1<<SHIFT)-1); + x.hi &= ~((1<<SHIFT)-1); + x.hi |= z>>bits; + deltae -= bits; + } + } + x.hi &= ~(MASK << SHIFT); + x.hi |= (long)e << SHIFT; + return x.x; +} + +double +math·modf(double d, double *ip) +{ + sys·DoubleWord 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 = math·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; +} |