scalar_ulp.inl 6.9 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284
  1. /// Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
  2. ///
  3. /// Developed at SunPro, a Sun Microsystems, Inc. business.
  4. /// Permission to use, copy, modify, and distribute this
  5. /// software is freely granted, provided that this notice
  6. /// is preserved.
  7. #include "../detail/type_float.hpp"
  8. #include "../ext/scalar_constants.hpp"
  9. #include <cmath>
  10. #include <cfloat>
  11. #if(GLM_COMPILER & GLM_COMPILER_VC)
  12. # pragma warning(push)
  13. # pragma warning(disable : 4127)
  14. #endif
  15. typedef union
  16. {
  17. float value;
  18. /* FIXME: Assumes 32 bit int. */
  19. unsigned int word;
  20. } ieee_float_shape_type;
  21. typedef union
  22. {
  23. double value;
  24. struct
  25. {
  26. int lsw;
  27. int msw;
  28. } parts;
  29. } ieee_double_shape_type;
  30. #define GLM_EXTRACT_WORDS(ix0,ix1,d) \
  31. do { \
  32. ieee_double_shape_type ew_u; \
  33. ew_u.value = (d); \
  34. (ix0) = ew_u.parts.msw; \
  35. (ix1) = ew_u.parts.lsw; \
  36. } while (0)
  37. #define GLM_GET_FLOAT_WORD(i,d) \
  38. do { \
  39. ieee_float_shape_type gf_u; \
  40. gf_u.value = (d); \
  41. (i) = gf_u.word; \
  42. } while (0)
  43. #define GLM_SET_FLOAT_WORD(d,i) \
  44. do { \
  45. ieee_float_shape_type sf_u; \
  46. sf_u.word = (i); \
  47. (d) = sf_u.value; \
  48. } while (0)
  49. #define GLM_INSERT_WORDS(d,ix0,ix1) \
  50. do { \
  51. ieee_double_shape_type iw_u; \
  52. iw_u.parts.msw = (ix0); \
  53. iw_u.parts.lsw = (ix1); \
  54. (d) = iw_u.value; \
  55. } while (0)
  56. namespace glm{
  57. namespace detail
  58. {
  59. GLM_FUNC_QUALIFIER float nextafterf(float x, float y)
  60. {
  61. volatile float t;
  62. int hx, hy, ix, iy;
  63. GLM_GET_FLOAT_WORD(hx, x);
  64. GLM_GET_FLOAT_WORD(hy, y);
  65. ix = hx & 0x7fffffff; // |x|
  66. iy = hy & 0x7fffffff; // |y|
  67. if((ix > 0x7f800000) || // x is nan
  68. (iy > 0x7f800000)) // y is nan
  69. return x + y;
  70. if(abs(y - x) <= epsilon<float>())
  71. return y; // x=y, return y
  72. if(ix == 0)
  73. { // x == 0
  74. GLM_SET_FLOAT_WORD(x, (hy & 0x80000000) | 1);// return +-minsubnormal
  75. t = x * x;
  76. if(abs(t - x) <= epsilon<float>())
  77. return t;
  78. else
  79. return x; // raise underflow flag
  80. }
  81. if(hx >= 0)
  82. { // x > 0
  83. if(hx > hy) // x > y, x -= ulp
  84. hx -= 1;
  85. else // x < y, x += ulp
  86. hx += 1;
  87. }
  88. else
  89. { // x < 0
  90. if(hy >= 0 || hx > hy) // x < y, x -= ulp
  91. hx -= 1;
  92. else // x > y, x += ulp
  93. hx += 1;
  94. }
  95. hy = hx & 0x7f800000;
  96. if(hy >= 0x7f800000)
  97. return x + x; // overflow
  98. if(hy < 0x00800000) // underflow
  99. {
  100. t = x * x;
  101. if(abs(t - x) > epsilon<float>())
  102. { // raise underflow flag
  103. GLM_SET_FLOAT_WORD(y, hx);
  104. return y;
  105. }
  106. }
  107. GLM_SET_FLOAT_WORD(x, hx);
  108. return x;
  109. }
  110. GLM_FUNC_QUALIFIER double nextafter(double x, double y)
  111. {
  112. volatile double t;
  113. int hx, hy, ix, iy;
  114. unsigned int lx, ly;
  115. GLM_EXTRACT_WORDS(hx, lx, x);
  116. GLM_EXTRACT_WORDS(hy, ly, y);
  117. ix = hx & 0x7fffffff; // |x|
  118. iy = hy & 0x7fffffff; // |y|
  119. if(((ix >= 0x7ff00000) && ((ix - 0x7ff00000) | lx) != 0) || // x is nan
  120. ((iy >= 0x7ff00000) && ((iy - 0x7ff00000) | ly) != 0)) // y is nan
  121. return x + y;
  122. if(abs(y - x) <= epsilon<double>())
  123. return y; // x=y, return y
  124. if((ix | lx) == 0)
  125. { // x == 0
  126. GLM_INSERT_WORDS(x, hy & 0x80000000, 1); // return +-minsubnormal
  127. t = x * x;
  128. if(abs(t - x) <= epsilon<double>())
  129. return t;
  130. else
  131. return x; // raise underflow flag
  132. }
  133. if(hx >= 0) { // x > 0
  134. if(hx > hy || ((hx == hy) && (lx > ly))) { // x > y, x -= ulp
  135. if(lx == 0) hx -= 1;
  136. lx -= 1;
  137. }
  138. else { // x < y, x += ulp
  139. lx += 1;
  140. if(lx == 0) hx += 1;
  141. }
  142. }
  143. else { // x < 0
  144. if(hy >= 0 || hx > hy || ((hx == hy) && (lx > ly))){// x < y, x -= ulp
  145. if(lx == 0) hx -= 1;
  146. lx -= 1;
  147. }
  148. else { // x > y, x += ulp
  149. lx += 1;
  150. if(lx == 0) hx += 1;
  151. }
  152. }
  153. hy = hx & 0x7ff00000;
  154. if(hy >= 0x7ff00000)
  155. return x + x; // overflow
  156. if(hy < 0x00100000)
  157. { // underflow
  158. t = x * x;
  159. if(abs(t - x) > epsilon<double>())
  160. { // raise underflow flag
  161. GLM_INSERT_WORDS(y, hx, lx);
  162. return y;
  163. }
  164. }
  165. GLM_INSERT_WORDS(x, hx, lx);
  166. return x;
  167. }
  168. }//namespace detail
  169. }//namespace glm
  170. #if(GLM_COMPILER & GLM_COMPILER_VC)
  171. # pragma warning(pop)
  172. #endif
  173. namespace glm
  174. {
  175. template<>
  176. GLM_FUNC_QUALIFIER float nextFloat(float x)
  177. {
  178. # if GLM_HAS_CXX11_STL
  179. return std::nextafter(x, std::numeric_limits<float>::max());
  180. # elif((GLM_COMPILER & GLM_COMPILER_VC) || ((GLM_COMPILER & GLM_COMPILER_INTEL) && (GLM_PLATFORM & GLM_PLATFORM_WINDOWS)))
  181. return detail::nextafterf(x, FLT_MAX);
  182. # elif(GLM_PLATFORM & GLM_PLATFORM_ANDROID)
  183. return __builtin_nextafterf(x, FLT_MAX);
  184. # else
  185. return nextafterf(x, FLT_MAX);
  186. # endif
  187. }
  188. template<>
  189. GLM_FUNC_QUALIFIER double nextFloat(double x)
  190. {
  191. # if GLM_HAS_CXX11_STL
  192. return std::nextafter(x, std::numeric_limits<double>::max());
  193. # elif((GLM_COMPILER & GLM_COMPILER_VC) || ((GLM_COMPILER & GLM_COMPILER_INTEL) && (GLM_PLATFORM & GLM_PLATFORM_WINDOWS)))
  194. return detail::nextafter(x, std::numeric_limits<double>::max());
  195. # elif(GLM_PLATFORM & GLM_PLATFORM_ANDROID)
  196. return __builtin_nextafter(x, DBL_MAX);
  197. # else
  198. return nextafter(x, DBL_MAX);
  199. # endif
  200. }
  201. template<typename T>
  202. GLM_FUNC_QUALIFIER T nextFloat(T x, int ULPs)
  203. {
  204. GLM_STATIC_ASSERT(std::numeric_limits<T>::is_iec559, "'next_float' only accept floating-point input");
  205. assert(ULPs >= 0);
  206. T temp = x;
  207. for(int i = 0; i < ULPs; ++i)
  208. temp = nextFloat(temp);
  209. return temp;
  210. }
  211. GLM_FUNC_QUALIFIER float prevFloat(float x)
  212. {
  213. # if GLM_HAS_CXX11_STL
  214. return std::nextafter(x, std::numeric_limits<float>::min());
  215. # elif((GLM_COMPILER & GLM_COMPILER_VC) || ((GLM_COMPILER & GLM_COMPILER_INTEL) && (GLM_PLATFORM & GLM_PLATFORM_WINDOWS)))
  216. return detail::nextafterf(x, FLT_MIN);
  217. # elif(GLM_PLATFORM & GLM_PLATFORM_ANDROID)
  218. return __builtin_nextafterf(x, FLT_MIN);
  219. # else
  220. return nextafterf(x, FLT_MIN);
  221. # endif
  222. }
  223. GLM_FUNC_QUALIFIER double prevFloat(double x)
  224. {
  225. # if GLM_HAS_CXX11_STL
  226. return std::nextafter(x, std::numeric_limits<double>::min());
  227. # elif((GLM_COMPILER & GLM_COMPILER_VC) || ((GLM_COMPILER & GLM_COMPILER_INTEL) && (GLM_PLATFORM & GLM_PLATFORM_WINDOWS)))
  228. return _nextafter(x, DBL_MIN);
  229. # elif(GLM_PLATFORM & GLM_PLATFORM_ANDROID)
  230. return __builtin_nextafter(x, DBL_MIN);
  231. # else
  232. return nextafter(x, DBL_MIN);
  233. # endif
  234. }
  235. template<typename T>
  236. GLM_FUNC_QUALIFIER T prevFloat(T x, int ULPs)
  237. {
  238. GLM_STATIC_ASSERT(std::numeric_limits<T>::is_iec559, "'prev_float' only accept floating-point input");
  239. assert(ULPs >= 0);
  240. T temp = x;
  241. for(int i = 0; i < ULPs; ++i)
  242. temp = prevFloat(temp);
  243. return temp;
  244. }
  245. GLM_FUNC_QUALIFIER int floatDistance(float x, float y)
  246. {
  247. detail::float_t<float> const a(x);
  248. detail::float_t<float> const b(y);
  249. return abs(a.i - b.i);
  250. }
  251. GLM_FUNC_QUALIFIER int64 floatDistance(double x, double y)
  252. {
  253. detail::float_t<double> const a(x);
  254. detail::float_t<double> const b(y);
  255. return abs(a.i - b.i);
  256. }
  257. }//namespace glm