source: libcfa/src/math.hfa

Last change on this file was 658f3179, checked in by Andrew Beach <ajbeach@…>, 12 days ago

Moved massive function log2_u32_32 out of header.
  • Property mode set to 100644
File size: 21.7 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
150unsigned long long log2_u32_32( unsigned long long val );
151
152//---------------------- Trigonometric ----------------------
153
154static inline __attribute__((always_inline)) {
155        float sin( float x ) { return sinf( x ); }
156        // extern "C" { double sin( double ); }
157        long double sin( long double x ) { return sinl( x ); }
158        float _Complex sin( float _Complex x ) { return csinf( x ); }
159        double _Complex sin( double _Complex x ) { return csin( x ); }
160        long double _Complex sin( long double _Complex x ) { return csinl( x ); }
161
162        float cos( float x ) { return cosf( x ); }
163        // extern "C" { double cos( double ); }
164        long double cos( long double x ) { return cosl( x ); }
165        float _Complex cos( float _Complex x ) { return ccosf( x ); }
166        double _Complex cos( double _Complex x ) { return ccos( x ); }
167        long double _Complex cos( long double _Complex x ) { return ccosl( x ); }
168
169        float tan( float x ) { return tanf( x ); }
170        // extern "C" { double tan( double ); }
171        long double tan( long double x ) { return tanl( x ); }
172        float _Complex tan( float _Complex x ) { return ctanf( x ); }
173        double _Complex tan( double _Complex x ) { return ctan( x ); }
174        long double _Complex tan( long double _Complex x ) { return ctanl( x ); }
175
176        float asin( float x ) { return asinf( x ); }
177        // extern "C" { double asin( double ); }
178        long double asin( long double x ) { return asinl( x ); }
179        float _Complex asin( float _Complex x ) { return casinf( x ); }
180        double _Complex asin( double _Complex x ) { return casin( x ); }
181        long double _Complex asin( long double _Complex x ) { return casinl( x ); }
182
183        float acos( float x ) { return acosf( x ); }
184        // extern "C" { double acos( double ); }
185        long double acos( long double x ) { return acosl( x ); }
186        float _Complex acos( float _Complex x ) { return cacosf( x ); }
187        double _Complex acos( double _Complex x ) { return cacos( x ); }
188        long double _Complex acos( long double _Complex x ) { return cacosl( x ); }
189
190        float atan( float x ) { return atanf( x ); }
191        // extern "C" { double atan( double ); }
192        long double atan( long double x ) { return atanl( x ); }
193        float _Complex atan( float _Complex x ) { return catanf( x ); }
194        double _Complex atan( double _Complex x ) { return catan( x ); }
195        long double _Complex atan( long double _Complex x ) { return catanl( x ); }
196
197        float atan2( float x, float y ) { return atan2f( x, y ); }
198        // extern "C" { double atan2( double, double ); }
199        long double atan2( long double x, long double y ) { return atan2l( x, y ); }
200
201        // alternative name for atan2
202        float atan( float x, float y ) { return atan2f( x, y ); }
203        double atan( double x, double y ) { return atan2( x, y ); }
204        long double atan( long double x, long double y ) { return atan2l( x, y ); }
205} // distribution
206
207//---------------------- Hyperbolic ----------------------
208
209static inline __attribute__((always_inline)) {
210        float sinh( float x ) { return sinhf( x ); }
211        // extern "C" { double sinh( double ); }
212        long double sinh( long double x ) { return sinhl( x ); }
213        float _Complex sinh( float _Complex x ) { return csinhf( x ); }
214        double _Complex sinh( double _Complex x ) { return csinh( x ); }
215        long double _Complex sinh( long double _Complex x ) { return csinhl( x ); }
216
217        float cosh( float x ) { return coshf( x ); }
218        // extern "C" { double cosh( double ); }
219        long double cosh( long double x ) { return coshl( x ); }
220        float _Complex cosh( float _Complex x ) { return ccoshf( x ); }
221        double _Complex cosh( double _Complex x ) { return ccosh( x ); }
222        long double _Complex cosh( long double _Complex x ) { return ccoshl( x ); }
223
224        float tanh( float x ) { return tanhf( x ); }
225        // extern "C" { double tanh( double ); }
226        long double tanh( long double x ) { return tanhl( x ); }
227        float _Complex tanh( float _Complex x ) { return ctanhf( x ); }
228        double _Complex tanh( double _Complex x ) { return ctanh( x ); }
229        long double _Complex tanh( long double _Complex x ) { return ctanhl( x ); }
230
231        float asinh( float x ) { return asinhf( x ); }
232        // extern "C" { double asinh( double ); }
233        long double asinh( long double x ) { return asinhl( x ); }
234        float _Complex asinh( float _Complex x ) { return casinhf( x ); }
235        double _Complex asinh( double _Complex x ) { return casinh( x ); }
236        long double _Complex asinh( long double _Complex x ) { return casinhl( x ); }
237
238        float acosh( float x ) { return acoshf( x ); }
239        // extern "C" { double acosh( double ); }
240        long double acosh( long double x ) { return acoshl( x ); }
241        float _Complex acosh( float _Complex x ) { return cacoshf( x ); }
242        double _Complex acosh( double _Complex x ) { return cacosh( x ); }
243        long double _Complex acosh( long double _Complex x ) { return cacoshl( x ); }
244
245        float atanh( float x ) { return atanhf( x ); }
246        // extern "C" { double atanh( double ); }
247        long double atanh( long double x ) { return atanhl( x ); }
248        float _Complex atanh( float _Complex x ) { return catanhf( x ); }
249        double _Complex atanh( double _Complex x ) { return catanh( x ); }
250        long double _Complex atanh( long double _Complex x ) { return catanhl( x ); }
251} // distribution
252
253//---------------------- Error / Gamma ----------------------
254
255static inline __attribute__((always_inline)) {
256        float erf( float x ) { return erff( x ); }
257        // extern "C" { double erf( double ); }
258        long double erf( long double x ) { return erfl( x ); }
259        // float _Complex erf( float _Complex );
260        // double _Complex erf( double _Complex );
261        // long double _Complex erf( long double _Complex );
262
263        float erfc( float x ) { return erfcf( x ); }
264        // extern "C" { double erfc( double ); }
265        long double erfc( long double x ) { return erfcl( x ); }
266        // float _Complex erfc( float _Complex );
267        // double _Complex erfc( double _Complex );
268        // long double _Complex erfc( long double _Complex );
269
270        float lgamma( float x ) { return lgammaf( x ); }
271        // extern "C" { double lgamma( double ); }
272        long double lgamma( long double x ) { return lgammal( x ); }
273        float lgamma( float x, int * sign ) { return lgammaf_r( x, sign ); }
274        double lgamma( double x, int * sign ) { return lgamma_r( x, sign ); }
275        long double lgamma( long double x, int * sign ) { return lgammal_r( x, sign ); }
276
277        float tgamma( float x ) { return tgammaf( x ); }
278        // extern "C" { double tgamma( double ); }
279        long double tgamma( long double x ) { return tgammal( x ); }
280} // distribution
281
282//---------------------- Nearest Integer ----------------------
283
284inline __attribute__((always_inline)) static {
285        // force divide before multiply
286        signed char floor( signed char n, signed char align ) { return (n / align) * align; }
287        unsigned char floor( unsigned char n, unsigned char align ) { return (n / align) * align; }
288        short int floor( short int n, short int align ) { return (n / align) * align; }
289        unsigned short int floor( unsigned short int n, unsigned short int align ) { return (n / align) * align; }
290        int floor( int n, int align ) { return (n / align) * align; }
291        unsigned int floor( unsigned int n, unsigned int align ) { return (n / align) * align; }
292        long int floor( long int n, long int align ) { return (n / align) * align; }
293        unsigned long int floor( unsigned long int n, unsigned long int align ) { return (n / align) * align; }
294        long long int floor( long long int n, long long int align ) { return (n / align) * align; }
295        unsigned long long int floor( unsigned long long int n, unsigned long long int align ) { return (n / align) * align; }
296
297        // forall( T | { T ?/?( T, T ); T ?*?( T, T ); } )
298        // T floor( T n, T align ) { return (n / align) * align; }
299
300        signed char ceiling_div( signed char n, char align ) { return (n + (align - 1hh)) / align; }
301        unsigned char ceiling_div( unsigned char n, unsigned char align ) { return (n + (align - 1hhu)) / align; }
302        short int ceiling_div( short int n, short int align ) { return (n + (align - 1h)) / align; }
303        unsigned short int ceiling_div( unsigned short int n, unsigned short int align ) { return (n + (align - 1hu)) / align; }
304        int ceiling_div( int n, int align ) { return (n + (align - 1n)) / align; }
305        unsigned int ceiling_div( unsigned int n, unsigned int align ) { return (n + (align - 1nu)) / align; }
306        long int ceiling_div( long int n, long int align ) { return (n + (align - 1l)) / align; }
307        unsigned long int ceiling_div( unsigned long int n, unsigned long int align ) { return (n + (align - 1lu)) / align; }
308        long long int ceiling_div( long long int n, long long int align ) { return (n + (align - 1ll)) / align; }
309        unsigned long long int ceiling_div( unsigned long long int n, unsigned long long int align ) { return (n + (align - 1llu)) / align; }
310
311        signed char ceiling( signed char n, char align ) {
312                typeof(n) trunc = floor( n, align );
313                return n < 0 || n == trunc ? trunc : trunc + align;
314        }
315        unsigned char ceiling( unsigned char n, unsigned char align ) {
316                typeof(n) trunc = floor( n, align );
317                return n == trunc ? trunc : trunc + align;
318        }
319        short int ceiling( short int n, short int align ) {
320                typeof(n) trunc = floor( n, align );
321                return n < 0 || n == trunc ? trunc : trunc + align;
322        }
323        unsigned short int ceiling( unsigned short int n, unsigned short int align ) {
324                typeof(n) trunc = floor( n, align );
325                return n == trunc ? trunc : trunc + align;
326        }
327        int ceiling( int n, int align ) {
328                typeof(n) trunc = floor( n, align );
329                return n < 0 || n == trunc ? trunc : trunc + align;
330        }
331        unsigned int ceiling( unsigned int n, unsigned int align ) {
332                typeof(n) trunc = floor( n, align );
333                return n == trunc ? trunc : trunc + align;
334        }
335        long int ceiling( long int n, long int align ) {
336                typeof(n) trunc = floor( n, align );
337                return n < 0 || n == trunc ? trunc : trunc + align;
338        }
339        unsigned long int ceiling( unsigned long int n, unsigned long int align ) {
340                typeof(n) trunc = floor( n, align );
341                return n == trunc ? trunc : trunc + align;
342        }
343        long long int ceiling( long long int n, signed long long int align ) {
344                typeof(n) trunc = floor( n, align );
345                return n < 0 || n == trunc ? trunc : trunc + align;
346        }
347        unsigned long long int ceiling( unsigned long long int n, unsigned long long int align ) {
348                typeof(n) trunc = floor( n, align );
349                return n == trunc ? trunc : trunc + align;
350        }
351
352        float floor( float x ) { return floorf( x ); }
353        // extern "C" { double floor( double ); }
354        long double floor( long double x ) { return floorl( x ); }
355
356        float ceil( float x ) { return ceilf( x ); }
357        // extern "C" { double ceil( double ); }
358        long double ceil( long double x ) { return ceill( x ); }
359
360        float trunc( float x ) { return truncf( x ); }
361        // extern "C" { double trunc( double ); }
362        long double trunc( long double x ) { return truncl( x ); }
363
364        float rint( float x ) { return rintf( x ); }
365        // extern "C" { double rint( double x ); }
366        long double rint( long double x ) { return rintl( x ); }
367        long int rint( float x ) { return lrintf( x ); }
368        long int rint( double x ) { return lrint( x ); }
369        long int rint( long double x ) { return lrintl( x ); }
370        long long int rint( float x ) { return llrintf( x ); }
371        long long int rint( double x ) { return llrint( x ); }
372        long long int rint( long double x ) { return llrintl( x ); }
373
374        long int lrint( float x ) { return lrintf( x ); }
375        // extern "C" { long int lrint( double ); }
376        long int lrint( long double x ) { return lrintl( x ); }
377        long long int llrint( float x ) { return llrintf( x ); }
378        // extern "C" { long long int llrint( double ); }
379        long long int llrint( long double x ) { return llrintl( x ); }
380
381        float nearbyint( float x ) { return nearbyintf( x ); }
382        // extern "C" { double nearbyint( double ); }
383        long double nearbyint( long double x ) { return nearbyintl( x ); }
384
385        float round( float x ) { return roundf( x ); }
386        // extern "C" { double round( double x ); }
387        long double round( long double x ) { return roundl( x ); }
388        long int round( float x ) { return lroundf( x ); }
389        long int round( double x ) { return lround( x ); }
390        long int round( long double x ) { return lroundl( x ); }
391        long long int round( float x ) { return llroundf( x ); }
392        long long int round( double x ) { return llround( x ); }
393        long long int round( long double x ) { return llroundl( x ); }
394
395        long int lround( float x ) { return lroundf( x ); }
396        // extern "C" { long int lround( double ); }
397        long int lround( long double x ) { return lroundl( x ); }
398        long long int llround( float x ) { return llroundf( x ); }
399        // extern "C" { long long int llround( double ); }
400        long long int llround( long double x ) { return llroundl( x ); }
401} // distribution
402
403//---------------------- Manipulation ----------------------
404
405static inline __attribute__((always_inline)) {
406        float copysign( float x, float y ) { return copysignf( x, y ); }
407        // extern "C" { double copysign( double, double ); }
408        long double copysign( long double x, long double y ) { return copysignl( x, y ); }
409
410        float frexp( float x, int * ip ) { return frexpf( x, ip ); }
411        // extern "C" { double frexp( double, int * ); }
412        long double frexp( long double x, int * ip ) { return frexpl( x, ip ); }
413
414        float ldexp( float x, int exp2 ) { return ldexpf( x, exp2 ); }
415        // extern "C" { double ldexp( double, int ); }
416        long double ldexp( long double x, int exp2 ) { return ldexpl( x, exp2 ); }
417
418        [ float, float ] modf( float x ) { float i; x = modff( x, &i ); return [ i, x ]; }
419        float modf( float x, float * i ) { return modff( x, i ); }
420        [ double, double ] modf( double x ) { double i; x = modf( x, &i ); return [ i, x ]; }
421        // extern "C" { double modf( double, double * ); }
422        [ long double, long double ] modf( long double x ) { long double i; x = modfl( x, &i ); return [ i, x ]; }
423        long double modf( long double x, long double * i ) { return modfl( x, i ); }
424
425        float nextafter( float x, float y ) { return nextafterf( x, y ); }
426        // extern "C" { double nextafter( double, double ); }
427        long double nextafter( long double x, long double y ) { return nextafterl( x, y ); }
428
429        float nexttoward( float x, long double y ) { return nexttowardf( x, y ); }
430        // extern "C" { double nexttoward( double, long double ); }
431        long double nexttoward( long double x, long double y ) { return nexttowardl( x, y ); }
432
433        float scalbn( float x, int exp ) { return scalbnf( x, exp ); }
434        // extern "C" { double scalbn( double, int ); }
435        long double scalbn( long double x, int exp ) { return scalbnl( x, exp ); }
436        float scalbn( float x, long int exp ) { return scalblnf( x, exp ); }
437        double scalbn( double x, long int exp ) { return scalbln( x, exp ); }
438        long double scalbn( long double x, long int exp ) { return scalblnl( x, exp ); }
439
440        float scalbln( float x, long int exp ) { return scalblnf( x, exp ); }
441        // extern "C" { double scalbln( double, long int ); }
442        long double scalbln( long double x, long int exp ) { return scalblnl( x, exp ); }
443} // distribution
444
445//---------------------------------------
446
447static inline __attribute__((always_inline)) {
448        forall( T | { void ?{}( T &, one_t ); T ?+?( T, T ); T ?-?( T, T );T ?*?( T, T ); } )
449        T lerp( T x, T y, T a ) { return x * ((T){1} - a) + y * a; }
450
451        forall( T | { void ?{}( T &, zero_t ); void ?{}( T &, one_t ); int ?<?( T, T ); } )
452        T step( T edge, T x ) { return x < edge ? (T){0} : (T){1}; }
453
454        forall( T | { void ?{}( T &, int ); T clamp( T, T, T ); T ?-?( T, T ); T ?*?( T, T ); T ?/?( T, T ); } )
455        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); }
456} // distribution
457
458// Local Variables: //
459// mode: c //
460// tab-width: 4 //
461// End: //
Note: See TracBrowser for help on using the repository browser.