OpenImageIO
 All Classes Files Friends Macros Pages
fmath.h
Go to the documentation of this file.
1 /*
2  Copyright 2008 Larry Gritz and the other authors and contributors.
3  All Rights Reserved.
4 
5  Redistribution and use in source and binary forms, with or without
6  modification, are permitted provided that the following conditions are
7  met:
8  * Redistributions of source code must retain the above copyright
9  notice, this list of conditions and the following disclaimer.
10  * Redistributions in binary form must reproduce the above copyright
11  notice, this list of conditions and the following disclaimer in the
12  documentation and/or other materials provided with the distribution.
13  * Neither the name of the software's owners nor the names of its
14  contributors may be used to endorse or promote products derived from
15  this software without specific prior written permission.
16  THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
17  "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
18  LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
19  A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
20  OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
21  SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
22  LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
23  DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
24  THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
25  (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
26  OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
27 
28  (This is the Modified BSD License)
29 
30  A few bits here are based upon code from NVIDIA that was also released
31  under the same modified BSD license, and marked as:
32  Copyright 2004 NVIDIA Corporation. All Rights Reserved.
33 */
34 
35 
41 
42 #ifndef OPENIMAGEIO_FMATH_H
43 #define OPENIMAGEIO_FMATH_H
44 
45 #include <cmath>
46 #include <limits>
47 #include <typeinfo>
48 #include <algorithm>
49 
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;
57 # ifndef _UINT64_T
58  typedef unsigned __int32 uint32_t;
59  typedef unsigned __int64 uint64_t;
60 # define _UINT32_T
61 # define _UINT64_T
62 # endif
63 #else
64 # include <stdint.h>
65 #endif
66 
67 #if defined(__FreeBSD__)
68 #include <sys/param.h>
69 #endif
70 
71 #include "version.h"
72 
73 OIIO_NAMESPACE_ENTER
74 {
75 
76 #ifndef M_PI
77 
78 
79 # define M_PI 3.1415926535897932
80 #endif
81 
82 #ifndef M_PI_2
83 
84 
85 # define M_PI_2 1.5707963267948966
86 #endif
87 
88 #ifndef M_TWO_PI
89 
90 
91 # define M_TWO_PI (M_PI * 2.0)
92 #endif
93 
94 #ifndef M_1_PI
95 
96 
97 # define M_1_PI 0.318309886183790671538
98 #endif
99 
100 #ifndef M_2_PI
101 
102 
103 # define M_2_PI 0.63661977236758134
104 #endif
105 
106 #ifndef M_SQRT2
107 
108 
109 # define M_SQRT2 1.414135623730950
110 #endif
111 
112 #ifndef M_SQRT1_2
113 
114 
115 # define M_SQRT1_2 0.7071067811865475
116 #endif
117 
118 #ifndef M_LN2
119 
120 
121 # define M_LN2 0.6931471805599453
122 #endif
123 
124 #ifndef M_LN10
125 
126 
127 # define M_LN10 2.3025850929940457
128 #endif
129 
130 
131 
132 // Some stuff we know easily if we're running on an Intel chip
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__
138 # endif
139 #endif
140 
141 
142 
145 #define HUGE_FLOAT ((float)1.0e38)
146 
149 inline bool huge (float f) { return (f >= HUGE_FLOAT/2); }
150 
153 #define UNINITIALIZED_FLOAT (- std::numeric_limits<float>::max())
154 
155 
156 
157 
160 inline bool
161 ispow2 (int x)
162 {
163  // Numerous references for this bit trick are on the web. The
164  // principle is that x is a power of 2 <=> x == 1<<b <=> x-1 is
165  // all 1 bits for bits < b.
166  return (x & (x-1)) == 0 && (x >= 0);
167 }
168 
169 
170 
173 inline bool
174 ispow2 (unsigned int x)
175 {
176  // Numerous references for this bit trick are on the web. The
177  // principle is that x is a power of 2 <=> x == 1<<b <=> x-1 is
178  // all 1 bits for bits < b.
179  return (x & (x-1)) == 0;
180 }
181 
182 
183 
186 inline int
187 pow2roundup (int x)
188 {
189  // Here's a version with no loops.
190  if (x < 0)
191  return 0;
192  // Subtract 1, then round up to a power of 2, that way if we are
193  // already a power of 2, we end up with the same number.
194  --x;
195  // Make all bits past the first 1 also be 1, i.e. 0001xxxx -> 00011111
196  x |= x >> 1;
197  x |= x >> 2;
198  x |= x >> 4;
199  x |= x >> 8;
200  x |= x >> 16;
201  // Now we have 2^n-1, by adding 1, we make it a power of 2 again
202  return x+1;
203 }
204 
205 
206 
209 inline int
210 pow2rounddown (int x)
211 {
212  // Make all bits past the first 1 also be 1, i.e. 0001xxxx -> 00011111
213  x |= x >> 1;
214  x |= x >> 2;
215  x |= x >> 4;
216  x |= x >> 8;
217  x |= x >> 16;
218  // Strip off all but the high bit, i.e. 00011111 -> 00010000
219  // That's the power of two <= the original x
220  return x & ~(x >> 1);
221 }
222 
223 
224 
227 inline int
228 round_to_multiple (int x, int m)
229 {
230  return ((x + m - 1) / m) * m;
231 }
232 
233 
234 
237 inline bool littleendian (void)
238 {
239 #if defined(__BIG_ENDIAN__)
240  return false;
241 #elif defined(__LITTLE_ENDIAN__)
242  return true;
243 #else
244  // Otherwise, do something quick to compute it
245  int i = 1;
246  return *((char *) &i);
247 #endif
248 }
249 
250 
251 
254 inline bool bigendian (void)
255 {
256  return ! littleendian();
257 }
258 
259 
260 
264 template<class T>
265 inline void
266 swap_endian (T *f, int len=1)
267 {
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]);
279  }
280  }
281 }
282 
283 
284 
287 template <class T=float>
288 class EightBitConverter {
289 public:
290  EightBitConverter () { init(); }
291  T operator() (unsigned char c) const { return val[c]; }
292 private:
293  T val[256];
294  void init () {
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);
300  }
301 };
302 
303 
304 
307 template <class T>
308 inline T
309 clamp (T a, T l, T h)
310 {
311  return (a < l)? l : ((a > h)? h : a);
312 }
313 
314 
315 
318 inline uint32_t
319 clamped_mult32 (uint32_t a, uint32_t b)
320 {
321  const uint32_t Err = std::numeric_limits<uint32_t>::max();
322  uint64_t r = (uint64_t)a * (uint64_t)b; // Multiply into a bigger int
323  return r < Err ? (uint32_t)r : Err;
324 }
325 
326 
327 
330 inline uint64_t
331 clamped_mult64 (uint64_t a, uint64_t b)
332 {
333  uint64_t ab = a*b;
334  if (b && ab/b != a)
335  return std::numeric_limits<uint64_t>::max();
336  else
337  return ab;
338 }
339 
340 
341 
342 // Helper template to let us tell if two types are the same.
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; };
345 
346 
347 
351 template<typename S, typename D, typename F>
352 D scaled_conversion (const S &src, F scale, F min, F max)
353 {
354  if (std::numeric_limits<S>::is_signed) {
355  F s = src * scale;
356  s += (s < 0 ? (F)-0.5 : (F)0.5);
357  return (D) clamp (s, min, max);
358  } else {
359  return (D) clamp ((F)src * scale + (F)0.5, min, max);
360  }
361 }
362 
363 
364 
371 //
372 // FIXME: make table-based specializations for common types with only a
373 // few possible src values (like unsigned char -> float).
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())
378 {
379  if (is_same<S,D>::value) {
380  // They must be the same type. Just memcpy.
381  memcpy (dst, src, n*sizeof(D));
382  return;
383  }
384  typedef double F;
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) {
388  // Converting to an integer-like type.
389  F min = (F)_min; // std::numeric_limits<D>::min();
390  F max = (F)_max; // std::numeric_limits<D>::max();
391  scale *= _max;
392  // Unroll loop for speed
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);
410  }
411  while (n--)
412  *dst++ = scaled_conversion<S,D,F> (*src++, scale, min, max);
413  } else {
414  // Converting to a float-like type, so we don't need to remap
415  // the range
416  // Unroll loop for speed
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);
434  }
435  while (n--)
436  *dst++ = (D)((*src++) * scale);
437  }
438 }
439 
440 
441 
447 template<typename S, typename D>
448 D convert_type (const S &src)
449 {
450  if (is_same<S,D>::value) {
451  // They must be the same type. Just return it.
452  return (D)src;
453  }
454  typedef double F;
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) {
458  // Converting to an integer-like type.
459  typedef double F;
460  F min = (F) std::numeric_limits<D>::min();
461  F max = (F) std::numeric_limits<D>::max();
462  scale *= max;
463  return scaled_conversion<S,D,F> (src, scale, min, max);
464  } else {
465  // Converting to a float-like type, so we don't need to remap
466  // the range
467  return (D)((F)src * scale);
468  }
469 }
470 
471 
472 
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)
489  out |= in << shift;
490  out |= in >> -shift;
491  return out;
492 }
493 
494 
495 
496 // non-templated version. Slow but general
497 inline unsigned int
498 bit_range_convert(unsigned int in, unsigned int FROM_BITS, unsigned int TO_BITS)
499 {
500  unsigned int out = 0;
501  int shift = TO_BITS - FROM_BITS;
502  for (; shift > 0; shift -= FROM_BITS)
503  out |= in << shift;
504  out |= in >> -shift;
505  return out;
506 }
507 
508 
509 
513 template<typename I, typename E>
514 struct DataProxy {
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); }
518 private:
519  DataProxy& operator = (const DataProxy&); // Do not implement
520  I &m_data;
521 };
522 
523 
524 
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); }
532 private:
533  const I &m_data;
534 };
535 
536 
537 
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;
551  }
552 private:
553  I *m_data;
554 };
555 
556 
557 
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;
570  }
571 private:
572  const I *m_data;
573 };
574 
575 
576 
577 
580 template <class T, class Q>
581 inline T
582 lerp (T v0, T v1, Q x)
583 {
584  // NOTE: a*(1-x) + b*x is much more numerically stable than a+x*(b-a)
585  return v0*(Q(1)-x) + v1*x;
586 }
587 
588 
589 
593 template <class T, class Q>
594 inline T
595 bilerp (T v0, T v1, T v2, T v3, Q s, Q t)
596 {
597  // NOTE: a*(t-1) + b*t is much more numerically stable than a+t*(b-a)
598  Q s1 = (Q)1 - s;
599  return (T) (((Q)1-t)*(v0*s1 + v1*s) + t*(v2*s1 + v3*s));
600 }
601 
602 
603 
608 template <class T, class Q>
609 inline void
610 bilerp (const T *v0, const T *v1,
611  const T *v2, const T *v3,
612  Q s, Q t, int n, T *result)
613 {
614  Q s1 = (Q)1 - s;
615  Q t1 = (Q)1 - t;
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));
618 }
619 
620 
621 
627 template <class T, class Q>
628 inline void
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)
632 {
633  Q s1 = (Q)1 - s;
634  Q t1 = (Q)1 - t;
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)));
638 }
639 
640 
641 
645 template <class T, class Q>
646 inline T
647 trilerp (T v0, T v1, T v2, T v3, T v4, T v5, T v6, T v7, Q s, Q t, Q r)
648 {
649  // NOTE: a*(t-1) + b*t is much more numerically stable than a+t*(b-a)
650  Q s1 = (Q)1 - s;
651  Q t1 = (Q)1 - t;
652  Q r1 = (Q)1 - 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)));
655 }
656 
657 
658 
663 template <class T, class Q>
664 inline void
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)
668 {
669  Q s1 = (Q)1 - s;
670  Q t1 = (Q)1 - t;
671  Q r1 = (Q)1 - r;
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)));
675 }
676 
677 
678 
684 template <class T, class Q>
685 inline void
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)
689 {
690  Q r1 = (Q)1 - r;
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);
693 }
694 
695 
696 
700 inline int
701 RoundToInt (double val)
702 {
703 #ifdef USE_INTEL_MATH_SHORTCUTS
704  union { int i; double d; } myunion;
705  myunion.d = val;
706  const double doublemagic = double (6755399441055744.0);
707  // 2^52 * 1.5, uses limited precisicion to floor
708  myunion.d += doublemagic;
709  return myunion.i;
710 #else
711  return round (val);
712 #endif
713 }
714 
715 
716 inline int
717 RoundToInt (float val)
718 {
719 #ifdef USE_INTEL_MATH_SHORTCUTS
720  return RoundToInt (double(val));
721 #else
722  return roundf (val);
723 #endif
724 }
725 
726 
730 inline int
731 FloorToInt (double val)
732 {
733 #ifdef USE_INTEL_MATH_SHORTCUTS
734  const double doublemagicdelta = (1.5e-8);
735  // almost .5f = .5f + 1e^(number of exp bit)
736  const double doublemagicroundeps = (0.5f-doublemagicdelta);
737  // almost .5f = .5f - 1e^(number of exp bit)
738  return RoundToInt (val - doublemagicroundeps);
739 #else
740  return (int) floor (val);
741 #endif
742 }
743 
744 
745 inline int
746 FloorToInt (float val)
747 {
748 #ifdef USE_INTEL_MATH_SHORTCUTS
749  return FloorToInt (double(val));
750 #else
751  return (int) floorf (val);
752 #endif
753 }
754 
755 
759 inline int
760 CeilToInt (double val)
761 {
762 #ifdef USE_INTEL_MATH_SHORTCUTS
763  const double doublemagicdelta = (1.5e-8);
764  // almost .5f = .5f + 1e^(number of exp bit)
765  const double doublemagicroundeps = (0.5f-doublemagicdelta);
766  // almost .5f = .5f - 1e^(number of exp bit)
767  return RoundToInt (val + doublemagicroundeps);
768 #else
769  return (int) ceil (val);
770 #endif
771 }
772 
773 inline int
774 CeilToInt (float val)
775 {
776 #ifdef USE_INTEL_MATH_SHORTCUTS
777  return CeilToInt (double(val));
778 #else
779  return (int) ceilf (val);
780 #endif
781 }
782 
786 inline int
787 FloatToInt (double val)
788 {
789 #ifdef USE_INTEL_MATH_SHORTCUTS
790  return (val<0) ? CeilToInt(val) : FloorToInt(val);
791 #else
792  return (int) val;
793 #endif
794 }
795 
796 
797 inline int
798 FloatToInt (float val)
799 {
800 #ifdef USE_INTEL_MATH_SHORTCUTS
801  return FloatToInt (double(val));
802 #else
803  return (int) val;
804 #endif
805 }
806 
807 
808 
809 
814 inline float
815 floorfrac (float x, int *xint)
816 {
817  // Find the greatest whole number <= x. This cast is faster than
818  // calling floorf.
819 #if 1
820  int i = (int) x - (x < 0.0f ? 1 : 0);
821 #else
822  // Why isn't this faster than the cast?
823  int i = FloorToInt (x);
824 #endif
825  *xint = i;
826  return x - static_cast<float>(i); // Return the fraction left over
827 }
828 
829 
830 
833 inline float radians (float deg) { return deg * (float)(M_PI / 180.0); }
834 
837 inline float degrees (float rad) { return rad * (float)(180.0 / M_PI); }
838 
839 
840 
842 inline float fast_expf(float x)
843 {
844 #if defined(__x86_64__) && defined(__GNU_LIBRARY__) && defined(__GLIBC__ ) && defined(__GLIBC_MINOR__) && __GLIBC__ <= 2 && __GLIBC_MINOR__ < 16
845  // On x86_64, versions of glibc < 2.16 have an issue where expf is
846  // much slower than the double version. This was fixed in glibc 2.16.
847  return static_cast<float>(std::exp(static_cast<double>(x)));
848 #else
849  return std::exp(x);
850 #endif
851 }
852 
853 
854 
855 #ifdef _WIN32
856 // Windows doesn't define these functions from math.h
857 #define hypotf _hypotf
858 #define copysign(x,y) _copysign(x,y)
859 #define copysignf(x,y) copysign(x,y)
860 #ifdef _MSC_VER
861 #define isnan(x) _isnan(x)
862 #define isfinite(x) _finite(x)
863 #endif
864 
865 #define M_E 2.71828182845904523536
866 #define M_LOG2E 1.44269504088896340736
867 //#define M_LOG10E 0.434294481903251827651
868 //#define M_LN2 0.693147180559945309417
869 //#define M_LN10 2.30258509299404568402
870 //#define M_PI 3.14159265358979323846
871 //#define M_PI_2 1.57079632679489661923
872 #define M_PI_4 0.785398163397448309616
873 //#define M_1_PI 0.318309886183790671538
874 //#define M_2_PI 0.636619772367581343076
875 //#define M_2_SQRTPI 1.12837916709551257390
876 //#define M_SQRT2 1.41421356237309504880
877 //#define M_SQRT1_2 0.707106781186547524401
878 
879 inline double
880 round (float val) {
881  return floor (val + 0.5);
882 }
883 
884 inline float
885 roundf (float val) {
886  return static_cast<float>(round (val));
887 }
888 
889 #ifdef _MSC_VER
890 template<class T>
891 inline int isinf (T x) {
892  return (isfinite(x)||isnan(x)) ? 0 : static_cast<int>(copysign(T(1.0), x));
893 }
894 #endif
895 
896 inline float
897 log2f (float val) {
898  return logf (val)/static_cast<float>(M_LN2);
899 }
900 
901 
902 inline float
903 logbf(float val) {
904  // please see http://www.kernel.org/doc/man-pages/online/pages/man3/logb.3.html
905  return logf(val)/logf(FLT_RADIX);
906 }
907 
908 inline float
909 exp2f(float val) {
910  // 2^val = e^(val*ln(2))
911  return exp( val*log(2.0f) );
912 }
913 
914 // from http://www.johndcook.com/cpp_expm1.html
915 inline double
916 expm1(double val)
917 {
918  // exp(x) - 1 without loss of precision for small values of x.
919  if (fabs(val) < 1e-5)
920  return val + 0.5*val*val;
921  else
922  return exp(val) - 1.0;
923 }
924 
925 inline float
926 expm1f(float val)
927 {
928  return (float)expm1(val);
929 }
930 
931 // from http://www.johndcook.com/cpp_erf.html
932 inline double
933 erf(double x)
934 {
935  // constants
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;
942 
943  // Save the sign of x
944  int sign = 1;
945  if (x < 0)
946  sign = -1;
947  x = fabs(x);
948 
949  // A&S formula 7.1.26
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);
952 
953  return sign*y;
954 }
955 
956 inline float
957 erff(float val)
958 {
959  return (float)erf(val);
960 }
961 
962 inline double
963 erfc(double val)
964 {
965  return 1.0 - erf(val);
966 }
967 
968 inline float
969 erfcf(float val)
970 {
971  return (float)erfc(val);
972 }
973 
974 inline float
975 truncf(float val)
976 {
977  return (float)(int)val;
978 }
979 
980 using OIIO::roundf;
981 using OIIO::truncf;
982 using OIIO::expm1f;
983 using OIIO::erff;
984 using OIIO::erfcf;
985 using OIIO::log2f;
986 using OIIO::logbf;
987 using OIIO::exp2f;
988 
989 #endif /* _WIN32 */
990 
991 
992 // Some systems have isnan, isinf and isfinite in the std namespace.
993 #ifndef _MSC_VER
994  using std::isnan;
995  using std::isinf;
996  using std::isfinite;
997 #endif
998 
999 
1000 
1001 // Functions missing from FreeBSD
1002 #if (defined(__FreeBSD__) && (__FreeBSD_version < 803000))
1003 
1004 inline float
1005 log2f (float val) {
1006  return logf (val)/static_cast<float>(M_LN2);
1007 }
1008 
1009 using OIIO::log2f;
1010 #endif
1011 
1012 
1013 
1019 inline void
1020 float_to_rational (float f, unsigned int &num, unsigned int &den)
1021 {
1022  if (f <= 0) { // Trivial case of zero, and handle all negative values
1023  num = 0;
1024  den = 1;
1025  } else if ((int)(1.0/f) == (1.0/f)) { // Exact results for perfect inverses
1026  num = 1;
1027  den = (int)f;
1028  } else {
1029  num = (int)f;
1030  den = 1;
1031  while (fabsf(f-static_cast<float>(num)) > 0.00001f && den < 1000000) {
1032  den *= 10;
1033  f *= 10;
1034  num = (int)f;
1035  }
1036  }
1037 }
1038 
1039 
1040 
1046 inline void
1047 float_to_rational (float f, int &num, int &den)
1048 {
1049  unsigned int n, d;
1050  float_to_rational (fabsf(f), n, d);
1051  num = (f >= 0) ? (int)n : -(int)n;
1052  den = (int) d;
1053 }
1054 
1055 
1056 inline void
1057 sincos(float x, float* sine, float* cosine)
1058 {
1059 #if defined(__GNUC__) && defined(__linux__) && !defined(__clang__)
1060  __builtin_sincosf(x, sine, cosine);
1061 #else
1062  *sine = std::sin(x);
1063  *cosine = std::cos(x);
1064 #endif
1065 }
1066 
1067 inline void
1068 sincos(double x, double* sine, double* cosine)
1069 {
1070 #if defined(__GNUC__) && defined(__linux__) && !defined(__clang__)
1071  __builtin_sincos(x, sine, cosine);
1072 #else
1073  *sine = std::sin(x);
1074  *cosine = std::cos(x);
1075 #endif
1076 }
1077 
1078 
1079 
1082 inline float
1083 safe_asinf (float x)
1084 {
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);
1088 }
1089 
1090 
1093 inline float
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);
1098 }
1099 
1100 
1103 inline double
1104 safe_sqrt (double x)
1105 {
1106  return (x > 0.0) ? sqrt(x) : 0.0;
1107 }
1108 
1109 inline float
1110 safe_sqrtf (float x)
1111 {
1112  return (x > 0.0f) ? sqrtf(x) : 0.0f;
1113 }
1114 
1116 inline float safe_inversesqrt (float x) {
1117  return (x > 0.0f) ? 1.0f/std::sqrt(x) : 0.0f;
1118 }
1119 
1120 inline float safe_log (float x) {
1121  return (x > 0.0f) ? logf(x) : -std::numeric_limits<float>::max();
1122 }
1123 
1124 inline float safe_log2(float x) {
1125  return (x > 0.0f) ? log2f(x) : -std::numeric_limits<float>::max();
1126 }
1127 
1128 inline float safe_log10(float x) {
1129  return (x > 0.0f) ? log10f(x) : -std::numeric_limits<float>::max();
1130 }
1131 
1132 inline float safe_logb (float x) {
1133  return (x != 0.0f) ? logbf(x) : -std::numeric_limits<float>::max();
1134 }
1135 
1136 
1137 
1138 
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)
1151 {
1152  // Use the Regula Falsi method, falling back to bisection if it
1153  // hasn't converged after 3/4 of the maximum number of iterations.
1154  // See, e.g., Numerical Recipes for the basic ideas behind both
1155  // methods.
1156  T v0 = func(xmin), v1 = func(xmax);
1157  T x = xmin, v = v0;
1158  bool increasing = (v0 < v1);
1159  T vmin = increasing ? v0 : v1;
1160  T vmax = increasing ? v1 : v0;
1161  bool bracketed = (y >= vmin && y <= vmax);
1162  if (brack)
1163  *brack = bracketed;
1164  if (! bracketed) {
1165  // If our bounds don't bracket the zero, just give up, and
1166  // return the approprate "edge" of the interval
1167  return ((y < vmin) == increasing) ? xmin : xmax;
1168  }
1169  if (fabs(v0-v1) < eps) // already close enough
1170  return x;
1171  int rfiters = (3*maxiters)/4; // how many times to try regula falsi
1172  for (int iters = 0; iters < maxiters; ++iters) {
1173  T t; // interpolation factor
1174  if (iters < rfiters) {
1175  // Regula falsi
1176  t = (y-v0)/(v1-v0);
1177  if (t <= T(0) || t >= T(1))
1178  t = T(0.5); // RF convergence failure -- bisect instead
1179  } else {
1180  t = T(0.5); // bisection
1181  }
1182  x = lerp (xmin, xmax, t);
1183  v = func(x);
1184  if ((v < y) == increasing) {
1185  xmin = x; v0 = v;
1186  } else {
1187  xmax = x; v1 = v;
1188  }
1189  if (fabs(xmax-xmin) < eps || fabs(v-y) < eps)
1190  return x; // converged
1191  }
1192  return x;
1193 }
1194 
1195 
1196 
1197 }
1198 OIIO_NAMESPACE_EXIT
1199 
1200 #endif // OPENIMAGEIO_FMATH_H