1/// @ref gtc_ulp
2/// @file glm/gtc/ulp.inl
3///
4/// Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
5///
6/// Developed at SunPro, a Sun Microsystems, Inc. business.
7/// Permission to use, copy, modify, and distribute this
8/// software is freely granted, provided that this notice
9/// is preserved.
10
11#include "../detail/type_int.hpp"
12#include <cmath>
13#include <cfloat>
14#include <limits>
15
16#if(GLM_COMPILER & GLM_COMPILER_VC)
17#	pragma warning(push)
18#	pragma warning(disable : 4127)
19#endif
20
21typedef union
22{
23	float value;
24	/* FIXME: Assumes 32 bit int.  */
25	unsigned int word;
26} ieee_float_shape_type;
27
28typedef union
29{
30	double value;
31	struct
32	{
33		glm::detail::int32 lsw;
34		glm::detail::int32 msw;
35	} parts;
36} ieee_double_shape_type;
37
38#define GLM_EXTRACT_WORDS(ix0,ix1,d)		\
39	do {									\
40		ieee_double_shape_type ew_u;		\
41		ew_u.value = (d);					\
42		(ix0) = ew_u.parts.msw;				\
43		(ix1) = ew_u.parts.lsw;				\
44	} while (0)
45
46#define GLM_GET_FLOAT_WORD(i,d)				\
47	do {									\
48		ieee_float_shape_type gf_u;			\
49		gf_u.value = (d);					\
50		(i) = gf_u.word;					\
51	} while (0)
52
53#define GLM_SET_FLOAT_WORD(d,i)				\
54	do {									\
55		ieee_float_shape_type sf_u;			\
56		sf_u.word = (i);					\
57		(d) = sf_u.value;					\
58	} while (0)
59
60#define GLM_INSERT_WORDS(d,ix0,ix1)			\
61	do {									\
62		ieee_double_shape_type iw_u;		\
63		iw_u.parts.msw = (ix0);				\
64		iw_u.parts.lsw = (ix1);				\
65		(d) = iw_u.value;					\
66	} while (0)
67
68namespace glm{
69namespace detail
70{
71	GLM_FUNC_QUALIFIER float nextafterf(float x, float y)
72	{
73		volatile float t;
74		glm::detail::int32 hx, hy, ix, iy;
75
76		GLM_GET_FLOAT_WORD(hx, x);
77		GLM_GET_FLOAT_WORD(hy, y);
78		ix = hx&0x7fffffff;		// |x|
79		iy = hy&0x7fffffff;		// |y|
80
81		if((ix>0x7f800000) ||	// x is nan
82			(iy>0x7f800000))	// y is nan
83			return x+y;
84		if(x==y) return y;		// x=y, return y
85		if(ix==0) {				// x == 0
86			GLM_SET_FLOAT_WORD(x,(hy&0x80000000)|1);// return +-minsubnormal
87			t = x*x;
88			if(t==x) return t; else return x;	// raise underflow flag
89		}
90		if(hx>=0) {				// x > 0
91			if(hx>hy) {			// x > y, x -= ulp
92				hx -= 1;
93			} else {			// x < y, x += ulp
94				hx += 1;
95			}
96		} else {				// x < 0
97			if(hy>=0||hx>hy){	// x < y, x -= ulp
98				hx -= 1;
99			} else {			// x > y, x += ulp
100				hx += 1;
101			}
102		}
103		hy = hx&0x7f800000;
104		if(hy>=0x7f800000) return x+x;  // overflow
105		if(hy<0x00800000) {             // underflow
106			t = x*x;
107			if(t!=x) {          // raise underflow flag
108				GLM_SET_FLOAT_WORD(y,hx);
109				return y;
110			}
111		}
112		GLM_SET_FLOAT_WORD(x,hx);
113		return x;
114	}
115
116	GLM_FUNC_QUALIFIER double nextafter(double x, double y)
117	{
118		volatile double t;
119		glm::detail::int32 hx, hy, ix, iy;
120		glm::detail::uint32 lx, ly;
121
122		GLM_EXTRACT_WORDS(hx, lx, x);
123		GLM_EXTRACT_WORDS(hy, ly, y);
124		ix = hx & 0x7fffffff;             // |x|
125		iy = hy & 0x7fffffff;             // |y|
126
127		if(((ix>=0x7ff00000)&&((ix-0x7ff00000)|lx)!=0) ||   // x is nan
128			((iy>=0x7ff00000)&&((iy-0x7ff00000)|ly)!=0))     // y is nan
129			return x+y;
130		if(x==y) return y;              // x=y, return y
131		if((ix|lx)==0) {                        // x == 0
132			GLM_INSERT_WORDS(x, hy & 0x80000000, 1);    // return +-minsubnormal
133			t = x*x;
134			if(t==x) return t; else return x;   // raise underflow flag
135		}
136		if(hx>=0) {                             // x > 0
137			if(hx>hy||((hx==hy)&&(lx>ly))) {    // x > y, x -= ulp
138				if(lx==0) hx -= 1;
139				lx -= 1;
140			} else {                            // x < y, x += ulp
141				lx += 1;
142				if(lx==0) hx += 1;
143			}
144		} else {                                // x < 0
145			if(hy>=0||hx>hy||((hx==hy)&&(lx>ly))){// x < y, x -= ulp
146				if(lx==0) hx -= 1;
147				lx -= 1;
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) return x+x;  // overflow
155		if(hy<0x00100000) {             // underflow
156			t = x*x;
157			if(t!=x) {          // raise underflow flag
158				GLM_INSERT_WORDS(y,hx,lx);
159				return y;
160			}
161		}
162		GLM_INSERT_WORDS(x,hx,lx);
163		return x;
164	}
165}//namespace detail
166}//namespace glm
167
168#if(GLM_COMPILER & GLM_COMPILER_VC)
169#	pragma warning(pop)
170#endif
171
172namespace glm
173{
174	template <>
175	GLM_FUNC_QUALIFIER float next_float(float const & x)
176	{
177#		if GLM_HAS_CXX11_STL
178			return std::nextafter(x, std::numeric_limits<float>::max());
179#		elif((GLM_COMPILER & GLM_COMPILER_VC) || ((GLM_COMPILER & GLM_COMPILER_INTEL) && (GLM_PLATFORM & GLM_PLATFORM_WINDOWS)))
180			return detail::nextafterf(x, FLT_MAX);
181#		elif(GLM_PLATFORM & GLM_PLATFORM_ANDROID)
182			return __builtin_nextafterf(x, FLT_MAX);
183#		else
184			return nextafterf(x, FLT_MAX);
185#		endif
186	}
187
188	template <>
189	GLM_FUNC_QUALIFIER double next_float(double const & 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, FLT_MAX);
197#		else
198			return nextafter(x, DBL_MAX);
199#		endif
200	}
201
202	template<typename T, precision P, template<typename, precision> class vecType>
203	GLM_FUNC_QUALIFIER vecType<T, P> next_float(vecType<T, P> const & x)
204	{
205		vecType<T, P> Result(uninitialize);
206		for(length_t i = 0, n = Result.length(); i < n; ++i)
207			Result[i] = next_float(x[i]);
208		return Result;
209	}
210
211	GLM_FUNC_QUALIFIER float prev_float(float const & 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
224	GLM_FUNC_QUALIFIER double prev_float(double const & x)
225	{
226#		if GLM_HAS_CXX11_STL
227			return std::nextafter(x, std::numeric_limits<double>::min());
228#		elif((GLM_COMPILER & GLM_COMPILER_VC) || ((GLM_COMPILER & GLM_COMPILER_INTEL) && (GLM_PLATFORM & GLM_PLATFORM_WINDOWS)))
229			return _nextafter(x, DBL_MIN);
230#		elif(GLM_PLATFORM & GLM_PLATFORM_ANDROID)
231			return __builtin_nextafter(x, DBL_MIN);
232#		else
233			return nextafter(x, DBL_MIN);
234#		endif
235	}
236
237	template<typename T, precision P, template<typename, precision> class vecType>
238	GLM_FUNC_QUALIFIER vecType<T, P> prev_float(vecType<T, P> const & x)
239	{
240		vecType<T, P> Result(uninitialize);
241		for(length_t i = 0, n = Result.length(); i < n; ++i)
242			Result[i] = prev_float(x[i]);
243		return Result;
244	}
245
246	template <typename T>
247	GLM_FUNC_QUALIFIER T next_float(T const & x, uint const & ulps)
248	{
249		T temp = x;
250		for(uint i = 0; i < ulps; ++i)
251			temp = next_float(temp);
252		return temp;
253	}
254
255	template<typename T, precision P, template<typename, precision> class vecType>
256	GLM_FUNC_QUALIFIER vecType<T, P> next_float(vecType<T, P> const & x, vecType<uint, P> const & ulps)
257	{
258		vecType<T, P> Result(uninitialize);
259		for(length_t i = 0, n = Result.length(); i < n; ++i)
260			Result[i] = next_float(x[i], ulps[i]);
261		return Result;
262	}
263
264	template <typename T>
265	GLM_FUNC_QUALIFIER T prev_float(T const & x, uint const & ulps)
266	{
267		T temp = x;
268		for(uint i = 0; i < ulps; ++i)
269			temp = prev_float(temp);
270		return temp;
271	}
272
273	template<typename T, precision P, template<typename, precision> class vecType>
274	GLM_FUNC_QUALIFIER vecType<T, P> prev_float(vecType<T, P> const & x, vecType<uint, P> const & ulps)
275	{
276		vecType<T, P> Result(uninitialize);
277		for(length_t i = 0, n = Result.length(); i < n; ++i)
278			Result[i] = prev_float(x[i], ulps[i]);
279		return Result;
280	}
281
282	template <typename T>
283	GLM_FUNC_QUALIFIER uint float_distance(T const & x, T const & y)
284	{
285		uint ulp = 0;
286
287		if(x < y)
288		{
289			T temp = x;
290			while(temp != y)// && ulp < std::numeric_limits<std::size_t>::max())
291			{
292				++ulp;
293				temp = next_float(temp);
294			}
295		}
296		else if(y < x)
297		{
298			T temp = y;
299			while(temp != x)// && ulp < std::numeric_limits<std::size_t>::max())
300			{
301				++ulp;
302				temp = next_float(temp);
303			}
304		}
305		else // ==
306		{
307
308		}
309
310		return ulp;
311	}
312
313	template<typename T, precision P, template<typename, precision> class vecType>
314	GLM_FUNC_QUALIFIER vecType<uint, P> float_distance(vecType<T, P> const & x, vecType<T, P> const & y)
315	{
316		vecType<uint, P> Result(uninitialize);
317		for(length_t i = 0, n = Result.length(); i < n; ++i)
318			Result[i] = float_distance(x[i], y[i]);
319		return Result;
320	}
321}//namespace glm
322