overby.me
this repo has no description
1
2/*
3 * logb.c
4 *
5 * by Ian Ollmann
6 *
7 * Copyright (c) 2007, Apple Inc. All Rights Reserved.
8 *
9 * Unix 2003 implementation of logb()
10 */
11
12#include <math.h>
13#include <stdint.h>
14
15#ifdef ARMLIBM_SET_FLAGS
16
17#include "required_arithmetic.h"
18
19double logb( double x )
20{
21 union{ double d; uint64_t u;}u = {x};
22
23 u.u &= 0x7fffffffffffffffULL;
24 int32_t exp = (int32_t) (u.u >> 52);
25
26 if( __builtin_expect( (uint32_t) exp - 1U >= 2046, 0 ) )
27 { // +-0, +-denorm, +-inf, NaN
28 if( 0.0 == x )
29 return required_divide_double( -1.0, 0.0 ); //set div/0, return -Inf
30
31 if( x != x )
32 return x + x;
33
34 if( u.u == 0x7ff0000000000000ULL )
35 return __builtin_inf();
36
37 u.u |= 0x3ff0000000000000ULL;
38 u.d -= 1.0;
39 exp = (int32_t) (u.u >> 52);
40
41 return (double) (exp - (1023+1022));
42 }
43
44 return (double) (exp - 1023);
45}
46
47#else
48
49double logb( double x )
50{
51 union{ double d; uint64_t u;}u = {x};
52
53 u.u &= 0x7fffffffffffffffULL;
54 int32_t exp = (int32_t) (u.u >> 52);
55
56 if( __builtin_expect( (uint32_t) exp - 1U >= 2046, 0 ) )
57 { // +-0, +-denorm, +-inf, NaN
58 if( 0.0 == x )
59 return -__builtin_inf();
60
61 if( x != x )
62 return x + x;
63
64 if( u.u == 0x7ff0000000000000ULL )
65 return __builtin_inf();
66
67 u.u |= 0x3ff0000000000000ULL;
68 u.d -= 1.0;
69 exp = (int32_t) (u.u >> 52);
70
71 return (double) (exp - (1023+1022));
72 }
73
74 return (double) (exp - 1023);
75}
76
77#endif