source: libcfa/src/math.hfa @ fa5e1aa5

Last change on this file since fa5e1aa5 was 95bda0a, checked in by Peter A. Buhr <pabuhr@…>, 2 years ago

add attribute always_inline to many CFA-library cover-routines to force inlining

  • Property mode set to 100644
File size: 25.5 KB
Line 
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 Oct  8 08:40:42 2022
13// Update Count     : 136
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
28static 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
66static 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
95static 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
150static 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}
239
240//---------------------- Trigonometric ----------------------
241
242static 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
297static 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
343static 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
372inline __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 - 1)) / align; }
388        unsigned char ceiling_div( unsigned char n, unsigned char align ) { return (n + (align - 1)) / align; }
389        short int ceiling_div( short int n, short int align ) { return (n + (align - 1)) / align; }
390        unsigned short int ceiling_div( unsigned short int n, unsigned short int align ) { return (n + (align - 1)) / align; }
391        int ceiling_div( int n, int align ) { return (n + (align - 1)) / align; }
392        unsigned int ceiling_div( unsigned int n, unsigned int align ) { return (n + (align - 1)) / align; }
393        long int ceiling_div( long int n, long int align ) { return (n + (align - 1)) / align; }
394        unsigned long int ceiling_div( unsigned long int n, unsigned long int align ) { return (n + (align - 1)) / align; }
395        long long int ceiling_div( long long int n, long long int align ) { return (n + (align - 1)) / align; }
396        unsigned long long int ceiling_div( unsigned long long int n, unsigned long long int align ) { return (n + (align - 1)) / align; }
397
398        // forall( T | { T ?+?( T, T ); T ?-?( T, T ); T ?%?( T, T ); } )
399        // T ceiling_div( T n, T align ) { verify( is_pow2( align ) );return (n + (align - 1)) / align; }
400
401        // gcc notices the div/mod pair and saves both so only one div.
402        signed char ceiling( signed char n, signed char align ) { return floor( n + (n % align != 0 ? align - 1 : 0), align ); }
403        unsigned char ceiling( unsigned char n, unsigned char align ) { return floor( n + (n % align != 0 ? align - 1 : 0), align ); }
404        short int ceiling( short int n, short int align ) { return floor( n + (n % align != 0 ? align - 1 : 0), align ); }
405        unsigned short int ceiling( unsigned short int n, unsigned short int align ) { return floor( n + (n % align != 0 ? align - 1 : 0), align ); }
406        int ceiling( int n, int align ) { return floor( n + (n % align != 0 ? align - 1 : 0), align ); }
407        unsigned int ceiling( unsigned int n, unsigned int align ) { return floor( n + (n % align != 0 ? align - 1 : 0), align ); }
408        long int ceiling( long int n, long int align ) { return floor( n + (n % align != 0 ? align - 1 : 0), align ); }
409        unsigned long int ceiling( unsigned long int n, unsigned long int align ) { return floor( n + (n % align != 0 ? align - 1 : 0) , align); }
410        long long int ceiling( long long int n, long long int align ) { return floor( n + (n % align != 0 ? align - 1 : 0), align ); }
411        unsigned long long int ceiling( unsigned long long int n, unsigned long long int align ) { return floor( n + (n % align != 0 ? align - 1 : 0), align ); }
412
413        // forall( T | { void ?{}( T &, one_t ); T ?+?( T, T ); T ?-?( T, T ); T ?/?( T, T ); } )
414        // T ceiling( T n, T align ) { return return floor( n + (n % align != 0 ? align - 1 : 0), align ); *}
415
416        float floor( float x ) { return floorf( x ); }
417        // extern "C" { double floor( double ); }
418        long double floor( long double x ) { return floorl( x ); }
419
420        float ceil( float x ) { return ceilf( x ); }
421        // extern "C" { double ceil( double ); }
422        long double ceil( long double x ) { return ceill( x ); }
423
424        float trunc( float x ) { return truncf( x ); }
425        // extern "C" { double trunc( double ); }
426        long double trunc( long double x ) { return truncl( x ); }
427
428        float rint( float x ) { return rintf( x ); }
429        // extern "C" { double rint( double x ); }
430        long double rint( long double x ) { return rintl( x ); }
431        long int rint( float x ) { return lrintf( x ); }
432        long int rint( double x ) { return lrint( x ); }
433        long int rint( long double x ) { return lrintl( x ); }
434        long long int rint( float x ) { return llrintf( x ); }
435        long long int rint( double x ) { return llrint( x ); }
436        long long int rint( long double x ) { return llrintl( x ); }
437
438        long int lrint( float x ) { return lrintf( x ); }
439        // extern "C" { long int lrint( double ); }
440        long int lrint( long double x ) { return lrintl( x ); }
441        long long int llrint( float x ) { return llrintf( x ); }
442        // extern "C" { long long int llrint( double ); }
443        long long int llrint( long double x ) { return llrintl( x ); }
444
445        float nearbyint( float x ) { return nearbyintf( x ); }
446        // extern "C" { double nearbyint( double ); }
447        long double nearbyint( long double x ) { return nearbyintl( x ); }
448
449        float round( float x ) { return roundf( x ); }
450        // extern "C" { double round( double x ); }
451        long double round( long double x ) { return roundl( x ); }
452        long int round( float x ) { return lroundf( x ); }
453        long int round( double x ) { return lround( x ); }
454        long int round( long double x ) { return lroundl( x ); }
455        long long int round( float x ) { return llroundf( x ); }
456        long long int round( double x ) { return llround( x ); }
457        long long int round( long double x ) { return llroundl( x ); }
458
459        long int lround( float x ) { return lroundf( x ); }
460        // extern "C" { long int lround( double ); }
461        long int lround( long double x ) { return lroundl( x ); }
462        long long int llround( float x ) { return llroundf( x ); }
463        // extern "C" { long long int llround( double ); }
464        long long int llround( long double x ) { return llroundl( x ); }
465} // distribution
466
467//---------------------- Manipulation ----------------------
468
469static inline __attribute__((always_inline)) {
470        float copysign( float x, float y ) { return copysignf( x, y ); }
471        // extern "C" { double copysign( double, double ); }
472        long double copysign( long double x, long double y ) { return copysignl( x, y ); }
473
474        float frexp( float x, int * ip ) { return frexpf( x, ip ); }
475        // extern "C" { double frexp( double, int * ); }
476        long double frexp( long double x, int * ip ) { return frexpl( x, ip ); }
477
478        float ldexp( float x, int exp2 ) { return ldexpf( x, exp2 ); }
479        // extern "C" { double ldexp( double, int ); }
480        long double ldexp( long double x, int exp2 ) { return ldexpl( x, exp2 ); }
481
482        [ float, float ] modf( float x ) { float i; x = modff( x, &i ); return [ i, x ]; }
483        float modf( float x, float * i ) { return modff( x, i ); }
484        [ double, double ] modf( double x ) { double i; x = modf( x, &i ); return [ i, x ]; }
485        // extern "C" { double modf( double, double * ); }
486        [ long double, long double ] modf( long double x ) { long double i; x = modfl( x, &i ); return [ i, x ]; }
487        long double modf( long double x, long double * i ) { return modfl( x, i ); }
488
489        float nextafter( float x, float y ) { return nextafterf( x, y ); }
490        // extern "C" { double nextafter( double, double ); }
491        long double nextafter( long double x, long double y ) { return nextafterl( x, y ); }
492
493        float nexttoward( float x, long double y ) { return nexttowardf( x, y ); }
494        // extern "C" { double nexttoward( double, long double ); }
495        long double nexttoward( long double x, long double y ) { return nexttowardl( x, y ); }
496
497        float scalbn( float x, int exp ) { return scalbnf( x, exp ); }
498        // extern "C" { double scalbn( double, int ); }
499        long double scalbn( long double x, int exp ) { return scalbnl( x, exp ); }
500        float scalbn( float x, long int exp ) { return scalblnf( x, exp ); }
501        double scalbn( double x, long int exp ) { return scalbln( x, exp ); }
502        long double scalbn( long double x, long int exp ) { return scalblnl( x, exp ); }
503
504        float scalbln( float x, long int exp ) { return scalblnf( x, exp ); }
505        // extern "C" { double scalbln( double, long int ); }
506        long double scalbln( long double x, long int exp ) { return scalblnl( x, exp ); }
507} // distribution
508
509//---------------------------------------
510
511static inline __attribute__((always_inline)) {
512        forall( T | { void ?{}( T &, one_t ); T ?+?( T, T ); T ?-?( T, T );T ?*?( T, T ); } )
513        T lerp( T x, T y, T a ) { return x * ((T){1} - a) + y * a; }
514
515        forall( T | { void ?{}( T &, zero_t ); void ?{}( T &, one_t ); int ?<?( T, T ); } )
516        T step( T edge, T x ) { return x < edge ? (T){0} : (T){1}; }
517
518        forall( T | { void ?{}( T &, int ); T clamp( T, T, T ); T ?-?( T, T ); T ?*?( T, T ); T ?/?( T, T ); } )
519        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); }
520} // distribution
521
522// Local Variables: //
523// mode: c //
524// tab-width: 4 //
525// End: //
Note: See TracBrowser for help on using the repository browser.