1 1.1 jtc /* @(#)s_logb.c 5.1 93/09/24 */ 2 1.1 jtc /* 3 1.1 jtc * ==================================================== 4 1.1 jtc * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. 5 1.1 jtc * 6 1.1 jtc * Developed at SunPro, a Sun Microsystems, Inc. business. 7 1.1 jtc * Permission to use, copy, modify, and distribute this 8 1.1 jtc * software is freely granted, provided that this notice 9 1.1 jtc * is preserved. 10 1.1 jtc * ==================================================== 11 1.1 jtc */ 12 1.3 jtc 13 1.3 jtc #ifndef lint 14 1.4 jtc static char rcsid[] = "$Id: s_logb.c,v 1.4 1994/03/03 17:04:41 jtc Exp $"; 15 1.3 jtc #endif 16 1.1 jtc 17 1.1 jtc /* 18 1.1 jtc * double logb(x) 19 1.1 jtc * IEEE 754 logb. Included to pass IEEE test suite. Not recommend. 20 1.1 jtc * Use ilogb instead. 21 1.1 jtc */ 22 1.1 jtc 23 1.2 jtc #include <math.h> 24 1.4 jtc #include <machine/endian.h> 25 1.1 jtc 26 1.4 jtc #if BYTE_ORDER == LITTLE_ENDIAN 27 1.4 jtc #define n0 1 28 1.1 jtc #else 29 1.4 jtc #define n0 0 30 1.1 jtc #endif 31 1.1 jtc 32 1.1 jtc #ifdef __STDC__ 33 1.1 jtc double logb(double x) 34 1.1 jtc #else 35 1.1 jtc double logb(x) 36 1.1 jtc double x; 37 1.1 jtc #endif 38 1.1 jtc { 39 1.4 jtc int lx,ix; 40 1.1 jtc ix = (*(n0+(int*)&x))&0x7fffffff; /* high |x| */ 41 1.1 jtc lx = *(1-n0+(int*)&x); /* low x */ 42 1.1 jtc if((ix|lx)==0) return -1.0/fabs(x); 43 1.1 jtc if(ix>=0x7ff00000) return x*x; 44 1.1 jtc if((ix>>=20)==0) /* IEEE 754 logb */ 45 1.1 jtc return -1022.0; 46 1.1 jtc else 47 1.1 jtc return (double) (ix-1023); 48 1.1 jtc } 49