maint: update copyright
[gnulib.git] / lib / integer_length.c
1 /* integer_length - find most significant bit in an 'unsigned int'.
2    Copyright (C) 2011-2014 Free Software Foundation, Inc.
3
4    This program is free software: you can redistribute it and/or modify
5    it under the terms of the GNU General Public License as published by
6    the Free Software Foundation; either version 3 of the License, or
7    (at your option) any later version.
8
9    This program is distributed in the hope that it will be useful,
10    but WITHOUT ANY WARRANTY; without even the implied warranty of
11    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
12    GNU General Public License for more details.
13
14    You should have received a copy of the GNU General Public License
15    along with this program.  If not, see <http://www.gnu.org/licenses/>.  */
16
17 /* Written by Bruno Haible <bruno@clisp.org>, 2011.  */
18
19 #include <config.h>
20
21 /* Specification.  */
22 #include "integer_length.h"
23
24 #include <limits.h>
25
26 #include "float+.h"
27
28 /* MSVC with option -fp:strict refuses to compile constant initializers that
29    contain floating-point operations.  Pacify this compiler.  */
30 #ifdef _MSC_VER
31 # pragma fenv_access (off)
32 #endif
33
34 #define NBITS (sizeof (unsigned int) * CHAR_BIT)
35
36 int
37 integer_length (unsigned int x)
38 {
39 #if __GNUC__ > 3 || (__GNUC__ == 3 && __GNUC_MINOR__ >= 4)
40   if (x == 0)
41     return 0;
42   else
43     return NBITS - __builtin_clz (x);
44 #else
45 # if defined DBL_EXPBIT0_WORD && defined DBL_EXPBIT0_BIT
46   if (NBITS <= DBL_MANT_BIT)
47     {
48       /* Use 'double' operations.
49          Assumes an IEEE 754 'double' implementation.  */
50 #  define DBL_EXP_MASK ((DBL_MAX_EXP - DBL_MIN_EXP) | 7)
51 #  define DBL_EXP_BIAS (DBL_EXP_MASK / 2 - 1)
52 #  define NWORDS \
53     ((sizeof (double) + sizeof (unsigned int) - 1) / sizeof (unsigned int))
54       typedef union { double value; unsigned int word[NWORDS]; }
55               memory_double;
56
57       if (x == 0)
58         return 0;
59       else
60         {
61           memory_double m;
62           unsigned int exponent;
63
64           if (1)
65             {
66               /* Use a single integer to floating-point conversion.  */
67               m.value = x;
68             }
69           else
70             {
71               /* Use a single floating-point subtraction.  */
72               /* 2^(DBL_MANT_DIG-1).  */
73               static const double TWO_DBL_MANT_DIG =
74                 /* Assume DBL_MANT_DIG <= 5 * 31.
75                    Use the identity
76                    n = floor(n/5) + floor((n+1)/5) + ... + floor((n+4)/5).  */
77                 (double) (1U << ((DBL_MANT_DIG - 1) / 5))
78                 * (double) (1U << ((DBL_MANT_DIG - 1 + 1) / 5))
79                 * (double) (1U << ((DBL_MANT_DIG - 1 + 2) / 5))
80                 * (double) (1U << ((DBL_MANT_DIG - 1 + 3) / 5))
81                 * (double) (1U << ((DBL_MANT_DIG - 1 + 4) / 5));
82
83               /* Construct 2^(DBL_MANT_DIG-1) + x by hand.  */
84               m.word[DBL_EXPBIT0_WORD] =
85                 (DBL_MANT_DIG + DBL_EXP_BIAS) << DBL_EXPBIT0_BIT;
86               m.word[1 - DBL_EXPBIT0_WORD] = x;
87
88               /* Subtract 2^(DBL_MANT_DIG-1).  */
89               m.value = m.value - TWO_DBL_MANT_DIG;
90             }
91
92           exponent =
93             (m.word[DBL_EXPBIT0_WORD] >> DBL_EXPBIT0_BIT) & DBL_EXP_MASK;
94           return exponent - DBL_EXP_BIAS;
95         }
96     }
97   else
98 # endif
99     if (NBITS == 32)
100       {
101         /* 6 comparisons.  */
102         if (x != 0)
103           {
104             int result = 1;
105             if (x >= 0x10000)
106               {
107                 x = x >> 16;
108                 result += 16;
109               }
110             if (x >= 0x100)
111               {
112                 x = x >> 8;
113                 result += 8;
114               }
115             if (x >= 0x10)
116               {
117                 x = x >> 4;
118                 result += 4;
119               }
120             if (x >= 0x4)
121               {
122                 x = x >> 2;
123                 result += 2;
124               }
125             if (x >= 0x2)
126               {
127                 x = x >> 1;
128                 result += 1;
129               }
130             return result;
131           }
132         else
133           return 0;
134       }
135     else
136       {
137         /* Naive loop.
138            Works for any value of NBITS.  */
139         int j;
140
141         for (j = NBITS - 1; j >= 0; j--)
142           if (x & (1U << j))
143             return j + 1;
144         return 0;
145       }
146 #endif
147 }