NFFT 3.5.3alpha
int.c
1/*
2 * Copyright (c) 2002, 2017 Jens Keiner, Stefan Kunis, Daniel Potts
3 *
4 * This program is free software; you can redistribute it and/or modify it under
5 * the terms of the GNU General Public License as published by the Free Software
6 * Foundation; either version 2 of the License, or (at your option) any later
7 * version.
8 *
9 * This program is distributed in the hope that it will be useful, but WITHOUT
10 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
11 * FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
12 * details.
13 *
14 * You should have received a copy of the GNU General Public License along with
15 * this program; if not, write to the Free Software Foundation, Inc., 51
16 * Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
17 */
18
19#include "infft.h"
20
21INT Y(exp2i)(const INT a)
22{
23 return (1U << a);
24}
25
32INT Y(log2i)(const INT m)
33{
34 /* Special case, zero or negative input. */
35 if (m <= 0)
36 return -1;
37
38#if SIZEOF_PTRDIFF_T == 8
39 /* Hash map with return values based on De Bruijn sequence. */
40 static INT debruijn[64] =
41 {
42 0, 58, 1, 59, 47, 53, 2, 60, 39, 48, 27, 54, 33, 42, 3, 61, 51, 37, 40, 49,
43 18, 28, 20, 55, 30, 34, 11, 43, 14, 22, 4, 62, 57, 46, 52, 38, 26, 32, 41,
44 50, 36, 17, 19, 29, 10, 13, 21, 56, 45, 25, 31, 35, 16, 9, 12, 44, 24, 15,
45 8, 23, 7, 6, 5, 63
46 };
47
48 register uint64_t v = (uint64_t)(m);
49
50 /* Round down to one less than a power of 2. */
51 v |= v >> 1;
52 v |= v >> 2;
53 v |= v >> 4;
54 v |= v >> 8;
55 v |= v >> 16;
56 v |= v >> 32;
57
58 /* 0x03f6eaf2cd271461 is a hexadecimal representation of a De Bruijn
59 * sequence for binary words of length 6. The binary representation
60 * starts with 000000111111. This is required to make it work with one less
61 * than a power of 2 instead of an actual power of 2.
62 */
63 return debruijn[(uint64_t)(v * 0x03f6eaf2cd271461LU) >> 58];
64#elif SIZEOF_PTRDIFF_T == 4
65 /* Hash map with return values based on De Bruijn sequence. */
66 static INT debruijn[32] =
67 {
68 0, 9, 1, 10, 13, 21, 2, 29, 11, 14, 16, 18, 22, 25, 3, 30, 8, 12, 20, 28,
69 15, 17, 24, 7, 19, 27, 23, 6, 26, 5, 4, 31
70 };
71
72 register uint32_t v = (uint32_t)(m);
73
74 /* Round down to one less than a power of 2. */
75 v |= v >> 1;
76 v |= v >> 2;
77 v |= v >> 4;
78 v |= v >> 8;
79 v |= v >> 16;
80
81 /* 0x07C4ACDD is a hexadecimal representation of a De Bruijn sequence for
82 * binary words of length 5. The binary representation starts with
83 * 0000011111. This is required to make it work with one less than a power of
84 * 2 instead of an actual power of 2.
85 */
86 return debruijn[(uint32_t)(v * 0x07C4ACDDU) >> 27];
87#else
88#error Incompatible size of ptrdiff_t
89#endif
90}
91
101INT Y(next_power_of_2)(const INT x)
102{
103 if (x < 0)
104 return -1;
105 else if (x < 2)
106 return x + 1;
107 else
108 {
109 uint64_t v = (uint64_t)x;
110
111 /* Subtract one, so we return the input if it is a power of two. */
112 v--;
113
114 /* Round down to one less than a power of two. */
115 v |= v >> 1;
116 v |= v >> 2;
117 v |= v >> 4;
118#if SIZEOF_PTRDIFF_T >= 2
119 v |= v >> 8;
120#endif
121#if SIZEOF_PTRDIFF_T >= 4
122 v |= v >> 16;
123#endif
124#if SIZEOF_PTRDIFF_T >= 8
125 v |= v >> 32;
126#endif
127 /* Add one to get power of two. */
128 v++;
129
130 return v;
131 }
132}
133
136void Y(next_power_of_2_exp)(const INT N, INT *N2, INT *t)
137{
138 INT n,i,logn;
139 INT N_is_not_power_of_2=0;
140
141 if (N == 0)
142 {
143 *N2 = 1;
144 *t = 0;
145 }
146 else
147 {
148 n = N;
149 logn = 0;
150 while (n != 1)
151 {
152 if (n%2 == 1)
153 {
154 N_is_not_power_of_2=1;
155 }
156 n = n/2;
157 logn++;
158 }
159
160 if (!N_is_not_power_of_2)
161 {
162 logn--;
163 }
164
165 for (i = 0; i <= logn; i++)
166 {
167 n = n*2;
168 }
169
170 *N2 = n;
171 *t = logn+1;
172 }
173}
174
175void Y(next_power_of_2_exp_int)(const int N, int *N2, int *t)
176{
177 int n,i,logn;
178 int N_is_not_power_of_2=0;
179
180 if (N == 0)
181 {
182 *N2 = 1;
183 *t = 0;
184 }
185 else
186 {
187 n = N;
188 logn = 0;
189 while (n != 1)
190 {
191 if (n%2 == 1)
192 {
193 N_is_not_power_of_2=1;
194 }
195 n = n/2;
196 logn++;
197 }
198
199 if (!N_is_not_power_of_2)
200 {
201 logn--;
202 }
203
204 for (i = 0; i <= logn; i++)
205 {
206 n = n*2;
207 }
208
209 *N2 = n;
210 *t = logn+1;
211 }
212}
Internal header file for auxiliary definitions and functions.