source: libcfa/src/math.hfa@ d8c96a9

ADT ast-experimental
Last change on this file since d8c96a9 was 0deeaad, checked in by Thierry Delisle <tdelisle@…>, 3 years ago

Added fixed point log2 calculation, which is not that useful but kind of cool.
  • Property mode set to 100644
File size: 25.2 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 : Thu Apr 15 11:47:56 2021
13// Update Count : 132
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 {
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 {
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 {
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//---------------------- Trigonometric ----------------------
151
152static inline {
153 float sin( float x ) { return sinf( x ); }
154 // extern "C" { double sin( double ); }
155 long double sin( long double x ) { return sinl( x ); }
156 float _Complex sin( float _Complex x ) { return csinf( x ); }
157 double _Complex sin( double _Complex x ) { return csin( x ); }
158 long double _Complex sin( long double _Complex x ) { return csinl( x ); }
159
160 float cos( float x ) { return cosf( x ); }
161 // extern "C" { double cos( double ); }
162 long double cos( long double x ) { return cosl( x ); }
163 float _Complex cos( float _Complex x ) { return ccosf( x ); }
164 double _Complex cos( double _Complex x ) { return ccos( x ); }
165 long double _Complex cos( long double _Complex x ) { return ccosl( x ); }
166
167 float tan( float x ) { return tanf( x ); }
168 // extern "C" { double tan( double ); }
169 long double tan( long double x ) { return tanl( x ); }
170 float _Complex tan( float _Complex x ) { return ctanf( x ); }
171 double _Complex tan( double _Complex x ) { return ctan( x ); }
172 long double _Complex tan( long double _Complex x ) { return ctanl( x ); }
173
174 float asin( float x ) { return asinf( x ); }
175 // extern "C" { double asin( double ); }
176 long double asin( long double x ) { return asinl( x ); }
177 float _Complex asin( float _Complex x ) { return casinf( x ); }
178 double _Complex asin( double _Complex x ) { return casin( x ); }
179 long double _Complex asin( long double _Complex x ) { return casinl( x ); }
180
181 float acos( float x ) { return acosf( x ); }
182 // extern "C" { double acos( double ); }
183 long double acos( long double x ) { return acosl( x ); }
184 float _Complex acos( float _Complex x ) { return cacosf( x ); }
185 double _Complex acos( double _Complex x ) { return cacos( x ); }
186 long double _Complex acos( long double _Complex x ) { return cacosl( x ); }
187
188 float atan( float x ) { return atanf( x ); }
189 // extern "C" { double atan( double ); }
190 long double atan( long double x ) { return atanl( x ); }
191 float _Complex atan( float _Complex x ) { return catanf( x ); }
192 double _Complex atan( double _Complex x ) { return catan( x ); }
193 long double _Complex atan( long double _Complex x ) { return catanl( x ); }
194
195 float atan2( float x, float y ) { return atan2f( x, y ); }
196 // extern "C" { double atan2( double, double ); }
197 long double atan2( long double x, long double y ) { return atan2l( x, y ); }
198
199 // alternative name for atan2
200 float atan( float x, float y ) { return atan2f( x, y ); }
201 double atan( double x, double y ) { return atan2( x, y ); }
202 long double atan( long double x, long double y ) { return atan2l( x, y ); }
203} // distribution
204
205//---------------------- Hyperbolic ----------------------
206
207static inline {
208 float sinh( float x ) { return sinhf( x ); }
209 // extern "C" { double sinh( double ); }
210 long double sinh( long double x ) { return sinhl( x ); }
211 float _Complex sinh( float _Complex x ) { return csinhf( x ); }
212 double _Complex sinh( double _Complex x ) { return csinh( x ); }
213 long double _Complex sinh( long double _Complex x ) { return csinhl( x ); }
214
215 float cosh( float x ) { return coshf( x ); }
216 // extern "C" { double cosh( double ); }
217 long double cosh( long double x ) { return coshl( x ); }
218 float _Complex cosh( float _Complex x ) { return ccoshf( x ); }
219 double _Complex cosh( double _Complex x ) { return ccosh( x ); }
220 long double _Complex cosh( long double _Complex x ) { return ccoshl( x ); }
221
222 float tanh( float x ) { return tanhf( x ); }
223 // extern "C" { double tanh( double ); }
224 long double tanh( long double x ) { return tanhl( x ); }
225 float _Complex tanh( float _Complex x ) { return ctanhf( x ); }
226 double _Complex tanh( double _Complex x ) { return ctanh( x ); }
227 long double _Complex tanh( long double _Complex x ) { return ctanhl( x ); }
228
229 float asinh( float x ) { return asinhf( x ); }
230 // extern "C" { double asinh( double ); }
231 long double asinh( long double x ) { return asinhl( x ); }
232 float _Complex asinh( float _Complex x ) { return casinhf( x ); }
233 double _Complex asinh( double _Complex x ) { return casinh( x ); }
234 long double _Complex asinh( long double _Complex x ) { return casinhl( x ); }
235
236 float acosh( float x ) { return acoshf( x ); }
237 // extern "C" { double acosh( double ); }
238 long double acosh( long double x ) { return acoshl( x ); }
239 float _Complex acosh( float _Complex x ) { return cacoshf( x ); }
240 double _Complex acosh( double _Complex x ) { return cacosh( x ); }
241 long double _Complex acosh( long double _Complex x ) { return cacoshl( x ); }
242
243 float atanh( float x ) { return atanhf( x ); }
244 // extern "C" { double atanh( double ); }
245 long double atanh( long double x ) { return atanhl( x ); }
246 float _Complex atanh( float _Complex x ) { return catanhf( x ); }
247 double _Complex atanh( double _Complex x ) { return catanh( x ); }
248 long double _Complex atanh( long double _Complex x ) { return catanhl( x ); }
249} // distribution
250
251//---------------------- Error / Gamma ----------------------
252
253static inline {
254 float erf( float x ) { return erff( x ); }
255 // extern "C" { double erf( double ); }
256 long double erf( long double x ) { return erfl( x ); }
257 // float _Complex erf( float _Complex );
258 // double _Complex erf( double _Complex );
259 // long double _Complex erf( long double _Complex );
260
261 float erfc( float x ) { return erfcf( x ); }
262 // extern "C" { double erfc( double ); }
263 long double erfc( long double x ) { return erfcl( x ); }
264 // float _Complex erfc( float _Complex );
265 // double _Complex erfc( double _Complex );
266 // long double _Complex erfc( long double _Complex );
267
268 float lgamma( float x ) { return lgammaf( x ); }
269 // extern "C" { double lgamma( double ); }
270 long double lgamma( long double x ) { return lgammal( x ); }
271 float lgamma( float x, int * sign ) { return lgammaf_r( x, sign ); }
272 double lgamma( double x, int * sign ) { return lgamma_r( x, sign ); }
273 long double lgamma( long double x, int * sign ) { return lgammal_r( x, sign ); }
274
275 float tgamma( float x ) { return tgammaf( x ); }
276 // extern "C" { double tgamma( double ); }
277 long double tgamma( long double x ) { return tgammal( x ); }
278} // distribution
279
280//---------------------- Nearest Integer ----------------------
281
282static inline {
283 signed char floor( signed char n, signed char align ) { return n / align * align; }
284 unsigned char floor( unsigned char n, unsigned char align ) { return n / align * align; }
285 short int floor( short int n, short int align ) { return n / align * align; }
286 unsigned short int floor( unsigned short int n, unsigned short int align ) { return n / align * align; }
287 int floor( int n, int align ) { return n / align * align; }
288 unsigned int floor( unsigned int n, unsigned int align ) { return n / align * align; }
289 long int floor( long int n, long int align ) { return n / align * align; }
290 unsigned long int floor( unsigned long int n, unsigned long int align ) { return n / align * align; }
291 long long int floor( long long int n, long long int align ) { return n / align * align; }
292 unsigned long long int floor( unsigned long long int n, unsigned long long int align ) { return n / align * align; }
293
294 // forall( T | { T ?/?( T, T ); T ?*?( T, T ); } )
295 // T floor( T n, T align ) { return n / align * align; }
296
297 signed char ceiling_div( signed char n, char align ) { return (n + (align - 1)) / align; }
298 unsigned char ceiling_div( unsigned char n, unsigned char align ) { return (n + (align - 1)) / align; }
299 short int ceiling_div( short int n, short int align ) { return (n + (align - 1)) / align; }
300 unsigned short int ceiling_div( unsigned short int n, unsigned short int align ) { return (n + (align - 1)) / align; }
301 int ceiling_div( int n, int align ) { return (n + (align - 1)) / align; }
302 unsigned int ceiling_div( unsigned int n, unsigned int align ) { return (n + (align - 1)) / align; }
303 long int ceiling_div( long int n, long int align ) { return (n + (align - 1)) / align; }
304 unsigned long int ceiling_div( unsigned long int n, unsigned long int align ) { return (n + (align - 1)) / align; }
305 long long int ceiling_div( long long int n, long long int align ) { return (n + (align - 1)) / align; }
306 unsigned long long int ceiling_div( unsigned long long int n, unsigned long long int align ) { return (n + (align - 1)) / align; }
307
308 // forall( T | { T ?+?( T, T ); T ?-?( T, T ); T ?%?( T, T ); } )
309 // T ceiling_div( T n, T align ) { verify( is_pow2( align ) );return (n + (align - 1)) / align; }
310
311 // gcc notices the div/mod pair and saves both so only one div.
312 signed char ceiling( signed char n, signed char align ) { return floor( n + (n % align != 0 ? align - 1 : 0), align ); }
313 unsigned char ceiling( unsigned char n, unsigned char align ) { return floor( n + (n % align != 0 ? align - 1 : 0), align ); }
314 short int ceiling( short int n, short int align ) { return floor( n + (n % align != 0 ? align - 1 : 0), align ); }
315 unsigned short int ceiling( unsigned short int n, unsigned short int align ) { return floor( n + (n % align != 0 ? align - 1 : 0), align ); }
316 int ceiling( int n, int align ) { return floor( n + (n % align != 0 ? align - 1 : 0), align ); }
317 unsigned int ceiling( unsigned int n, unsigned int align ) { return floor( n + (n % align != 0 ? align - 1 : 0), align ); }
318 long int ceiling( long int n, long int align ) { return floor( n + (n % align != 0 ? align - 1 : 0), align ); }
319 unsigned long int ceiling( unsigned long int n, unsigned long int align ) { return floor( n + (n % align != 0 ? align - 1 : 0) , align); }
320 long long int ceiling( long long int n, long long int align ) { return floor( n + (n % align != 0 ? align - 1 : 0), align ); }
321 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 ); }
322
323 // forall( T | { void ?{}( T &, one_t ); T ?+?( T, T ); T ?-?( T, T ); T ?/?( T, T ); } )
324 // T ceiling( T n, T align ) { return return floor( n + (n % align != 0 ? align - 1 : 0), align ); *}
325
326 float floor( float x ) { return floorf( x ); }
327 // extern "C" { double floor( double ); }
328 long double floor( long double x ) { return floorl( x ); }
329
330 float ceil( float x ) { return ceilf( x ); }
331 // extern "C" { double ceil( double ); }
332 long double ceil( long double x ) { return ceill( x ); }
333
334 float trunc( float x ) { return truncf( x ); }
335 // extern "C" { double trunc( double ); }
336 long double trunc( long double x ) { return truncl( x ); }
337
338 float rint( float x ) { return rintf( x ); }
339 // extern "C" { double rint( double x ); }
340 long double rint( long double x ) { return rintl( x ); }
341 long int rint( float x ) { return lrintf( x ); }
342 long int rint( double x ) { return lrint( x ); }
343 long int rint( long double x ) { return lrintl( x ); }
344 long long int rint( float x ) { return llrintf( x ); }
345 long long int rint( double x ) { return llrint( x ); }
346 long long int rint( long double x ) { return llrintl( x ); }
347
348 long int lrint( float x ) { return lrintf( x ); }
349 // extern "C" { long int lrint( double ); }
350 long int lrint( long double x ) { return lrintl( x ); }
351 long long int llrint( float x ) { return llrintf( x ); }
352 // extern "C" { long long int llrint( double ); }
353 long long int llrint( long double x ) { return llrintl( x ); }
354
355 float nearbyint( float x ) { return nearbyintf( x ); }
356 // extern "C" { double nearbyint( double ); }
357 long double nearbyint( long double x ) { return nearbyintl( x ); }
358
359 float round( float x ) { return roundf( x ); }
360 // extern "C" { double round( double x ); }
361 long double round( long double x ) { return roundl( x ); }
362 long int round( float x ) { return lroundf( x ); }
363 long int round( double x ) { return lround( x ); }
364 long int round( long double x ) { return lroundl( x ); }
365 long long int round( float x ) { return llroundf( x ); }
366 long long int round( double x ) { return llround( x ); }
367 long long int round( long double x ) { return llroundl( x ); }
368
369 long int lround( float x ) { return lroundf( x ); }
370 // extern "C" { long int lround( double ); }
371 long int lround( long double x ) { return lroundl( x ); }
372 long long int llround( float x ) { return llroundf( x ); }
373 // extern "C" { long long int llround( double ); }
374 long long int llround( long double x ) { return llroundl( x ); }
375} // distribution
376
377//---------------------- Manipulation ----------------------
378
379static inline {
380 float copysign( float x, float y ) { return copysignf( x, y ); }
381 // extern "C" { double copysign( double, double ); }
382 long double copysign( long double x, long double y ) { return copysignl( x, y ); }
383
384 float frexp( float x, int * ip ) { return frexpf( x, ip ); }
385 // extern "C" { double frexp( double, int * ); }
386 long double frexp( long double x, int * ip ) { return frexpl( x, ip ); }
387
388 float ldexp( float x, int exp2 ) { return ldexpf( x, exp2 ); }
389 // extern "C" { double ldexp( double, int ); }
390 long double ldexp( long double x, int exp2 ) { return ldexpl( x, exp2 ); }
391
392 [ float, float ] modf( float x ) { float i; x = modff( x, &i ); return [ i, x ]; }
393 float modf( float x, float * i ) { return modff( x, i ); }
394 [ double, double ] modf( double x ) { double i; x = modf( x, &i ); return [ i, x ]; }
395 // extern "C" { double modf( double, double * ); }
396 [ long double, long double ] modf( long double x ) { long double i; x = modfl( x, &i ); return [ i, x ]; }
397 long double modf( long double x, long double * i ) { return modfl( x, i ); }
398
399 float nextafter( float x, float y ) { return nextafterf( x, y ); }
400 // extern "C" { double nextafter( double, double ); }
401 long double nextafter( long double x, long double y ) { return nextafterl( x, y ); }
402
403 float nexttoward( float x, long double y ) { return nexttowardf( x, y ); }
404 // extern "C" { double nexttoward( double, long double ); }
405 long double nexttoward( long double x, long double y ) { return nexttowardl( x, y ); }
406
407 float scalbn( float x, int exp ) { return scalbnf( x, exp ); }
408 // extern "C" { double scalbn( double, int ); }
409 long double scalbn( long double x, int exp ) { return scalbnl( x, exp ); }
410 float scalbn( float x, long int exp ) { return scalblnf( x, exp ); }
411 double scalbn( double x, long int exp ) { return scalbln( x, exp ); }
412 long double scalbn( long double x, long int exp ) { return scalblnl( x, exp ); }
413
414 float scalbln( float x, long int exp ) { return scalblnf( x, exp ); }
415 // extern "C" { double scalbln( double, long int ); }
416 long double scalbln( long double x, long int exp ) { return scalblnl( x, exp ); }
417} // distribution
418
419//---------------------------------------
420
421static inline {
422 forall( T | { void ?{}( T &, one_t ); T ?+?( T, T ); T ?-?( T, T );T ?*?( T, T ); } )
423 T lerp( T x, T y, T a ) { return x * ((T){1} - a) + y * a; }
424
425 forall( T | { void ?{}( T &, zero_t ); void ?{}( T &, one_t ); int ?<?( T, T ); } )
426 T step( T edge, T x ) { return x < edge ? (T){0} : (T){1}; }
427
428 forall( T | { void ?{}( T &, int ); T clamp( T, T, T ); T ?-?( T, T ); T ?*?( T, T ); T ?/?( T, T ); } )
429 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); }
430} // distribution
431
432static inline unsigned long long log2_u32_32(unsigned long long val) {
433 enum {
434 TABLE_BITS = 6,
435 TABLE_SIZE = (1 << TABLE_BITS) + 2,
436 };
437 // for(i; TABLE_SIZE) {
438 // table[i] = (unsigned long long)(log2(1.0 + i / pow(2, TABLE_BITS)) * pow(2, 32)));
439 // }
440 static const unsigned long long table[] = {
441 0x0000000000, 0x0005b9e5a1, 0x000b5d69ba, 0x0010eb389f,
442 0x001663f6fa, 0x001bc84240, 0x002118b119, 0x002655d3c4,
443 0x002b803473, 0x00309857a0, 0x00359ebc5b, 0x003a93dc98,
444 0x003f782d72, 0x00444c1f6b, 0x0049101eac, 0x004dc4933a,
445 0x005269e12f, 0x00570068e7, 0x005b888736, 0x006002958c,
446 0x00646eea24, 0x0068cdd829, 0x006d1fafdc, 0x007164beb4,
447 0x00759d4f80, 0x0079c9aa87, 0x007dea15a3, 0x0081fed45c,
448 0x0086082806, 0x008a064fd5, 0x008df988f4, 0x0091e20ea1,
449 0x0095c01a39, 0x009993e355, 0x009d5d9fd5, 0x00a11d83f4,
450 0x00a4d3c25e, 0x00a8808c38, 0x00ac241134, 0x00afbe7fa0,
451 0x00b3500472, 0x00b6d8cb53, 0x00ba58feb2, 0x00bdd0c7c9,
452 0x00c1404ead, 0x00c4a7ba58, 0x00c80730b0, 0x00cb5ed695,
453 0x00ceaecfea, 0x00d1f73f9c, 0x00d53847ac, 0x00d8720935,
454 0x00dba4a47a, 0x00ded038e6, 0x00e1f4e517, 0x00e512c6e5,
455 0x00e829fb69, 0x00eb3a9f01, 0x00ee44cd59, 0x00f148a170,
456 0x00f446359b, 0x00f73da38d, 0x00fa2f045e, 0x00fd1a708b,
457 0x0100000000, 0x0102dfca16,
458 };
459 _Static_assert((sizeof(table) / sizeof(table[0])) == TABLE_SIZE, "TABLE_SIZE should be accurate");
460 // starting from val = (2 ** i)*(1 + f) where 0 <= f < 1
461 // log identities mean log2(val) = log2((2 ** i)*(1 + f)) = log2(2**i) + log2(1+f)
462 //
463 // getting i is easy to do using builtin_clz (count leading zero)
464 //
465 // we want to calculate log2(1+f) independently to have a many bits of precision as possible.
466 // val = (2 ** i)*(1 + f) = 2 ** i + f * 2 ** i
467 // isolating f we get
468 // val - 2 ** i = f * 2 ** i
469 // (val - 2 ** i) / 2 ** i = f
470 //
471 // we want to interpolate from the table to get the values
472 // and compromise by doing quadratic interpolation (rather than higher degree interpolation)
473 //
474 // for the interpolation we want to shift everything the fist sample point
475 // so our parabola becomes x = 0
476 // this further simplifies the equations
477 //
478 // the consequence is that we need f in 2 forms:
479 // - finding the index of x0
480 // - finding the distance between f and x0
481 //
482 // since sample points are equidistant we can significantly simplify the equations
483
484 // get i
485 const unsigned long long bits = sizeof(val) * __CHAR_BIT__;
486 const unsigned long long lz = __builtin_clzl(val);
487 const unsigned long long i = bits - 1 - lz;
488
489 // get the fractinal part as a u32.32
490 const unsigned long long frac = (val << (lz + 1)) >> 32;
491
492 // get high order bits for the index into the table
493 const unsigned long long idx0 = frac >> (32 - TABLE_BITS);
494
495 // get the x offset, i.e., the difference between the first sample point and the actual fractional part
496 const long long udx = frac - (idx0 << (32 - TABLE_BITS));
497 /* paranoid */ verify((idx0 + 2) < TABLE_SIZE);
498
499 const long long y0 = table[idx0 + 0];
500 const long long y1 = table[idx0 + 1];
501 const long long y2 = table[idx0 + 2];
502
503 // from there we can quadraticly interpolate to get the data, using the lagrange polynomial
504 // normally it would look like:
505 // double r0 = y0 * ((x - x1) / (x0 - x1)) * ((x - x2) / (x0 - x2));
506 // double r1 = y1 * ((x - x0) / (x1 - x0)) * ((x - x2) / (x1 - x2));
507 // double r2 = y2 * ((x - x0) / (x2 - x0)) * ((x - x1) / (x2 - x1));
508 // but since the spacing between sample points is fixed, we can simplify it and extract common expressions
509 const long long f1 = (y1 - y0);
510 const long long f2 = (y2 - y0);
511 const long long a = f2 - (f1 * 2l);
512 const long long b = (f1 * 2l) - a;
513
514 // Now we can compute it in the form (ax + b)x + c (which avoid repeating steps)
515 long long sum = ((a*udx) >> (32 - TABLE_BITS)) + b;
516 sum = (sum*udx) >> (32 - TABLE_BITS + 1);
517 sum = y0 + sum;
518
519 return (i << 32) + (sum);
520}
521
522// Local Variables: //
523// mode: c //
524// tab-width: 4 //
525// End: //
Note: See TracBrowser for help on using the repository browser.