42 #ifndef OPENIMAGEIO_FMATH_H
43 #define OPENIMAGEIO_FMATH_H
50 #if defined(_MSC_VER) && _MSC_VER < 1600
51 typedef __int8 int8_t;
52 typedef __int16 int16_t;
53 typedef __int32 int32_t;
54 typedef __int64 int64_t;
55 typedef unsigned __int8 uint8_t;
56 typedef unsigned __int16 uint16_t;
58 typedef unsigned __int32 uint32_t;
59 typedef unsigned __int64 uint64_t;
67 #if defined(__FreeBSD__)
68 #include <sys/param.h>
79 # define M_PI 3.1415926535897932
85 # define M_PI_2 1.5707963267948966
91 # define M_TWO_PI (M_PI * 2.0)
97 # define M_1_PI 0.318309886183790671538
103 # define M_2_PI 0.63661977236758134
109 # define M_SQRT2 1.414135623730950
115 # define M_SQRT1_2 0.7071067811865475
121 # define M_LN2 0.6931471805599453
127 # define M_LN10 2.3025850929940457
133 #if (defined(_WIN32) || defined(__i386__) || defined(__x86_64__))
134 # define USE_INTEL_MATH_SHORTCUTS 1
135 # ifndef __LITTLE_ENDIAN__
136 # define __LITTLE_ENDIAN__ 1
137 # undef __BIG_ENDIAN__
145 #define HUGE_FLOAT ((float)1.0e38)
149 inline bool huge (
float f) {
return (f >= HUGE_FLOAT/2); }
153 #define UNINITIALIZED_FLOAT (- std::numeric_limits<float>::max())
166 return (x & (x-1)) == 0 && (x >= 0);
174 ispow2 (
unsigned int x)
179 return (x & (x-1)) == 0;
210 pow2rounddown (
int x)
220 return x & ~(x >> 1);
228 round_to_multiple (
int x,
int m)
230 return ((x + m - 1) / m) * m;
237 inline bool littleendian (
void)
239 #if defined(__BIG_ENDIAN__)
241 #elif defined(__LITTLE_ENDIAN__)
246 return *((
char *) &i);
254 inline bool bigendian (
void)
256 return ! littleendian();
266 swap_endian (T *f,
int len=1)
268 for (
char *c = (
char *) f; len--; c +=
sizeof(T)) {
269 if (
sizeof(T) == 2) {
270 std::swap (c[0], c[1]);
271 }
else if (
sizeof(T) == 4) {
272 std::swap (c[0], c[3]);
273 std::swap (c[1], c[2]);
274 }
else if (
sizeof(T) == 8) {
275 std::swap (c[0], c[7]);
276 std::swap (c[1], c[6]);
277 std::swap (c[2], c[5]);
278 std::swap (c[3], c[4]);
287 template <
class T=
float>
288 class EightBitConverter {
290 EightBitConverter () { init(); }
291 T operator() (
unsigned char c)
const {
return val[c]; }
295 float scale = 1.0f / 255.0f;
296 if (std::numeric_limits<T>::is_integer)
297 scale *= (float)std::numeric_limits<T>::max();
298 for (
int i = 0; i < 256; ++i)
299 val[i] = (T)(i * scale);
309 clamp (T a, T l, T h)
311 return (a < l)? l : ((a > h)? h : a);
319 clamped_mult32 (uint32_t a, uint32_t b)
321 const uint32_t Err = std::numeric_limits<uint32_t>::max();
322 uint64_t r = (uint64_t)a * (uint64_t)b;
323 return r < Err ? (uint32_t)r : Err;
331 clamped_mult64 (uint64_t a, uint64_t b)
335 return std::numeric_limits<uint64_t>::max();
343 template<
typename T,
typename U>
struct is_same {
static const bool value =
false; };
344 template<
typename T>
struct is_same<T,T> {
static const bool value =
true; };
351 template<
typename S,
typename D,
typename F>
352 D scaled_conversion (
const S &src, F scale, F min, F max)
354 if (std::numeric_limits<S>::is_signed) {
356 s += (s < 0 ? (F)-0.5 : (F)0.5);
357 return (D) clamp (s, min, max);
359 return (D) clamp ((F)src * scale + (F)0.5, min, max);
374 template<
typename S,
typename D>
375 void convert_type (
const S *src, D *dst,
size_t n, D _zero=0, D _one=1,
376 D _min=std::numeric_limits<D>::min(),
377 D _max=std::numeric_limits<D>::max())
379 if (is_same<S,D>::value) {
381 memcpy (dst, src, n*
sizeof(D));
385 F scale = std::numeric_limits<S>::is_integer ?
386 ((F)1.0)/std::numeric_limits<S>::max() : (F)1.0;
387 if (std::numeric_limits<D>::is_integer) {
393 for ( ; n >= 16; n -= 16) {
394 *dst++ = scaled_conversion<S,D,F> (*src++, scale, min, max);
395 *dst++ = scaled_conversion<S,D,F> (*src++, scale, min, max);
396 *dst++ = scaled_conversion<S,D,F> (*src++, scale, min, max);
397 *dst++ = scaled_conversion<S,D,F> (*src++, scale, min, max);
398 *dst++ = scaled_conversion<S,D,F> (*src++, scale, min, max);
399 *dst++ = scaled_conversion<S,D,F> (*src++, scale, min, max);
400 *dst++ = scaled_conversion<S,D,F> (*src++, scale, min, max);
401 *dst++ = scaled_conversion<S,D,F> (*src++, scale, min, max);
402 *dst++ = scaled_conversion<S,D,F> (*src++, scale, min, max);
403 *dst++ = scaled_conversion<S,D,F> (*src++, scale, min, max);
404 *dst++ = scaled_conversion<S,D,F> (*src++, scale, min, max);
405 *dst++ = scaled_conversion<S,D,F> (*src++, scale, min, max);
406 *dst++ = scaled_conversion<S,D,F> (*src++, scale, min, max);
407 *dst++ = scaled_conversion<S,D,F> (*src++, scale, min, max);
408 *dst++ = scaled_conversion<S,D,F> (*src++, scale, min, max);
409 *dst++ = scaled_conversion<S,D,F> (*src++, scale, min, max);
412 *dst++ = scaled_conversion<S,D,F> (*src++, scale, min, max);
417 for ( ; n >= 16; n -= 16) {
418 *dst++ = (D)((*src++) * scale);
419 *dst++ = (D)((*src++) * scale);
420 *dst++ = (D)((*src++) * scale);
421 *dst++ = (D)((*src++) * scale);
422 *dst++ = (D)((*src++) * scale);
423 *dst++ = (D)((*src++) * scale);
424 *dst++ = (D)((*src++) * scale);
425 *dst++ = (D)((*src++) * scale);
426 *dst++ = (D)((*src++) * scale);
427 *dst++ = (D)((*src++) * scale);
428 *dst++ = (D)((*src++) * scale);
429 *dst++ = (D)((*src++) * scale);
430 *dst++ = (D)((*src++) * scale);
431 *dst++ = (D)((*src++) * scale);
432 *dst++ = (D)((*src++) * scale);
433 *dst++ = (D)((*src++) * scale);
436 *dst++ = (D)((*src++) * scale);
447 template<
typename S,
typename D>
448 D convert_type (
const S &src)
450 if (is_same<S,D>::value) {
455 F scale = std::numeric_limits<S>::is_integer ?
456 ((F)1.0)/std::numeric_limits<S>::max() : (F)1.0;
457 if (std::numeric_limits<D>::is_integer) {
460 F min = (F) std::numeric_limits<D>::min();
461 F max = (F) std::numeric_limits<D>::max();
463 return scaled_conversion<S,D,F> (src, scale, min, max);
467 return (D)((F)src * scale);
484 template<
unsigned int FROM_BITS,
unsigned int TO_BITS>
485 inline unsigned int bit_range_convert(
unsigned int in) {
486 unsigned int out = 0;
487 int shift = TO_BITS - FROM_BITS;
488 for (; shift > 0; shift -= FROM_BITS)
498 bit_range_convert(
unsigned int in,
unsigned int FROM_BITS,
unsigned int TO_BITS)
500 unsigned int out = 0;
501 int shift = TO_BITS - FROM_BITS;
502 for (; shift > 0; shift -= FROM_BITS)
513 template<
typename I,
typename E>
515 DataProxy (I &data) : m_data(data) { }
516 E operator= (E newval) { m_data = convert_type<E,I>(newval);
return newval; }
517 operator E ()
const {
return convert_type<I,E>(m_data); }
519 DataProxy& operator = (
const DataProxy&);
528 template<
typename I,
typename E>
529 struct ConstDataProxy {
530 ConstDataProxy (
const I &data) : m_data(data) { }
531 operator E ()
const {
return convert_type<E,I>(*m_data); }
541 template<
typename I,
typename E>
542 struct DataArrayProxy {
543 DataArrayProxy (I *data=NULL) : m_data(data) { }
544 E operator* ()
const {
return convert_type<I,E>(*m_data); }
545 E operator[] (
int i)
const {
return convert_type<I,E>(m_data[i]); }
546 DataProxy<I,E> operator[] (
int i) {
return DataProxy<I,E> (m_data[i]); }
547 void set (I *data) { m_data = data; }
548 I *
get ()
const {
return m_data; }
549 const DataArrayProxy<I,E> & operator+= (
int i) {
550 m_data += i;
return *
this;
561 template<
typename I,
typename E>
562 struct ConstDataArrayProxy {
563 ConstDataArrayProxy (
const I *data=NULL) : m_data(data) { }
564 E operator* ()
const {
return convert_type<I,E>(*m_data); }
565 E operator[] (
int i)
const {
return convert_type<I,E>(m_data[i]); }
566 void set (
const I *data) { m_data = data; }
567 const I *
get ()
const {
return m_data; }
568 const ConstDataArrayProxy<I,E> & operator+= (
int i) {
569 m_data += i;
return *
this;
580 template <
class T,
class Q>
582 lerp (T v0, T v1, Q x)
585 return v0*(Q(1)-x) + v1*x;
593 template <
class T,
class Q>
595 bilerp (T v0, T v1, T v2, T v3, Q s, Q t)
599 return (T) (((Q)1-t)*(v0*s1 + v1*s) + t*(v2*s1 + v3*s));
608 template <
class T,
class Q>
610 bilerp (
const T *v0,
const T *v1,
611 const T *v2,
const T *v3,
612 Q s, Q t,
int n, T *result)
616 for (
int i = 0; i < n; ++i)
617 result[i] = (T) (t1*(v0[i]*s1 + v1[i]*s) + t*(v2[i]*s1 + v3[i]*s));
627 template <
class T,
class Q>
629 bilerp_mad (
const T *v0,
const T *v1,
630 const T *v2,
const T *v3,
631 Q s, Q t, Q scale,
int n, T *result)
635 for (
int i = 0; i < n; ++i)
636 result[i] += (T) (scale * (t1*(v0[i]*s1 + v1[i]*s) +
637 t*(v2[i]*s1 + v3[i]*s)));
645 template <
class T,
class Q>
647 trilerp (T v0, T v1, T v2, T v3, T v4, T v5, T v6, T v7, Q s, Q t, Q r)
653 return (T) (r1*(t1*(v0*s1 + v1*s) + t*(v2*s1 + v3*s)) +
654 r*(t1*(v4*s1 + v5*s) + t*(v6*s1 + v7*s)));
663 template <
class T,
class Q>
665 trilerp (
const T *v0,
const T *v1,
const T *v2,
const T *v3,
666 const T *v4,
const T *v5,
const T *v6,
const T *v7,
667 Q s, Q t, Q r,
int n, T *result)
672 for (
int i = 0; i < n; ++i)
673 result[i] = (T) (r1*(t1*(v0[i]*s1 + v1[i]*s) + t*(v2[i]*s1 + v3[i]*s)) +
674 r*(t1*(v4[i]*s1 + v5[i]*s) + t*(v6[i]*s1 + v7[i]*s)));
684 template <
class T,
class Q>
686 trilerp_mad (
const T *v0,
const T *v1,
const T *v2,
const T *v3,
687 const T *v4,
const T *v5,
const T *v6,
const T *v7,
688 Q s, Q t, Q r, Q scale,
int n, T *result)
691 bilerp_mad (v0, v1, v2, v3, s, t, scale*r1, n, result);
692 bilerp_mad (v4, v5, v6, v7, s, t, scale*r, n, result);
701 RoundToInt (
double val)
703 #ifdef USE_INTEL_MATH_SHORTCUTS
704 union {
int i;
double d; } myunion;
706 const double doublemagic = double (6755399441055744.0);
708 myunion.d += doublemagic;
717 RoundToInt (
float val)
719 #ifdef USE_INTEL_MATH_SHORTCUTS
720 return RoundToInt (
double(val));
731 FloorToInt (
double val)
733 #ifdef USE_INTEL_MATH_SHORTCUTS
734 const double doublemagicdelta = (1.5e-8);
736 const double doublemagicroundeps = (0.5f-doublemagicdelta);
738 return RoundToInt (val - doublemagicroundeps);
740 return (
int) floor (val);
746 FloorToInt (
float val)
748 #ifdef USE_INTEL_MATH_SHORTCUTS
749 return FloorToInt (
double(val));
751 return (
int) floorf (val);
760 CeilToInt (
double val)
762 #ifdef USE_INTEL_MATH_SHORTCUTS
763 const double doublemagicdelta = (1.5e-8);
765 const double doublemagicroundeps = (0.5f-doublemagicdelta);
767 return RoundToInt (val + doublemagicroundeps);
769 return (
int) ceil (val);
774 CeilToInt (
float val)
776 #ifdef USE_INTEL_MATH_SHORTCUTS
777 return CeilToInt (
double(val));
779 return (
int) ceilf (val);
787 FloatToInt (
double val)
789 #ifdef USE_INTEL_MATH_SHORTCUTS
790 return (val<0) ? CeilToInt(val) : FloorToInt(val);
798 FloatToInt (
float val)
800 #ifdef USE_INTEL_MATH_SHORTCUTS
801 return FloatToInt (
double(val));
815 floorfrac (
float x,
int *xint)
820 int i = (int) x - (x < 0.0f ? 1 : 0);
823 int i = FloorToInt (x);
826 return x -
static_cast<float>(i);
833 inline float radians (
float deg) {
return deg * (float)(M_PI / 180.0); }
837 inline float degrees (
float rad) {
return rad * (float)(180.0 / M_PI); }
842 inline float fast_expf(
float x)
844 #if defined(__x86_64__) && defined(__GNU_LIBRARY__) && defined(__GLIBC__ ) && defined(__GLIBC_MINOR__) && __GLIBC__ <= 2 && __GLIBC_MINOR__ < 16
847 return static_cast<float>(std::exp(static_cast<double>(x)));
857 #define hypotf _hypotf
858 #define copysign(x,y) _copysign(x,y)
859 #define copysignf(x,y) copysign(x,y)
861 #define isnan(x) _isnan(x)
862 #define isfinite(x) _finite(x)
865 #define M_E 2.71828182845904523536
866 #define M_LOG2E 1.44269504088896340736
872 #define M_PI_4 0.785398163397448309616
881 return floor (val + 0.5);
886 return static_cast<float>(round (val));
891 inline int isinf (T x) {
892 return (isfinite(x)||isnan(x)) ? 0 :
static_cast<int>(copysign(T(1.0), x));
898 return logf (val)/
static_cast<float>(M_LN2);
905 return logf(val)/logf(FLT_RADIX);
911 return exp( val*log(2.0f) );
919 if (fabs(val) < 1e-5)
920 return val + 0.5*val*val;
922 return exp(val) - 1.0;
928 return (
float)expm1(val);
936 double a1 = 0.254829592;
937 double a2 = -0.284496736;
938 double a3 = 1.421413741;
939 double a4 = -1.453152027;
940 double a5 = 1.061405429;
941 double p = 0.3275911;
950 double t = 1.0/(1.0 + p*x);
951 double y = 1.0 - (((((a5*t + a4)*t) + a3)*t + a2)*t + a1)*t*exp(-x*x);
959 return (
float)erf(val);
965 return 1.0 - erf(val);
971 return (
float)erfc(val);
977 return (
float)(int)val;
1002 #if (defined(__FreeBSD__) && (__FreeBSD_version < 803000))
1006 return logf (val)/
static_cast<float>(M_LN2);
1020 float_to_rational (
float f,
unsigned int &num,
unsigned int &den)
1025 }
else if ((
int)(1.0/f) == (1.0/f)) {
1031 while (fabsf(f-static_cast<float>(num)) > 0.00001f && den < 1000000) {
1047 float_to_rational (
float f,
int &num,
int &den)
1050 float_to_rational (fabsf(f), n, d);
1051 num = (f >= 0) ? (
int)n : -(int)n;
1057 sincos(
float x,
float* sine,
float* cosine)
1059 #if defined(__GNUC__) && defined(__linux__) && !defined(__clang__)
1060 __builtin_sincosf(x, sine, cosine);
1062 *sine = std::sin(x);
1063 *cosine = std::cos(x);
1068 sincos(
double x,
double* sine,
double* cosine)
1070 #if defined(__GNUC__) && defined(__linux__) && !defined(__clang__)
1071 __builtin_sincos(x, sine, cosine);
1073 *sine = std::sin(x);
1074 *cosine = std::cos(x);
1083 safe_asinf (
float x)
1085 if (x >= 1.0f)
return static_cast<float>(M_PI)/2;
1086 if (x <= -1.0f)
return -
static_cast<float>(M_PI)/2;
1087 return std::asin (x);
1094 safe_acosf (
float x) {
1095 if (x >= 1.0f)
return 0.0f;
1096 if (x <= -1.0f)
return static_cast<float>(M_PI);
1097 return std::acos (x);
1104 safe_sqrt (
double x)
1106 return (x > 0.0) ? sqrt(x) : 0.0;
1110 safe_sqrtf (
float x)
1112 return (x > 0.0f) ? sqrtf(x) : 0.0f;
1116 inline float safe_inversesqrt (
float x) {
1117 return (x > 0.0f) ? 1.0f/std::sqrt(x) : 0.0f;
1120 inline float safe_log (
float x) {
1121 return (x > 0.0f) ? logf(x) : -std::numeric_limits<float>::max();
1124 inline float safe_log2(
float x) {
1125 return (x > 0.0f) ? log2f(x) : -std::numeric_limits<float>::max();
1128 inline float safe_log10(
float x) {
1129 return (x > 0.0f) ? log10f(x) : -std::numeric_limits<float>::max();
1132 inline float safe_logb (
float x) {
1133 return (x != 0.0f) ? logbf(x) : -std::numeric_limits<float>::max();
1148 template<
class T,
class Func>
1149 T invert (Func &func, T y, T xmin=0.0, T xmax=1.0,
1150 int maxiters=32, T eps=1.0e-6,
bool *brack=0)
1156 T v0 = func(xmin), v1 = func(xmax);
1158 bool increasing = (v0 < v1);
1159 T vmin = increasing ? v0 : v1;
1160 T vmax = increasing ? v1 : v0;
1161 bool bracketed = (y >= vmin && y <= vmax);
1167 return ((y < vmin) == increasing) ? xmin : xmax;
1169 if (fabs(v0-v1) < eps)
1171 int rfiters = (3*maxiters)/4;
1172 for (
int iters = 0; iters < maxiters; ++iters) {
1174 if (iters < rfiters) {
1177 if (t <= T(0) || t >= T(1))
1182 x = lerp (xmin, xmax, t);
1184 if ((v < y) == increasing) {
1189 if (fabs(xmax-xmin) < eps || fabs(v-y) < eps)
1200 #endif // OPENIMAGEIO_FMATH_H