// ulps.h: Counts the FP intervals between two FP numbers

1. // ulps.h: Counts the FP intervals between two FP numbers

2. #include <algorithm>

3. #include <cassert>

4. #include <cstdlib>

5. #include <cmath>

6. #include <limits>

7. template<typename FType>

8. int sign(const FType& x) {

9. return x < FType(0) ? -1 : x > FType(0) ? 1 : 0;

10. }

11. template<typename FType>

12. int lowerexp(FType z, FType& glb) {

13. int base = std::numeric_limits<FType>::radix;

14. // Find largest exponent of base <= z ("glb" = "base"^"exp").

15. // z is bracketed by [glb, lub)

16. // First, search above 1

17. FType lub = 1;

18. int exp = 0;

19. while (lub <= z) {

20. lub *= base;

21. ++exp;

22. }

23. glb = lub/base;

24. --exp;

25. // In case z < 1, search there

26. while (glb > z) {

27. lub = glb;

28. glb /= base;

29. --exp;

30. }

31. assert(glb <= z && z < lub);

32. return exp;

33. }

34. template<typename FType>

35. long double ulps(FType x, FType y) {

36. // Handle special cases first

37. long double inf = std::numeric_limits<long double>::infinity();

38. if (x == y)

39. return 0;

40. else if (sign(x)*sign(y) <= 0 || x == inf || y == inf)

41. return inf;

42. // Ignore sign for the result; make x <= y

43. x = std::abs(x);

44. y = std::abs(y);

45. if (x > y)

46. std::swap(x, y);

47. // Get needed floating-point parameters

48. // (The cast is to force long double arithmetic below)

49. long double eps =

50. static_cast<long double>(std::numeric_limits<FType>::epsilon());