/glm/gtx/ulp.inl
C++ Header | 299 lines | 248 code | 32 blank | 19 comment | 54 complexity | cf25eb10af997fd49dae0cc839970cbf MD5 | raw file
1///////////////////////////////////////////////////////////////////////////////////////////////////
2// OpenGL Mathematics Copyright (c) 2005 - 2012 G-Truc Creation (www.g-truc.net)
3///////////////////////////////////////////////////////////////////////////////////////////////////
4// Created : 2011-03-07
5// Updated : 2011-04-26
6// Licence : This source is under MIT License
7// File : glm/gtx/ulp.inl
8///////////////////////////////////////////////////////////////////////////////////////////////////
9
10#include <cmath>
11#include <cfloat>
12
13/*
14 * ====================================================
15 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
16 *
17 * Developed at SunPro, a Sun Microsystems, Inc. business.
18 * Permission to use, copy, modify, and distribute this
19 * software is freely granted, provided that this notice
20 * is preserved.
21 * ====================================================
22 */
23
24#pragma warning(push)
25#pragma warning(disable : 4127)
26
27typedef union
28{
29 float value;
30 /* FIXME: Assumes 32 bit int. */
31 unsigned int word;
32} ieee_float_shape_type;
33
34typedef union
35{
36 double value;
37 struct
38 {
39 glm::detail::int32 lsw;
40 glm::detail::int32 msw;
41 } parts;
42} ieee_double_shape_type;
43
44#define GLM_EXTRACT_WORDS(ix0,ix1,d) \
45 do { \
46 ieee_double_shape_type ew_u; \
47 ew_u.value = (d); \
48 (ix0) = ew_u.parts.msw; \
49 (ix1) = ew_u.parts.lsw; \
50 } while (0)
51
52#define GLM_GET_FLOAT_WORD(i,d) \
53 do { \
54 ieee_float_shape_type gf_u; \
55 gf_u.value = (d); \
56 (i) = gf_u.word; \
57 } while (0)
58
59#define GLM_SET_FLOAT_WORD(d,i) \
60 do { \
61 ieee_float_shape_type sf_u; \
62 sf_u.word = (i); \
63 (d) = sf_u.value; \
64 } while (0)
65
66#define GLM_INSERT_WORDS(d,ix0,ix1) \
67 do { \
68 ieee_double_shape_type iw_u; \
69 iw_u.parts.msw = (ix0); \
70 iw_u.parts.lsw = (ix1); \
71 (d) = iw_u.value; \
72 } while (0)
73
74namespace glm{
75namespace detail
76{
77 GLM_FUNC_QUALIFIER float nextafterf(float x, float y)
78 {
79 volatile float t;
80 glm::detail::int32 hx, hy, ix, iy;
81
82 GLM_GET_FLOAT_WORD(hx, x);
83 GLM_GET_FLOAT_WORD(hy, y);
84 ix = hx&0x7fffffff; // |x|
85 iy = hy&0x7fffffff; // |y|
86
87 if((ix>0x7f800000) || // x is nan
88 (iy>0x7f800000)) // y is nan
89 return x+y;
90 if(x==y) return y; // x=y, return y
91 if(ix==0) { // x == 0
92 GLM_SET_FLOAT_WORD(x,(hy&0x80000000)|1);// return +-minsubnormal
93 t = x*x;
94 if(t==x) return t; else return x; // raise underflow flag
95 }
96 if(hx>=0) { // x > 0
97 if(hx>hy) { // x > y, x -= ulp
98 hx -= 1;
99 } else { // x < y, x += ulp
100 hx += 1;
101 }
102 } else { // x < 0
103 if(hy>=0||hx>hy){ // x < y, x -= ulp
104 hx -= 1;
105 } else { // x > y, x += ulp
106 hx += 1;
107 }
108 }
109 hy = hx&0x7f800000;
110 if(hy>=0x7f800000) return x+x; // overflow
111 if(hy<0x00800000) { // underflow
112 t = x*x;
113 if(t!=x) { // raise underflow flag
114 GLM_SET_FLOAT_WORD(y,hx);
115 return y;
116 }
117 }
118 GLM_SET_FLOAT_WORD(x,hx);
119 return x;
120 }
121
122 GLM_FUNC_QUALIFIER double nextafter(double x, double y)
123 {
124 volatile double t;
125 glm::detail::int32 hx, hy, ix, iy;
126 glm::detail::uint32 lx, ly;
127
128 GLM_EXTRACT_WORDS(hx, lx, x);
129 GLM_EXTRACT_WORDS(hy, ly, y);
130 ix = hx & 0x7fffffff; // |x|
131 iy = hy & 0x7fffffff; // |y|
132
133 if(((ix>=0x7ff00000)&&((ix-0x7ff00000)|lx)!=0) || // x is nan
134 ((iy>=0x7ff00000)&&((iy-0x7ff00000)|ly)!=0)) // y is nan
135 return x+y;
136 if(x==y) return y; // x=y, return y
137 if((ix|lx)==0) { // x == 0
138 GLM_INSERT_WORDS(x, hy & 0x80000000, 1); // return +-minsubnormal
139 t = x*x;
140 if(t==x) return t; else return x; // raise underflow flag
141 }
142 if(hx>=0) { // x > 0
143 if(hx>hy||((hx==hy)&&(lx>ly))) { // x > y, x -= ulp
144 if(lx==0) hx -= 1;
145 lx -= 1;
146 } else { // x < y, x += ulp
147 lx += 1;
148 if(lx==0) hx += 1;
149 }
150 } else { // x < 0
151 if(hy>=0||hx>hy||((hx==hy)&&(lx>ly))){// x < y, x -= ulp
152 if(lx==0) hx -= 1;
153 lx -= 1;
154 } else { // x > y, x += ulp
155 lx += 1;
156 if(lx==0) hx += 1;
157 }
158 }
159 hy = hx&0x7ff00000;
160 if(hy>=0x7ff00000) return x+x; // overflow
161 if(hy<0x00100000) { // underflow
162 t = x*x;
163 if(t!=x) { // raise underflow flag
164 GLM_INSERT_WORDS(y,hx,lx);
165 return y;
166 }
167 }
168 GLM_INSERT_WORDS(x,hx,lx);
169 return x;
170 }
171}//namespace detail
172}//namespace glm
173
174#pragma warning(pop)
175
176#if(GLM_COMPILER & GLM_COMPILER_VC || GLM_COMPILER & GLM_COMPILER_INTEL)
177# define GLM_NEXT_AFTER_FLT(x, toward) glm::detail::nextafterf((x), (toward))
178# define GLM_NEXT_AFTER_DBL(x, toward) _nextafter((x), (toward))
179#else
180# define GLM_NEXT_AFTER_FLT(x, toward) nextafterf((x), (toward))
181# define GLM_NEXT_AFTER_DBL(x, toward) nextafter((x), (toward))
182#endif
183
184namespace glm
185{
186 GLM_FUNC_QUALIFIER float next_float(float const & x)
187 {
188 return GLM_NEXT_AFTER_FLT(x, std::numeric_limits<float>::max());
189 }
190
191 GLM_FUNC_QUALIFIER double next_float(double const & x)
192 {
193 return GLM_NEXT_AFTER_DBL(x, std::numeric_limits<double>::max());
194 }
195
196 template<typename T, template<typename> class vecType>
197 GLM_FUNC_QUALIFIER vecType<T> next_float(vecType<T> const & x)
198 {
199 vecType<T> Result;
200 for(std::size_t i = 0; i < Result.length(); ++i)
201 Result[i] = next_float(x[i]);
202 return Result;
203 }
204
205 GLM_FUNC_QUALIFIER float prev_float(float const & x)
206 {
207 return GLM_NEXT_AFTER_FLT(x, std::numeric_limits<float>::min());
208 }
209
210 GLM_FUNC_QUALIFIER double prev_float(double const & x)
211 {
212 return GLM_NEXT_AFTER_DBL(x, std::numeric_limits<double>::min());
213 }
214
215 template<typename T, template<typename> class vecType>
216 GLM_FUNC_QUALIFIER vecType<T> prev_float(vecType<T> const & x)
217 {
218 vecType<T> Result;
219 for(std::size_t i = 0; i < Result.length(); ++i)
220 Result[i] = prev_float(x[i]);
221 return Result;
222 }
223
224 template <typename T>
225 GLM_FUNC_QUALIFIER T next_float(T const & x, uint const & ulps)
226 {
227 T temp = x;
228 for(std::size_t i = 0; i < ulps; ++i)
229 temp = next_float(temp);
230 return temp;
231 }
232
233 template<typename T, template<typename> class vecType>
234 GLM_FUNC_QUALIFIER vecType<T> next_float(vecType<T> const & x, vecType<uint> const & ulps)
235 {
236 vecType<T> Result;
237 for(std::size_t i = 0; i < Result.length(); ++i)
238 Result[i] = next_float(x[i], ulps[i]);
239 return Result;
240 }
241
242 template <typename T>
243 GLM_FUNC_QUALIFIER T prev_float(T const & x, uint const & ulps)
244 {
245 T temp = x;
246 for(std::size_t i = 0; i < ulps; ++i)
247 temp = prev_float(temp);
248 return temp;
249 }
250
251 template<typename T, template<typename> class vecType>
252 GLM_FUNC_QUALIFIER vecType<T> prev_float(vecType<T> const & x, vecType<uint> const & ulps)
253 {
254 vecType<T> Result;
255 for(std::size_t i = 0; i < Result.length(); ++i)
256 Result[i] = prev_float(x[i], ulps[i]);
257 return Result;
258 }
259
260 template <typename T>
261 GLM_FUNC_QUALIFIER uint float_distance(T const & x, T const & y)
262 {
263 uint ulp = 0;
264
265 if(x < y)
266 {
267 T temp = x;
268 while(temp != y && ulp < std::numeric_limits<std::size_t>::max())
269 {
270 ++ulp;
271 temp = next_float(temp);
272 }
273 }
274 else if(y < x)
275 {
276 T temp = y;
277 while(temp != x && ulp < std::numeric_limits<std::size_t>::max())
278 {
279 ++ulp;
280 temp = next_float(temp);
281 }
282 }
283 else // ==
284 {
285
286 }
287
288 return ulp;
289 }
290
291 template<typename T, template<typename> class vecType>
292 GLM_FUNC_QUALIFIER vecType<uint> float_distance(vecType<T> const & x, vecType<T> const & y)
293 {
294 vecType<uint> Result;
295 for(std::size_t i = 0; i < Result.length(); ++i)
296 Result[i] = float_distance(x[i], y[i]);
297 return Result;
298 }
299}//namespace glm