1/*	$OpenBSD: s_logbl.c,v 1.1 2008/12/09 20:00:35 martynas Exp $	*/
2/*
3 * From: @(#)s_ilogb.c 5.1 93/09/24
4 * ====================================================
5 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
6 *
7 * Developed at SunPro, a Sun Microsystems, Inc. business.
8 * Permission to use, copy, modify, and distribute this
9 * software is freely granted, provided that this notice
10 * is preserved.
11 * ====================================================
12 */
13
14#include <sys/types.h>
15#include <machine/ieee.h>
16#include <float.h>
17#include <limits.h>
18#include <math.h>
19
20long double
21logbl(long double x)
22{
23	union {
24		long double e;
25		struct ieee_ext bits;
26	} u;
27	unsigned long m;
28	int b;
29
30	u.e = x;
31	if (u.bits.ext_exp == 0) {
32		if ((u.bits.ext_fracl
33#ifdef EXT_FRACLMBITS
34			| u.bits.ext_fraclm
35#endif /* EXT_FRACLMBITS */
36#ifdef EXT_FRACHMBITS
37			| u.bits.ext_frachm
38#endif /* EXT_FRACHMBITS */
39			| u.bits.ext_frach) == 0) {	/* x == 0 */
40			u.bits.ext_sign = 1;
41			return (1.0L / u.e);
42		}
43		/* denormalized */
44		if (u.bits.ext_frach == 0
45#ifdef EXT_FRACHMBITS
46			&& u.bits.ext_frachm == 0
47#endif
48			) {
49			m = 1lu << (EXT_FRACLBITS - 1);
50			for (b = EXT_FRACHBITS; !(u.bits.ext_fracl & m); m >>= 1)
51				b++;
52#if defined(EXT_FRACHMBITS) && defined(EXT_FRACLMBITS)
53			m = 1lu << (EXT_FRACLMBITS - 1);
54			for (b += EXT_FRACHMBITS; !(u.bits.ext_fraclm & m);
55				m >>= 1)
56				b++;
57#endif /* defined(EXT_FRACHMBITS) && defined(EXT_FRACLMBITS) */
58		} else {
59			m = 1lu << (EXT_FRACHBITS - 1);
60			for (b = 0; !(u.bits.ext_frach & m); m >>= 1)
61				b++;
62#ifdef EXT_FRACHMBITS
63			m = 1lu << (EXT_FRACHMBITS - 1);
64			for (; !(u.bits.ext_frachm & m); m >>= 1)
65				b++;
66#endif /* EXT_FRACHMBITS */
67		}
68#ifdef EXT_IMPLICIT_NBIT
69		b++;
70#endif
71		return ((long double)(LDBL_MIN_EXP - b - 1));
72	}
73	if (u.bits.ext_exp < (LDBL_MAX_EXP << 1) - 1)	/* normal */
74		return ((long double)(u.bits.ext_exp - LDBL_MAX_EXP + 1));
75	else						/* +/- inf or nan */
76		return (x * x);
77}
78