12116Sjkh/* @(#)s_logb.c 5.1 93/09/24 */
22116Sjkh/*
32116Sjkh * ====================================================
42116Sjkh * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
52116Sjkh *
62116Sjkh * Developed at SunPro, a Sun Microsystems, Inc. business.
72116Sjkh * Permission to use, copy, modify, and distribute this
88870Srgrimes * software is freely granted, provided that this notice
92116Sjkh * is preserved.
102116Sjkh * ====================================================
112116Sjkh */
122116Sjkh
13176101Sbde#include <sys/cdefs.h>
14176101Sbde__FBSDID("$FreeBSD: releng/11.0/lib/msun/src/s_logb.c 176101 2008-02-08 01:22:13Z bde $");
152116Sjkh
162116Sjkh/*
172116Sjkh * double logb(x)
182116Sjkh * IEEE 754 logb. Included to pass IEEE test suite. Not recommend.
192116Sjkh * Use ilogb instead.
202116Sjkh */
212116Sjkh
22174698Sdas#include <float.h>
23174698Sdas
242116Sjkh#include "math.h"
252116Sjkh#include "math_private.h"
262116Sjkh
27153049Sbdestatic const double
28153049Sbdetwo54 = 1.80143985094819840000e+16;	/* 43500000 00000000 */
29153049Sbde
3097413Salfreddouble
31117912Speterlogb(double x)
322116Sjkh{
332116Sjkh	int32_t lx,ix;
342116Sjkh	EXTRACT_WORDS(ix,lx,x);
352116Sjkh	ix &= 0x7fffffff;			/* high |x| */
362116Sjkh	if((ix|lx)==0) return -1.0/fabs(x);
372116Sjkh	if(ix>=0x7ff00000) return x*x;
38153049Sbde	if(ix<0x00100000) {
39153049Sbde		x *= two54;		 /* convert subnormal x to normal */
40176101Sbde		GET_HIGH_WORD(ix,x);
41153049Sbde		ix &= 0x7fffffff;
42176101Sbde		return (double) ((ix>>20)-1023-54);
43153049Sbde	} else
44153049Sbde		return (double) ((ix>>20)-1023);
452116Sjkh}
46174698Sdas
47174698Sdas#if (LDBL_MANT_DIG == 53)
48174698Sdas__weak_reference(logb, logbl);
49174698Sdas#endif
50