source: libcfa/src/math.hfa

Last change on this file was 7c012e8, checked in by Peter A. Buhr <pabuhr@…>, 8 months ago

simplify computation for ceiling, add unsigned qualifiers on one_t constants in ceiling_div

  • 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 : Sun Jun 18 08:13:53 2023
13// Update Count     : 202
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} // log2_u32_32
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        // force divide before multiply
374        signed char floor( signed char n, signed char align ) { return (n / align) * align; }
375        unsigned char floor( unsigned char n, unsigned char align ) { return (n / align) * align; }
376        short int floor( short int n, short int align ) { return (n / align) * align; }
377        unsigned short int floor( unsigned short int n, unsigned short int align ) { return (n / align) * align; }
378        int floor( int n, int align ) { return (n / align) * align; }
379        unsigned int floor( unsigned int n, unsigned int align ) { return (n / align) * align; }
380        long int floor( long int n, long int align ) { return (n / align) * align; }
381        unsigned long int floor( unsigned long int n, unsigned long int align ) { return (n / align) * align; }
382        long long int floor( long long int n, long long int align ) { return (n / align) * align; }
383        unsigned long long int floor( unsigned long long int n, unsigned long long int align ) { return (n / align) * align; }
384
385        // forall( T | { T ?/?( T, T ); T ?*?( T, T ); } )
386        // T floor( T n, T align ) { return (n / align) * align; }
387
388        signed char ceiling_div( signed char n, char align ) { return (n + (align - 1hh)) / align; }
389        unsigned char ceiling_div( unsigned char n, unsigned char align ) { return (n + (align - 1hhu)) / align; }
390        short int ceiling_div( short int n, short int align ) { return (n + (align - 1h)) / align; }
391        unsigned short int ceiling_div( unsigned short int n, unsigned short int align ) { return (n + (align - 1hu)) / align; }
392        int ceiling_div( int n, int align ) { return (n + (align - 1n)) / align; }
393        unsigned int ceiling_div( unsigned int n, unsigned int align ) { return (n + (align - 1nu)) / align; }
394        long int ceiling_div( long int n, long int align ) { return (n + (align - 1l)) / align; }
395        unsigned long int ceiling_div( unsigned long int n, unsigned long int align ) { return (n + (align - 1lu)) / align; }
396        long long int ceiling_div( long long int n, long long int align ) { return (n + (align - 1ll)) / align; }
397        unsigned long long int ceiling_div( unsigned long long int n, unsigned long long int align ) { return (n + (align - 1llu)) / align; }
398
399        signed char ceiling( signed char n, char align ) {
400                typeof(n) trunc = floor( n, align );
401                return n < 0 || n == trunc ? trunc : trunc + align;
402        }
403        unsigned char ceiling( unsigned char n, unsigned char align ) {
404                typeof(n) trunc = floor( n, align );
405                return n == trunc ? trunc : trunc + align;
406        }
407        short int ceiling( short int n, short int align ) {
408                typeof(n) trunc = floor( n, align );
409                return n < 0 || n == trunc ? trunc : trunc + align;
410        }
411        unsigned short int ceiling( unsigned short int n, unsigned short int align ) {
412                typeof(n) trunc = floor( n, align );
413                return n == trunc ? trunc : trunc + align;
414        }
415        int ceiling( int n, int align ) {
416                typeof(n) trunc = floor( n, align );
417                return n < 0 || n == trunc ? trunc : trunc + align;
418        }
419        unsigned int ceiling( unsigned int n, unsigned int align ) {
420                typeof(n) trunc = floor( n, align );
421                return n == trunc ? trunc : trunc + align;
422        }
423        long int ceiling( long int n, long int align ) {
424                typeof(n) trunc = floor( n, align );
425                return n < 0 || n == trunc ? trunc : trunc + align;
426        }
427        unsigned long int ceiling( unsigned long int n, unsigned long int align ) {
428                typeof(n) trunc = floor( n, align );
429                return n == trunc ? trunc : trunc + align;
430        }
431        long long int ceiling( long long int n, signed long long int align ) {
432                typeof(n) trunc = floor( n, align );
433                return n < 0 || n == trunc ? trunc : trunc + align;
434        }
435        unsigned long long int ceiling( unsigned long long int n, unsigned long long int align ) {
436                typeof(n) trunc = floor( n, align );
437                return n == trunc ? trunc : trunc + align;
438        }
439
440        float floor( float x ) { return floorf( x ); }
441        // extern "C" { double floor( double ); }
442        long double floor( long double x ) { return floorl( x ); }
443
444        float ceil( float x ) { return ceilf( x ); }
445        // extern "C" { double ceil( double ); }
446        long double ceil( long double x ) { return ceill( x ); }
447
448        float trunc( float x ) { return truncf( x ); }
449        // extern "C" { double trunc( double ); }
450        long double trunc( long double x ) { return truncl( x ); }
451
452        float rint( float x ) { return rintf( x ); }
453        // extern "C" { double rint( double x ); }
454        long double rint( long double x ) { return rintl( x ); }
455        long int rint( float x ) { return lrintf( x ); }
456        long int rint( double x ) { return lrint( x ); }
457        long int rint( long double x ) { return lrintl( x ); }
458        long long int rint( float x ) { return llrintf( x ); }
459        long long int rint( double x ) { return llrint( x ); }
460        long long int rint( long double x ) { return llrintl( x ); }
461
462        long int lrint( float x ) { return lrintf( x ); }
463        // extern "C" { long int lrint( double ); }
464        long int lrint( long double x ) { return lrintl( x ); }
465        long long int llrint( float x ) { return llrintf( x ); }
466        // extern "C" { long long int llrint( double ); }
467        long long int llrint( long double x ) { return llrintl( x ); }
468
469        float nearbyint( float x ) { return nearbyintf( x ); }
470        // extern "C" { double nearbyint( double ); }
471        long double nearbyint( long double x ) { return nearbyintl( x ); }
472
473        float round( float x ) { return roundf( x ); }
474        // extern "C" { double round( double x ); }
475        long double round( long double x ) { return roundl( x ); }
476        long int round( float x ) { return lroundf( x ); }
477        long int round( double x ) { return lround( x ); }
478        long int round( long double x ) { return lroundl( x ); }
479        long long int round( float x ) { return llroundf( x ); }
480        long long int round( double x ) { return llround( x ); }
481        long long int round( long double x ) { return llroundl( x ); }
482
483        long int lround( float x ) { return lroundf( x ); }
484        // extern "C" { long int lround( double ); }
485        long int lround( long double x ) { return lroundl( x ); }
486        long long int llround( float x ) { return llroundf( x ); }
487        // extern "C" { long long int llround( double ); }
488        long long int llround( long double x ) { return llroundl( x ); }
489} // distribution
490
491//---------------------- Manipulation ----------------------
492
493static inline __attribute__((always_inline)) {
494        float copysign( float x, float y ) { return copysignf( x, y ); }
495        // extern "C" { double copysign( double, double ); }
496        long double copysign( long double x, long double y ) { return copysignl( x, y ); }
497
498        float frexp( float x, int * ip ) { return frexpf( x, ip ); }
499        // extern "C" { double frexp( double, int * ); }
500        long double frexp( long double x, int * ip ) { return frexpl( x, ip ); }
501
502        float ldexp( float x, int exp2 ) { return ldexpf( x, exp2 ); }
503        // extern "C" { double ldexp( double, int ); }
504        long double ldexp( long double x, int exp2 ) { return ldexpl( x, exp2 ); }
505
506        [ float, float ] modf( float x ) { float i; x = modff( x, &i ); return [ i, x ]; }
507        float modf( float x, float * i ) { return modff( x, i ); }
508        [ double, double ] modf( double x ) { double i; x = modf( x, &i ); return [ i, x ]; }
509        // extern "C" { double modf( double, double * ); }
510        [ long double, long double ] modf( long double x ) { long double i; x = modfl( x, &i ); return [ i, x ]; }
511        long double modf( long double x, long double * i ) { return modfl( x, i ); }
512
513        float nextafter( float x, float y ) { return nextafterf( x, y ); }
514        // extern "C" { double nextafter( double, double ); }
515        long double nextafter( long double x, long double y ) { return nextafterl( x, y ); }
516
517        float nexttoward( float x, long double y ) { return nexttowardf( x, y ); }
518        // extern "C" { double nexttoward( double, long double ); }
519        long double nexttoward( long double x, long double y ) { return nexttowardl( x, y ); }
520
521        float scalbn( float x, int exp ) { return scalbnf( x, exp ); }
522        // extern "C" { double scalbn( double, int ); }
523        long double scalbn( long double x, int exp ) { return scalbnl( x, exp ); }
524        float scalbn( float x, long int exp ) { return scalblnf( x, exp ); }
525        double scalbn( double x, long int exp ) { return scalbln( x, exp ); }
526        long double scalbn( long double x, long int exp ) { return scalblnl( x, exp ); }
527
528        float scalbln( float x, long int exp ) { return scalblnf( x, exp ); }
529        // extern "C" { double scalbln( double, long int ); }
530        long double scalbln( long double x, long int exp ) { return scalblnl( x, exp ); }
531} // distribution
532
533//---------------------------------------
534
535static inline __attribute__((always_inline)) {
536        forall( T | { void ?{}( T &, one_t ); T ?+?( T, T ); T ?-?( T, T );T ?*?( T, T ); } )
537        T lerp( T x, T y, T a ) { return x * ((T){1} - a) + y * a; }
538
539        forall( T | { void ?{}( T &, zero_t ); void ?{}( T &, one_t ); int ?<?( T, T ); } )
540        T step( T edge, T x ) { return x < edge ? (T){0} : (T){1}; }
541
542        forall( T | { void ?{}( T &, int ); T clamp( T, T, T ); T ?-?( T, T ); T ?*?( T, T ); T ?/?( T, T ); } )
543        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); }
544} // distribution
545
546// Local Variables: //
547// mode: c //
548// tab-width: 4 //
549// End: //
Note: See TracBrowser for help on using the repository browser.