From 137b9334f3909bb9333d50892492062745a2991c Mon Sep 17 00:00:00 2001 From: Alan Cox Date: Mon, 5 Mar 2018 20:27:47 +0000 Subject: [PATCH] frexp: replace MUSL with Sun version for portability --- Library/libs/frexp.c | 62 ++++++++++++++++++++++++++++++-------------- 1 file changed, 42 insertions(+), 20 deletions(-) diff --git a/Library/libs/frexp.c b/Library/libs/frexp.c index 509e92f9..4023ce4b 100644 --- a/Library/libs/frexp.c +++ b/Library/libs/frexp.c @@ -1,25 +1,47 @@ -/* From MUSL */ - #include -#include +#include "libm.h" -double frexp(double x, int *e) -{ - union { double d; uint64_t i; } y = { x }; - int ee = y.i>>52 & 0x7ff; +/* + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunPro, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + */ - if (!ee) { - if (x) { - x = frexp(x*0x1p64, e); - *e -= 64; - } else *e = 0; - return x; - } else if (ee == 0x7ff) { - return x; - } - *e = ee - 0x3fe; - y.i &= 0x800fffffffffffffull; - y.i |= 0x3fe0000000000000ull; - return y.d; +/* + * for non-zero x + * x = frexp(arg,&exp); + * return a double fp quantity x such that 0.5 <= |x| <1.0 + * and the corresponding binary exponent "exp". That is + * arg = x*2^exp. + * If arg is inf, 0.0, or NaN, then frexp(arg,&exp) returns arg + * with *exp=0. + */ + +static const double +two54 = 1.80143985094819840000e+16; /* 0x43500000, 0x00000000 */ + +double +frexp(double x, int *eptr) +{ + int32_t hx, ix, lx; + EXTRACT_WORDS(hx,lx,x); + ix = 0x7fffffff&hx; + *eptr = 0; + if(ix>=0x7ff00000||((ix|lx)==0)) return x; /* 0,inf,nan */ + if (ix<0x00100000) { /* subnormal */ + x *= two54; + GET_HIGH_WORD(hx,x); + ix = hx&0x7fffffff; + *eptr = -54; + } + *eptr += (ix>>20)-1022; + hx = (hx&0x800fffff)|0x3fe00000; + SET_HIGH_WORD(x,hx); + return x; } -- 2.34.1