1 /* @(#)s_frexp.c 5.1 93/09/24 */ 2 /* 3 * ==================================================== 4 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. 5 * 6 * Developed at SunPro, a Sun Microsystems, Inc. business. 7 * Permission to use, copy, modify, and distribute this 8 * software is freely granted, provided that this notice 9 * is preserved. 10 * ==================================================== 11 */ 12 13 #ifndef lint 14 static char rcsid[] = "$Id: s_frexp.c,v 1.4 1994/03/03 17:04:36 jtc Exp $"; 15 #endif 16 17 /* 18 * for non-zero x 19 * x = frexp(arg,&exp); 20 * return a double fp quantity x such that 0.5 <= |x| <1.0 21 * and the corresponding binary exponent "exp". That is 22 * arg = x*2^exp. 23 * If arg is inf, 0.0, or NaN, then frexp(arg,&exp) returns arg 24 * with *exp=0. 25 */ 26 27 #include <math.h> 28 #include <machine/endian.h> 29 30 #if BYTE_ORDER == LITTLE_ENDIAN 31 #define n0 1 32 #else 33 #define n0 0 34 #endif 35 36 #ifdef __STDC__ 37 static const double 38 #else 39 static double 40 #endif 41 one = 1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */ 42 two54 = 1.80143985094819840000e+16; /* 0x43500000, 0x00000000 */ 43 44 #ifdef __STDC__ 45 double frexp(double x, int *eptr) 46 #else 47 double frexp(x, eptr) 48 double x; int *eptr; 49 #endif 50 { 51 int hx, ix, lx; 52 hx = *(n0+(int*)&x); 53 ix = 0x7fffffff&hx; 54 lx = *(1-n0+(int*)&x); 55 *eptr = 0; 56 if(ix>=0x7ff00000||((ix|lx)==0)) return x; /* 0,inf,nan */ 57 if (ix<0x00100000) { /* subnormal */ 58 x *= two54; 59 hx = *(n0+(int*)&x); 60 ix = hx&0x7fffffff; 61 *eptr = -54; 62 } 63 *eptr += (ix>>20)-1022; 64 hx = (hx&0x800fffff)|0x3fe00000; 65 *(n0 + (int*)&x) = hx; 66 return x; 67 } 68