1 | // |
---|
2 | // Cforall Version 1.0.0 Copyright (C) 2016 University of Waterloo |
---|
3 | // |
---|
4 | // The contents of this file are covered under the licence agreement in the |
---|
5 | // file "LICENCE" distributed with Cforall. |
---|
6 | // |
---|
7 | // math.hfa -- |
---|
8 | // |
---|
9 | // Author : Peter A. Buhr |
---|
10 | // Created On : Mon Apr 18 23:37:04 2016 |
---|
11 | // Last Modified By : Peter A. Buhr |
---|
12 | // Last Modified On : Sat Jun 17 22:59:42 2023 |
---|
13 | // Update Count : 199 |
---|
14 | // |
---|
15 | |
---|
16 | #pragma once |
---|
17 | |
---|
18 | #include <math.h> |
---|
19 | #include <complex.h> |
---|
20 | |
---|
21 | //--------------------------------------- |
---|
22 | |
---|
23 | #include "common.hfa" |
---|
24 | #include "bits/debug.hfa" |
---|
25 | |
---|
26 | //---------------------- General ---------------------- |
---|
27 | |
---|
28 | static inline __attribute__((always_inline)) { |
---|
29 | float ?%?( float x, float y ) { return fmodf( x, y ); } |
---|
30 | float fmod( float x, float y ) { return fmodf( x, y ); } |
---|
31 | double ?%?( double x, double y ) { return fmod( x, y ); } |
---|
32 | // extern "C" { double fmod( double, double ); } |
---|
33 | long double ?%?( long double x, long double y ) { return fmodl( x, y ); } |
---|
34 | long double fmod( long double x, long double y ) { return fmodl( x, y ); } |
---|
35 | |
---|
36 | float remainder( float x, float y ) { return remainderf( x, y ); } |
---|
37 | // extern "C" { double remainder( double, double ); } |
---|
38 | long double remainder( long double x, long double y ) { return remainderl( x, y ); } |
---|
39 | |
---|
40 | float remquo( float x, float y, int * quo ) { return remquof( x, y, quo ); } |
---|
41 | // extern "C" { double remquo( double x, double y, int * quo ); } |
---|
42 | long double remquo( long double x, long double y, int * quo ) { return remquol( x, y, quo ); } |
---|
43 | [ int, float ] remquo( float x, float y ) { int quo; x = remquof( x, y, &quo ); return [ quo, x ]; } |
---|
44 | [ int, double ] remquo( double x, double y ) { int quo; x = remquo( x, y, &quo ); return [ quo, x ]; } |
---|
45 | [ int, long double ] remquo( long double x, long double y ) { int quo; x = remquol( x, y, &quo ); return [ quo, x ]; } |
---|
46 | |
---|
47 | [ float, float ] div( float x, float y ) { y = modff( x / y, &x ); return [ x, y ]; } |
---|
48 | [ double, double ] div( double x, double y ) { y = modf( x / y, &x ); return [ x, y ]; } |
---|
49 | [ long double, long double ] div( long double x, long double y ) { y = modfl( x / y, &x ); return [ x, y ]; } |
---|
50 | |
---|
51 | float fma( float x, float y, float z ) { return fmaf( x, y, z ); } |
---|
52 | // extern "C" { double fma( double, double, double ); } |
---|
53 | long double fma( long double x, long double y, long double z ) { return fmal( x, y, z ); } |
---|
54 | |
---|
55 | float fdim( float x, float y ) { return fdimf( x, y ); } |
---|
56 | // extern "C" { double fdim( double, double ); } |
---|
57 | long double fdim( long double x, long double y ) { return fdiml( x, y ); } |
---|
58 | |
---|
59 | float nan( const char tag[] ) { return nanf( tag ); } |
---|
60 | // extern "C" { double nan( const char [] ); } |
---|
61 | long double nan( const char tag[] ) { return nanl( tag ); } |
---|
62 | } // distribution |
---|
63 | |
---|
64 | //---------------------- Exponential ---------------------- |
---|
65 | |
---|
66 | static inline __attribute__((always_inline)) { |
---|
67 | float exp( float x ) { return expf( x ); } |
---|
68 | // extern "C" { double exp( double ); } |
---|
69 | long double exp( long double x ) { return expl( x ); } |
---|
70 | float _Complex exp( float _Complex x ) { return cexpf( x ); } |
---|
71 | double _Complex exp( double _Complex x ) { return cexp( x ); } |
---|
72 | long double _Complex exp( long double _Complex x ) { return cexpl( x ); } |
---|
73 | |
---|
74 | float exp2( float x ) { return exp2f( x ); } |
---|
75 | // extern "C" { double exp2( double ); } |
---|
76 | long double exp2( long double x ) { return exp2l( x ); } |
---|
77 | //float _Complex exp2( float _Complex x ) { return cexp2f( x ); } |
---|
78 | //double _Complex exp2( double _Complex x ) { return cexp2( x ); } |
---|
79 | //long double _Complex exp2( long double _Complex x ) { return cexp2l( x ); } |
---|
80 | |
---|
81 | float expm1( float x ) { return expm1f( x ); } |
---|
82 | // extern "C" { double expm1( double ); } |
---|
83 | long double expm1( long double x ) { return expm1l( x ); } |
---|
84 | |
---|
85 | float pow( float x, float y ) { return powf( x, y ); } |
---|
86 | // extern "C" { double pow( double, double ); } |
---|
87 | long double pow( long double x, long double y ) { return powl( x, y ); } |
---|
88 | float _Complex pow( float _Complex x, float _Complex y ) { return cpowf( x, y ); } |
---|
89 | double _Complex pow( double _Complex x, double _Complex y ) { return cpow( x, y ); } |
---|
90 | long double _Complex pow( long double _Complex x, long double _Complex y ) { return cpowl( x, y ); } |
---|
91 | } // distribution |
---|
92 | |
---|
93 | //---------------------- Logarithm ---------------------- |
---|
94 | |
---|
95 | static inline __attribute__((always_inline)) { |
---|
96 | float log( float x ) { return logf( x ); } |
---|
97 | // extern "C" { double log( double ); } |
---|
98 | long double log( long double x ) { return logl( x ); } |
---|
99 | float _Complex log( float _Complex x ) { return clogf( x ); } |
---|
100 | double _Complex log( double _Complex x ) { return clog( x ); } |
---|
101 | long double _Complex log( long double _Complex x ) { return clogl( x ); } |
---|
102 | |
---|
103 | // O(1) polymorphic integer log2, using clz, which returns the number of leading 0-bits, starting at the most |
---|
104 | // significant bit (single instruction on x86) |
---|
105 | int log2( unsigned int n ) { return n == 0 ? -1 : sizeof(n) * __CHAR_BIT__ - 1 - __builtin_clz( n ); } |
---|
106 | long int log2( unsigned long int n ) { return n == 0 ? -1 : sizeof(n) * __CHAR_BIT__ - 1 - __builtin_clzl( n ); } |
---|
107 | long long int log2( unsigned long long int n ) { return n == 0 ? -1 : sizeof(n) * __CHAR_BIT__ - 1 - __builtin_clzll( n ); } |
---|
108 | float log2( float x ) { return log2f( x ); } |
---|
109 | // extern "C" { double log2( double ); } |
---|
110 | long double log2( long double x ) { return log2l( x ); } |
---|
111 | // float _Complex log2( float _Complex x ) { return clog2f( x ); } |
---|
112 | // double _Complex log2( double _Complex x ) { return clog2( x ); } |
---|
113 | // long double _Complex log2( long double _Complex x ) { return clog2l( x ); } |
---|
114 | |
---|
115 | float log10( float x ) { return log10f( x ); } |
---|
116 | // extern "C" { double log10( double ); } |
---|
117 | long double log10( long double x ) { return log10l( x ); } |
---|
118 | // float _Complex log10( float _Complex x ) { return clog10f( x ); } |
---|
119 | // double _Complex log10( double _Complex x ) { return clog10( x ); } |
---|
120 | // long double _Complex log10( long double _Complex x ) { return clog10l( x ); } |
---|
121 | |
---|
122 | float log1p( float x ) { return log1pf( x ); } |
---|
123 | // extern "C" { double log1p( double ); } |
---|
124 | long double log1p( long double x ) { return log1pl( x ); } |
---|
125 | |
---|
126 | int ilogb( float x ) { return ilogbf( x ); } |
---|
127 | // extern "C" { int ilogb( double ); } |
---|
128 | int ilogb( long double x ) { return ilogbl( x ); } |
---|
129 | |
---|
130 | float logb( float x ) { return logbf( x ); } |
---|
131 | // extern "C" { double logb( double ); } |
---|
132 | long double logb( long double x ) { return logbl( x ); } |
---|
133 | |
---|
134 | float sqrt( float x ) { return sqrtf( x ); } |
---|
135 | // extern "C" { double sqrt( double ); } |
---|
136 | long double sqrt( long double x ) { return sqrtl( x ); } |
---|
137 | float _Complex sqrt( float _Complex x ) { return csqrtf( x ); } |
---|
138 | double _Complex sqrt( double _Complex x ) { return csqrt( x ); } |
---|
139 | long double _Complex sqrt( long double _Complex x ) { return csqrtl( x ); } |
---|
140 | |
---|
141 | float cbrt( float x ) { return cbrtf( x ); } |
---|
142 | // extern "C" { double cbrt( double ); } |
---|
143 | long double cbrt( long double x ) { return cbrtl( x ); } |
---|
144 | |
---|
145 | float hypot( float x, float y ) { return hypotf( x, y ); } |
---|
146 | // extern "C" { double hypot( double, double ); } |
---|
147 | long double hypot( long double x, long double y ) { return hypotl( x, y ); } |
---|
148 | } // distribution |
---|
149 | |
---|
150 | static inline unsigned long long log2_u32_32( unsigned long long val ) { |
---|
151 | enum { |
---|
152 | TABLE_BITS = 6, |
---|
153 | TABLE_SIZE = (1 << TABLE_BITS) + 2, |
---|
154 | }; |
---|
155 | // for(i; TABLE_SIZE) { |
---|
156 | // table[i] = (unsigned long long)(log2(1.0 + i / pow(2, TABLE_BITS)) * pow(2, 32))); |
---|
157 | // } |
---|
158 | static const unsigned long long table[] = { |
---|
159 | 0x0000000000, 0x0005b9e5a1, 0x000b5d69ba, 0x0010eb389f, |
---|
160 | 0x001663f6fa, 0x001bc84240, 0x002118b119, 0x002655d3c4, |
---|
161 | 0x002b803473, 0x00309857a0, 0x00359ebc5b, 0x003a93dc98, |
---|
162 | 0x003f782d72, 0x00444c1f6b, 0x0049101eac, 0x004dc4933a, |
---|
163 | 0x005269e12f, 0x00570068e7, 0x005b888736, 0x006002958c, |
---|
164 | 0x00646eea24, 0x0068cdd829, 0x006d1fafdc, 0x007164beb4, |
---|
165 | 0x00759d4f80, 0x0079c9aa87, 0x007dea15a3, 0x0081fed45c, |
---|
166 | 0x0086082806, 0x008a064fd5, 0x008df988f4, 0x0091e20ea1, |
---|
167 | 0x0095c01a39, 0x009993e355, 0x009d5d9fd5, 0x00a11d83f4, |
---|
168 | 0x00a4d3c25e, 0x00a8808c38, 0x00ac241134, 0x00afbe7fa0, |
---|
169 | 0x00b3500472, 0x00b6d8cb53, 0x00ba58feb2, 0x00bdd0c7c9, |
---|
170 | 0x00c1404ead, 0x00c4a7ba58, 0x00c80730b0, 0x00cb5ed695, |
---|
171 | 0x00ceaecfea, 0x00d1f73f9c, 0x00d53847ac, 0x00d8720935, |
---|
172 | 0x00dba4a47a, 0x00ded038e6, 0x00e1f4e517, 0x00e512c6e5, |
---|
173 | 0x00e829fb69, 0x00eb3a9f01, 0x00ee44cd59, 0x00f148a170, |
---|
174 | 0x00f446359b, 0x00f73da38d, 0x00fa2f045e, 0x00fd1a708b, |
---|
175 | 0x0100000000, 0x0102dfca16, |
---|
176 | }; |
---|
177 | _Static_assert((sizeof(table) / sizeof(table[0])) == TABLE_SIZE, "TABLE_SIZE should be accurate"); |
---|
178 | // starting from val = (2 ** i)*(1 + f) where 0 <= f < 1 |
---|
179 | // log identities mean log2(val) = log2((2 ** i)*(1 + f)) = log2(2**i) + log2(1+f) |
---|
180 | // |
---|
181 | // getting i is easy to do using builtin_clz (count leading zero) |
---|
182 | // |
---|
183 | // we want to calculate log2(1+f) independently to have a many bits of precision as possible. |
---|
184 | // val = (2 ** i)*(1 + f) = 2 ** i + f * 2 ** i |
---|
185 | // isolating f we get |
---|
186 | // val - 2 ** i = f * 2 ** i |
---|
187 | // (val - 2 ** i) / 2 ** i = f |
---|
188 | // |
---|
189 | // we want to interpolate from the table to get the values |
---|
190 | // and compromise by doing quadratic interpolation (rather than higher degree interpolation) |
---|
191 | // |
---|
192 | // for the interpolation we want to shift everything the fist sample point |
---|
193 | // so our parabola becomes x = 0 |
---|
194 | // this further simplifies the equations |
---|
195 | // |
---|
196 | // the consequence is that we need f in 2 forms: |
---|
197 | // - finding the index of x0 |
---|
198 | // - finding the distance between f and x0 |
---|
199 | // |
---|
200 | // since sample points are equidistant we can significantly simplify the equations |
---|
201 | |
---|
202 | // get i |
---|
203 | const unsigned long long bits = sizeof(val) * __CHAR_BIT__; |
---|
204 | const unsigned long long lz = __builtin_clzl(val); |
---|
205 | const unsigned long long i = bits - 1 - lz; |
---|
206 | |
---|
207 | // get the fractinal part as a u32.32 |
---|
208 | const unsigned long long frac = (val << (lz + 1)) >> 32; |
---|
209 | |
---|
210 | // get high order bits for the index into the table |
---|
211 | const unsigned long long idx0 = frac >> (32 - TABLE_BITS); |
---|
212 | |
---|
213 | // get the x offset, i.e., the difference between the first sample point and the actual fractional part |
---|
214 | const long long udx = frac - (idx0 << (32 - TABLE_BITS)); |
---|
215 | /* paranoid */ verify((idx0 + 2) < TABLE_SIZE); |
---|
216 | |
---|
217 | const long long y0 = table[idx0 + 0]; |
---|
218 | const long long y1 = table[idx0 + 1]; |
---|
219 | const long long y2 = table[idx0 + 2]; |
---|
220 | |
---|
221 | // from there we can quadraticly interpolate to get the data, using the lagrange polynomial |
---|
222 | // normally it would look like: |
---|
223 | // double r0 = y0 * ((x - x1) / (x0 - x1)) * ((x - x2) / (x0 - x2)); |
---|
224 | // double r1 = y1 * ((x - x0) / (x1 - x0)) * ((x - x2) / (x1 - x2)); |
---|
225 | // double r2 = y2 * ((x - x0) / (x2 - x0)) * ((x - x1) / (x2 - x1)); |
---|
226 | // but since the spacing between sample points is fixed, we can simplify it and extract common expressions |
---|
227 | const long long f1 = (y1 - y0); |
---|
228 | const long long f2 = (y2 - y0); |
---|
229 | const long long a = f2 - (f1 * 2l); |
---|
230 | const long long b = (f1 * 2l) - a; |
---|
231 | |
---|
232 | // Now we can compute it in the form (ax + b)x + c (which avoid repeating steps) |
---|
233 | long long sum = ((a*udx) >> (32 - TABLE_BITS)) + b; |
---|
234 | sum = (sum*udx) >> (32 - TABLE_BITS + 1); |
---|
235 | sum = y0 + sum; |
---|
236 | |
---|
237 | return (i << 32) + (sum); |
---|
238 | } // log2_u32_32 |
---|
239 | |
---|
240 | //---------------------- Trigonometric ---------------------- |
---|
241 | |
---|
242 | static inline __attribute__((always_inline)) { |
---|
243 | float sin( float x ) { return sinf( x ); } |
---|
244 | // extern "C" { double sin( double ); } |
---|
245 | long double sin( long double x ) { return sinl( x ); } |
---|
246 | float _Complex sin( float _Complex x ) { return csinf( x ); } |
---|
247 | double _Complex sin( double _Complex x ) { return csin( x ); } |
---|
248 | long double _Complex sin( long double _Complex x ) { return csinl( x ); } |
---|
249 | |
---|
250 | float cos( float x ) { return cosf( x ); } |
---|
251 | // extern "C" { double cos( double ); } |
---|
252 | long double cos( long double x ) { return cosl( x ); } |
---|
253 | float _Complex cos( float _Complex x ) { return ccosf( x ); } |
---|
254 | double _Complex cos( double _Complex x ) { return ccos( x ); } |
---|
255 | long double _Complex cos( long double _Complex x ) { return ccosl( x ); } |
---|
256 | |
---|
257 | float tan( float x ) { return tanf( x ); } |
---|
258 | // extern "C" { double tan( double ); } |
---|
259 | long double tan( long double x ) { return tanl( x ); } |
---|
260 | float _Complex tan( float _Complex x ) { return ctanf( x ); } |
---|
261 | double _Complex tan( double _Complex x ) { return ctan( x ); } |
---|
262 | long double _Complex tan( long double _Complex x ) { return ctanl( x ); } |
---|
263 | |
---|
264 | float asin( float x ) { return asinf( x ); } |
---|
265 | // extern "C" { double asin( double ); } |
---|
266 | long double asin( long double x ) { return asinl( x ); } |
---|
267 | float _Complex asin( float _Complex x ) { return casinf( x ); } |
---|
268 | double _Complex asin( double _Complex x ) { return casin( x ); } |
---|
269 | long double _Complex asin( long double _Complex x ) { return casinl( x ); } |
---|
270 | |
---|
271 | float acos( float x ) { return acosf( x ); } |
---|
272 | // extern "C" { double acos( double ); } |
---|
273 | long double acos( long double x ) { return acosl( x ); } |
---|
274 | float _Complex acos( float _Complex x ) { return cacosf( x ); } |
---|
275 | double _Complex acos( double _Complex x ) { return cacos( x ); } |
---|
276 | long double _Complex acos( long double _Complex x ) { return cacosl( x ); } |
---|
277 | |
---|
278 | float atan( float x ) { return atanf( x ); } |
---|
279 | // extern "C" { double atan( double ); } |
---|
280 | long double atan( long double x ) { return atanl( x ); } |
---|
281 | float _Complex atan( float _Complex x ) { return catanf( x ); } |
---|
282 | double _Complex atan( double _Complex x ) { return catan( x ); } |
---|
283 | long double _Complex atan( long double _Complex x ) { return catanl( x ); } |
---|
284 | |
---|
285 | float atan2( float x, float y ) { return atan2f( x, y ); } |
---|
286 | // extern "C" { double atan2( double, double ); } |
---|
287 | long double atan2( long double x, long double y ) { return atan2l( x, y ); } |
---|
288 | |
---|
289 | // alternative name for atan2 |
---|
290 | float atan( float x, float y ) { return atan2f( x, y ); } |
---|
291 | double atan( double x, double y ) { return atan2( x, y ); } |
---|
292 | long double atan( long double x, long double y ) { return atan2l( x, y ); } |
---|
293 | } // distribution |
---|
294 | |
---|
295 | //---------------------- Hyperbolic ---------------------- |
---|
296 | |
---|
297 | static inline __attribute__((always_inline)) { |
---|
298 | float sinh( float x ) { return sinhf( x ); } |
---|
299 | // extern "C" { double sinh( double ); } |
---|
300 | long double sinh( long double x ) { return sinhl( x ); } |
---|
301 | float _Complex sinh( float _Complex x ) { return csinhf( x ); } |
---|
302 | double _Complex sinh( double _Complex x ) { return csinh( x ); } |
---|
303 | long double _Complex sinh( long double _Complex x ) { return csinhl( x ); } |
---|
304 | |
---|
305 | float cosh( float x ) { return coshf( x ); } |
---|
306 | // extern "C" { double cosh( double ); } |
---|
307 | long double cosh( long double x ) { return coshl( x ); } |
---|
308 | float _Complex cosh( float _Complex x ) { return ccoshf( x ); } |
---|
309 | double _Complex cosh( double _Complex x ) { return ccosh( x ); } |
---|
310 | long double _Complex cosh( long double _Complex x ) { return ccoshl( x ); } |
---|
311 | |
---|
312 | float tanh( float x ) { return tanhf( x ); } |
---|
313 | // extern "C" { double tanh( double ); } |
---|
314 | long double tanh( long double x ) { return tanhl( x ); } |
---|
315 | float _Complex tanh( float _Complex x ) { return ctanhf( x ); } |
---|
316 | double _Complex tanh( double _Complex x ) { return ctanh( x ); } |
---|
317 | long double _Complex tanh( long double _Complex x ) { return ctanhl( x ); } |
---|
318 | |
---|
319 | float asinh( float x ) { return asinhf( x ); } |
---|
320 | // extern "C" { double asinh( double ); } |
---|
321 | long double asinh( long double x ) { return asinhl( x ); } |
---|
322 | float _Complex asinh( float _Complex x ) { return casinhf( x ); } |
---|
323 | double _Complex asinh( double _Complex x ) { return casinh( x ); } |
---|
324 | long double _Complex asinh( long double _Complex x ) { return casinhl( x ); } |
---|
325 | |
---|
326 | float acosh( float x ) { return acoshf( x ); } |
---|
327 | // extern "C" { double acosh( double ); } |
---|
328 | long double acosh( long double x ) { return acoshl( x ); } |
---|
329 | float _Complex acosh( float _Complex x ) { return cacoshf( x ); } |
---|
330 | double _Complex acosh( double _Complex x ) { return cacosh( x ); } |
---|
331 | long double _Complex acosh( long double _Complex x ) { return cacoshl( x ); } |
---|
332 | |
---|
333 | float atanh( float x ) { return atanhf( x ); } |
---|
334 | // extern "C" { double atanh( double ); } |
---|
335 | long double atanh( long double x ) { return atanhl( x ); } |
---|
336 | float _Complex atanh( float _Complex x ) { return catanhf( x ); } |
---|
337 | double _Complex atanh( double _Complex x ) { return catanh( x ); } |
---|
338 | long double _Complex atanh( long double _Complex x ) { return catanhl( x ); } |
---|
339 | } // distribution |
---|
340 | |
---|
341 | //---------------------- Error / Gamma ---------------------- |
---|
342 | |
---|
343 | static inline __attribute__((always_inline)) { |
---|
344 | float erf( float x ) { return erff( x ); } |
---|
345 | // extern "C" { double erf( double ); } |
---|
346 | long double erf( long double x ) { return erfl( x ); } |
---|
347 | // float _Complex erf( float _Complex ); |
---|
348 | // double _Complex erf( double _Complex ); |
---|
349 | // long double _Complex erf( long double _Complex ); |
---|
350 | |
---|
351 | float erfc( float x ) { return erfcf( x ); } |
---|
352 | // extern "C" { double erfc( double ); } |
---|
353 | long double erfc( long double x ) { return erfcl( x ); } |
---|
354 | // float _Complex erfc( float _Complex ); |
---|
355 | // double _Complex erfc( double _Complex ); |
---|
356 | // long double _Complex erfc( long double _Complex ); |
---|
357 | |
---|
358 | float lgamma( float x ) { return lgammaf( x ); } |
---|
359 | // extern "C" { double lgamma( double ); } |
---|
360 | long double lgamma( long double x ) { return lgammal( x ); } |
---|
361 | float lgamma( float x, int * sign ) { return lgammaf_r( x, sign ); } |
---|
362 | double lgamma( double x, int * sign ) { return lgamma_r( x, sign ); } |
---|
363 | long double lgamma( long double x, int * sign ) { return lgammal_r( x, sign ); } |
---|
364 | |
---|
365 | float tgamma( float x ) { return tgammaf( x ); } |
---|
366 | // extern "C" { double tgamma( double ); } |
---|
367 | long double tgamma( long double x ) { return tgammal( x ); } |
---|
368 | } // distribution |
---|
369 | |
---|
370 | //---------------------- Nearest Integer ---------------------- |
---|
371 | |
---|
372 | inline __attribute__((always_inline)) static { |
---|
373 | signed char floor( signed char n, signed char align ) { return n / align * align; } |
---|
374 | unsigned char floor( unsigned char n, unsigned char align ) { return n / align * align; } |
---|
375 | short int floor( short int n, short int align ) { return n / align * align; } |
---|
376 | unsigned short int floor( unsigned short int n, unsigned short int align ) { return n / align * align; } |
---|
377 | int floor( int n, int align ) { return n / align * align; } |
---|
378 | unsigned int floor( unsigned int n, unsigned int align ) { return n / align * align; } |
---|
379 | long int floor( long int n, long int align ) { return n / align * align; } |
---|
380 | unsigned long int floor( unsigned long int n, unsigned long int align ) { return n / align * align; } |
---|
381 | long long int floor( long long int n, long long int align ) { return n / align * align; } |
---|
382 | unsigned long long int floor( unsigned long long int n, unsigned long long int align ) { return n / align * align; } |
---|
383 | |
---|
384 | // forall( T | { T ?/?( T, T ); T ?*?( T, T ); } ) |
---|
385 | // T floor( T n, T align ) { return n / align * align; } |
---|
386 | |
---|
387 | signed char ceiling_div( signed char n, char align ) { return (n + (align - 1hh)) / align; } |
---|
388 | unsigned char ceiling_div( unsigned char n, unsigned char align ) { return (n + (align - 1hh)) / align; } |
---|
389 | short int ceiling_div( short int n, short int align ) { return (n + (align - 1h)) / align; } |
---|
390 | unsigned short int ceiling_div( unsigned short int n, unsigned short int align ) { return (n + (align - 1h)) / align; } |
---|
391 | int ceiling_div( int n, int align ) { return (n + (align - 1n)) / align; } |
---|
392 | unsigned int ceiling_div( unsigned int n, unsigned int align ) { return (n + (align - 1n)) / align; } |
---|
393 | long int ceiling_div( long int n, long int align ) { return (n + (align - 1l)) / align; } |
---|
394 | unsigned long int ceiling_div( unsigned long int n, unsigned long int align ) { return (n + (align - 1l)) / align; } |
---|
395 | long long int ceiling_div( long long int n, long long int align ) { return (n + (align - 1ll)) / align; } |
---|
396 | unsigned long long int ceiling_div( unsigned long long int n, unsigned long long int align ) { return (n + (align - 1ll)) / align; } |
---|
397 | |
---|
398 | signed char ceiling( signed char n, char align ) { |
---|
399 | align = align < 0 ? -align : align; // align must be signed |
---|
400 | typeof(n) trunc = floor( n, align ); |
---|
401 | return n >= 0 ? trunc + (trunc != n ? align : 0) : trunc; |
---|
402 | } |
---|
403 | unsigned char ceiling( unsigned char n, unsigned char align ) { |
---|
404 | typeof(n) trunc = floor( n, align ); |
---|
405 | return n >= 0 ? trunc + (trunc != n ? align : 0) : trunc; |
---|
406 | } |
---|
407 | short int ceiling( short int n, short int align ) { |
---|
408 | align = align < 0 ? -align : align; // align must be signed |
---|
409 | typeof(n) trunc = floor( n, align ); |
---|
410 | return n >= 0 ? trunc + (trunc != n ? align : 0) : trunc; |
---|
411 | } |
---|
412 | unsigned short int ceiling( unsigned short int n, unsigned short int align ) { |
---|
413 | typeof(n) trunc = floor( n, align ); |
---|
414 | return n >= 0 ? trunc + (trunc != n ? align : 0) : trunc; |
---|
415 | } |
---|
416 | int ceiling( int n, int align ) { |
---|
417 | align = align < 0 ? -align : align; // align must be signed |
---|
418 | typeof(n) trunc = floor( n, align ); |
---|
419 | return n >= 0 ? trunc + (trunc != n ? align : 0) : trunc; |
---|
420 | } |
---|
421 | unsigned int ceiling( unsigned int n, unsigned int align ) { |
---|
422 | typeof(n) trunc = floor( n, align ); |
---|
423 | return n >= 0 ? trunc + (trunc != n ? align : 0) : trunc; |
---|
424 | } |
---|
425 | long int ceiling( long int n, long int align ) { |
---|
426 | align = align < 0 ? -align : align; // align must be signed |
---|
427 | typeof(n) trunc = floor( n, align ); |
---|
428 | return n >= 0 ? trunc + (trunc != n ? align : 0) : trunc; |
---|
429 | } |
---|
430 | unsigned long int ceiling( unsigned long int n, unsigned long int align ) { |
---|
431 | typeof(n) trunc = floor( n, align ); |
---|
432 | return n >= 0 ? trunc + (trunc != n ? align : 0) : trunc; |
---|
433 | } |
---|
434 | long long int ceiling( long long int n, signed long long int align ) { |
---|
435 | align = align < 0 ? -align : align; // align must be signed |
---|
436 | typeof(n) trunc = floor( n, align ); |
---|
437 | return n >= 0 ? trunc + (trunc != n ? align : 0) : trunc; |
---|
438 | } |
---|
439 | unsigned long long int ceiling( unsigned long long int n, unsigned long long int align ) { |
---|
440 | typeof(n) trunc = floor( n, align ); |
---|
441 | return n >= 0 ? trunc + (trunc != n ? align : 0) : trunc; |
---|
442 | } |
---|
443 | |
---|
444 | float floor( float x ) { return floorf( x ); } |
---|
445 | // extern "C" { double floor( double ); } |
---|
446 | long double floor( long double x ) { return floorl( x ); } |
---|
447 | |
---|
448 | float ceil( float x ) { return ceilf( x ); } |
---|
449 | // extern "C" { double ceil( double ); } |
---|
450 | long double ceil( long double x ) { return ceill( x ); } |
---|
451 | |
---|
452 | float trunc( float x ) { return truncf( x ); } |
---|
453 | // extern "C" { double trunc( double ); } |
---|
454 | long double trunc( long double x ) { return truncl( x ); } |
---|
455 | |
---|
456 | float rint( float x ) { return rintf( x ); } |
---|
457 | // extern "C" { double rint( double x ); } |
---|
458 | long double rint( long double x ) { return rintl( x ); } |
---|
459 | long int rint( float x ) { return lrintf( x ); } |
---|
460 | long int rint( double x ) { return lrint( x ); } |
---|
461 | long int rint( long double x ) { return lrintl( x ); } |
---|
462 | long long int rint( float x ) { return llrintf( x ); } |
---|
463 | long long int rint( double x ) { return llrint( x ); } |
---|
464 | long long int rint( long double x ) { return llrintl( x ); } |
---|
465 | |
---|
466 | long int lrint( float x ) { return lrintf( x ); } |
---|
467 | // extern "C" { long int lrint( double ); } |
---|
468 | long int lrint( long double x ) { return lrintl( x ); } |
---|
469 | long long int llrint( float x ) { return llrintf( x ); } |
---|
470 | // extern "C" { long long int llrint( double ); } |
---|
471 | long long int llrint( long double x ) { return llrintl( x ); } |
---|
472 | |
---|
473 | float nearbyint( float x ) { return nearbyintf( x ); } |
---|
474 | // extern "C" { double nearbyint( double ); } |
---|
475 | long double nearbyint( long double x ) { return nearbyintl( x ); } |
---|
476 | |
---|
477 | float round( float x ) { return roundf( x ); } |
---|
478 | // extern "C" { double round( double x ); } |
---|
479 | long double round( long double x ) { return roundl( x ); } |
---|
480 | long int round( float x ) { return lroundf( x ); } |
---|
481 | long int round( double x ) { return lround( x ); } |
---|
482 | long int round( long double x ) { return lroundl( x ); } |
---|
483 | long long int round( float x ) { return llroundf( x ); } |
---|
484 | long long int round( double x ) { return llround( x ); } |
---|
485 | long long int round( long double x ) { return llroundl( x ); } |
---|
486 | |
---|
487 | long int lround( float x ) { return lroundf( x ); } |
---|
488 | // extern "C" { long int lround( double ); } |
---|
489 | long int lround( long double x ) { return lroundl( x ); } |
---|
490 | long long int llround( float x ) { return llroundf( x ); } |
---|
491 | // extern "C" { long long int llround( double ); } |
---|
492 | long long int llround( long double x ) { return llroundl( x ); } |
---|
493 | } // distribution |
---|
494 | |
---|
495 | //---------------------- Manipulation ---------------------- |
---|
496 | |
---|
497 | static inline __attribute__((always_inline)) { |
---|
498 | float copysign( float x, float y ) { return copysignf( x, y ); } |
---|
499 | // extern "C" { double copysign( double, double ); } |
---|
500 | long double copysign( long double x, long double y ) { return copysignl( x, y ); } |
---|
501 | |
---|
502 | float frexp( float x, int * ip ) { return frexpf( x, ip ); } |
---|
503 | // extern "C" { double frexp( double, int * ); } |
---|
504 | long double frexp( long double x, int * ip ) { return frexpl( x, ip ); } |
---|
505 | |
---|
506 | float ldexp( float x, int exp2 ) { return ldexpf( x, exp2 ); } |
---|
507 | // extern "C" { double ldexp( double, int ); } |
---|
508 | long double ldexp( long double x, int exp2 ) { return ldexpl( x, exp2 ); } |
---|
509 | |
---|
510 | [ float, float ] modf( float x ) { float i; x = modff( x, &i ); return [ i, x ]; } |
---|
511 | float modf( float x, float * i ) { return modff( x, i ); } |
---|
512 | [ double, double ] modf( double x ) { double i; x = modf( x, &i ); return [ i, x ]; } |
---|
513 | // extern "C" { double modf( double, double * ); } |
---|
514 | [ long double, long double ] modf( long double x ) { long double i; x = modfl( x, &i ); return [ i, x ]; } |
---|
515 | long double modf( long double x, long double * i ) { return modfl( x, i ); } |
---|
516 | |
---|
517 | float nextafter( float x, float y ) { return nextafterf( x, y ); } |
---|
518 | // extern "C" { double nextafter( double, double ); } |
---|
519 | long double nextafter( long double x, long double y ) { return nextafterl( x, y ); } |
---|
520 | |
---|
521 | float nexttoward( float x, long double y ) { return nexttowardf( x, y ); } |
---|
522 | // extern "C" { double nexttoward( double, long double ); } |
---|
523 | long double nexttoward( long double x, long double y ) { return nexttowardl( x, y ); } |
---|
524 | |
---|
525 | float scalbn( float x, int exp ) { return scalbnf( x, exp ); } |
---|
526 | // extern "C" { double scalbn( double, int ); } |
---|
527 | long double scalbn( long double x, int exp ) { return scalbnl( x, exp ); } |
---|
528 | float scalbn( float x, long int exp ) { return scalblnf( x, exp ); } |
---|
529 | double scalbn( double x, long int exp ) { return scalbln( x, exp ); } |
---|
530 | long double scalbn( long double x, long int exp ) { return scalblnl( x, exp ); } |
---|
531 | |
---|
532 | float scalbln( float x, long int exp ) { return scalblnf( x, exp ); } |
---|
533 | // extern "C" { double scalbln( double, long int ); } |
---|
534 | long double scalbln( long double x, long int exp ) { return scalblnl( x, exp ); } |
---|
535 | } // distribution |
---|
536 | |
---|
537 | //--------------------------------------- |
---|
538 | |
---|
539 | static inline __attribute__((always_inline)) { |
---|
540 | forall( T | { void ?{}( T &, one_t ); T ?+?( T, T ); T ?-?( T, T );T ?*?( T, T ); } ) |
---|
541 | T lerp( T x, T y, T a ) { return x * ((T){1} - a) + y * a; } |
---|
542 | |
---|
543 | forall( T | { void ?{}( T &, zero_t ); void ?{}( T &, one_t ); int ?<?( T, T ); } ) |
---|
544 | T step( T edge, T x ) { return x < edge ? (T){0} : (T){1}; } |
---|
545 | |
---|
546 | forall( T | { void ?{}( T &, int ); T clamp( T, T, T ); T ?-?( T, T ); T ?*?( T, T ); T ?/?( T, T ); } ) |
---|
547 | T smoothstep( T edge0, T edge1, T x ) { T t = clamp( (x - edge0) / (edge1 - edge0), (T){0}, (T){1} ); return t * t * ((T){3} - (T){2} * t); } |
---|
548 | } // distribution |
---|
549 | |
---|
550 | // Local Variables: // |
---|
551 | // mode: c // |
---|
552 | // tab-width: 4 // |
---|
553 | // End: // |
---|