code_junkie

Most effective way for float and double comparison

Nov 10th, 2011
135
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
text 11.61 KB | None | 0 0
  1. bool CompareDoubles1 (double A, double B)
  2. {
  3. return A == B;
  4. }
  5.  
  6. bool CompareDoubles2 (double A, double B)
  7. {
  8. diff = A - B;
  9. return (diff < EPSILON) && (-diff > EPSILON);
  10. }
  11.  
  12. bool AreSame(double a, double b)
  13. {
  14. return fabs(a - b) < EPSILON;
  15. }
  16.  
  17. // Usable AlmostEqual function
  18. bool AlmostEqual2sComplement(float A, float B, int maxUlps)
  19. {
  20. // Make sure maxUlps is non-negative and small enough that the
  21. // default NAN won't compare as equal to anything.
  22. assert(maxUlps > 0 && maxUlps < 4 * 1024 * 1024);
  23. int aInt = *(int*)&A;
  24. // Make aInt lexicographically ordered as a twos-complement int
  25. if (aInt < 0)
  26. aInt = 0x80000000 - aInt;
  27. // Make bInt lexicographically ordered as a twos-complement int
  28. int bInt = *(int*)&B;
  29. if (bInt < 0)
  30. bInt = 0x80000000 - bInt;
  31. int intDiff = abs(aInt - bInt);
  32. if (intDiff <= maxUlps)
  33. return true;
  34. return false;
  35. }
  36.  
  37. #include <limits>
  38. std::numeric_limits<double>::epsilon()
  39.  
  40. #include <cmath>
  41. #include <limits>
  42.  
  43. bool AreSame(double a, double b) {
  44. return std::fabs(a - b) < std::numeric_limits<double>::epsilon();
  45. }
  46.  
  47. double left = // something
  48. double right = // something
  49. const FloatingPoint<double> lhs(left), rhs(right);
  50.  
  51. if (lhs.AlmostEquals(rhs)) {
  52. //they're equal!
  53. }
  54.  
  55. // Copyright 2005, Google Inc.
  56. // All rights reserved.
  57. //
  58. // Redistribution and use in source and binary forms, with or without
  59. // modification, are permitted provided that the following conditions are
  60. // met:
  61. //
  62. // * Redistributions of source code must retain the above copyright
  63. // notice, this list of conditions and the following disclaimer.
  64. // * Redistributions in binary form must reproduce the above
  65. // copyright notice, this list of conditions and the following disclaimer
  66. // in the documentation and/or other materials provided with the
  67. // distribution.
  68. // * Neither the name of Google Inc. nor the names of its
  69. // contributors may be used to endorse or promote products derived from
  70. // this software without specific prior written permission.
  71. //
  72. // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
  73. // "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
  74. // LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
  75. // A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
  76. // OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
  77. // SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
  78. // LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
  79. // DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
  80. // THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
  81. // (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
  82. // OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
  83. //
  84. // Authors: [email protected] (Zhanyong Wan), [email protected] (Sean Mcafee)
  85. //
  86. // The Google C++ Testing Framework (Google Test)
  87.  
  88.  
  89. // This template class serves as a compile-time function from size to
  90. // type. It maps a size in bytes to a primitive type with that
  91. // size. e.g.
  92. //
  93. // TypeWithSize<4>::UInt
  94. //
  95. // is typedef-ed to be unsigned int (unsigned integer made up of 4
  96. // bytes).
  97. //
  98. // Such functionality should belong to STL, but I cannot find it
  99. // there.
  100. //
  101. // Google Test uses this class in the implementation of floating-point
  102. // comparison.
  103. //
  104. // For now it only handles UInt (unsigned int) as that's all Google Test
  105. // needs. Other types can be easily added in the future if need
  106. // arises.
  107. template <size_t size>
  108. class TypeWithSize {
  109. public:
  110. // This prevents the user from using TypeWithSize<N> with incorrect
  111. // values of N.
  112. typedef void UInt;
  113. };
  114.  
  115. // The specialization for size 4.
  116. template <>
  117. class TypeWithSize<4> {
  118. public:
  119. // unsigned int has size 4 in both gcc and MSVC.
  120. //
  121. // As base/basictypes.h doesn't compile on Windows, we cannot use
  122. // uint32, uint64, and etc here.
  123. typedef int Int;
  124. typedef unsigned int UInt;
  125. };
  126.  
  127. // The specialization for size 8.
  128. template <>
  129. class TypeWithSize<8> {
  130. public:
  131. #if GTEST_OS_WINDOWS
  132. typedef __int64 Int;
  133. typedef unsigned __int64 UInt;
  134. #else
  135. typedef long long Int; // NOLINT
  136. typedef unsigned long long UInt; // NOLINT
  137. #endif // GTEST_OS_WINDOWS
  138. };
  139.  
  140.  
  141. // This template class represents an IEEE floating-point number
  142. // (either single-precision or double-precision, depending on the
  143. // template parameters).
  144. //
  145. // The purpose of this class is to do more sophisticated number
  146. // comparison. (Due to round-off error, etc, it's very unlikely that
  147. // two floating-points will be equal exactly. Hence a naive
  148. // comparison by the == operation often doesn't work.)
  149. //
  150. // Format of IEEE floating-point:
  151. //
  152. // The most-significant bit being the leftmost, an IEEE
  153. // floating-point looks like
  154. //
  155. // sign_bit exponent_bits fraction_bits
  156. //
  157. // Here, sign_bit is a single bit that designates the sign of the
  158. // number.
  159. //
  160. // For float, there are 8 exponent bits and 23 fraction bits.
  161. //
  162. // For double, there are 11 exponent bits and 52 fraction bits.
  163. //
  164. // More details can be found at
  165. // http://en.wikipedia.org/wiki/IEEE_floating-point_standard.
  166. //
  167. // Template parameter:
  168. //
  169. // RawType: the raw floating-point type (either float or double)
  170. template <typename RawType>
  171. class FloatingPoint {
  172. public:
  173. // Defines the unsigned integer type that has the same size as the
  174. // floating point number.
  175. typedef typename TypeWithSize<sizeof(RawType)>::UInt Bits;
  176.  
  177. // Constants.
  178.  
  179. // # of bits in a number.
  180. static const size_t kBitCount = 8*sizeof(RawType);
  181.  
  182. // # of fraction bits in a number.
  183. static const size_t kFractionBitCount =
  184. std::numeric_limits<RawType>::digits - 1;
  185.  
  186. // # of exponent bits in a number.
  187. static const size_t kExponentBitCount = kBitCount - 1 - kFractionBitCount;
  188.  
  189. // The mask for the sign bit.
  190. static const Bits kSignBitMask = static_cast<Bits>(1) << (kBitCount - 1);
  191.  
  192. // The mask for the fraction bits.
  193. static const Bits kFractionBitMask =
  194. ~static_cast<Bits>(0) >> (kExponentBitCount + 1);
  195.  
  196. // The mask for the exponent bits.
  197. static const Bits kExponentBitMask = ~(kSignBitMask | kFractionBitMask);
  198.  
  199. // How many ULP's (Units in the Last Place) we want to tolerate when
  200. // comparing two numbers. The larger the value, the more error we
  201. // allow. A 0 value means that two numbers must be exactly the same
  202. // to be considered equal.
  203. //
  204. // The maximum error of a single floating-point operation is 0.5
  205. // units in the last place. On Intel CPU's, all floating-point
  206. // calculations are done with 80-bit precision, while double has 64
  207. // bits. Therefore, 4 should be enough for ordinary use.
  208. //
  209. // See the following article for more details on ULP:
  210. // http://www.cygnus-software.com/papers/comparingfloats/comparingfloats.htm.
  211. static const size_t kMaxUlps = 4;
  212.  
  213. // Constructs a FloatingPoint from a raw floating-point number.
  214. //
  215. // On an Intel CPU, passing a non-normalized NAN (Not a Number)
  216. // around may change its bits, although the new value is guaranteed
  217. // to be also a NAN. Therefore, don't expect this constructor to
  218. // preserve the bits in x when x is a NAN.
  219. explicit FloatingPoint(const RawType& x) { u_.value_ = x; }
  220.  
  221. // Static methods
  222.  
  223. // Reinterprets a bit pattern as a floating-point number.
  224. //
  225. // This function is needed to test the AlmostEquals() method.
  226. static RawType ReinterpretBits(const Bits bits) {
  227. FloatingPoint fp(0);
  228. fp.u_.bits_ = bits;
  229. return fp.u_.value_;
  230. }
  231.  
  232. // Returns the floating-point number that represent positive infinity.
  233. static RawType Infinity() {
  234. return ReinterpretBits(kExponentBitMask);
  235. }
  236.  
  237. // Non-static methods
  238.  
  239. // Returns the bits that represents this number.
  240. const Bits &bits() const { return u_.bits_; }
  241.  
  242. // Returns the exponent bits of this number.
  243. Bits exponent_bits() const { return kExponentBitMask & u_.bits_; }
  244.  
  245. // Returns the fraction bits of this number.
  246. Bits fraction_bits() const { return kFractionBitMask & u_.bits_; }
  247.  
  248. // Returns the sign bit of this number.
  249. Bits sign_bit() const { return kSignBitMask & u_.bits_; }
  250.  
  251. // Returns true iff this is NAN (not a number).
  252. bool is_nan() const {
  253. // It's a NAN if the exponent bits are all ones and the fraction
  254. // bits are not entirely zeros.
  255. return (exponent_bits() == kExponentBitMask) && (fraction_bits() != 0);
  256. }
  257.  
  258. // Returns true iff this number is at most kMaxUlps ULP's away from
  259. // rhs. In particular, this function:
  260. //
  261. // - returns false if either number is (or both are) NAN.
  262. // - treats really large numbers as almost equal to infinity.
  263. // - thinks +0.0 and -0.0 are 0 DLP's apart.
  264. bool AlmostEquals(const FloatingPoint& rhs) const {
  265. // The IEEE standard says that any comparison operation involving
  266. // a NAN must return false.
  267. if (is_nan() || rhs.is_nan()) return false;
  268.  
  269. return DistanceBetweenSignAndMagnitudeNumbers(u_.bits_, rhs.u_.bits_)
  270. <= kMaxUlps;
  271. }
  272.  
  273. private:
  274. // The data type used to store the actual floating-point number.
  275. union FloatingPointUnion {
  276. RawType value_; // The raw floating-point number.
  277. Bits bits_; // The bits that represent the number.
  278. };
  279.  
  280. // Converts an integer from the sign-and-magnitude representation to
  281. // the biased representation. More precisely, let N be 2 to the
  282. // power of (kBitCount - 1), an integer x is represented by the
  283. // unsigned number x + N.
  284. //
  285. // For instance,
  286. //
  287. // -N + 1 (the most negative number representable using
  288. // sign-and-magnitude) is represented by 1;
  289. // 0 is represented by N; and
  290. // N - 1 (the biggest number representable using
  291. // sign-and-magnitude) is represented by 2N - 1.
  292. //
  293. // Read http://en.wikipedia.org/wiki/Signed_number_representations
  294. // for more details on signed number representations.
  295. static Bits SignAndMagnitudeToBiased(const Bits &sam) {
  296. if (kSignBitMask & sam) {
  297. // sam represents a negative number.
  298. return ~sam + 1;
  299. } else {
  300. // sam represents a positive number.
  301. return kSignBitMask | sam;
  302. }
  303. }
  304.  
  305. // Given two numbers in the sign-and-magnitude representation,
  306. // returns the distance between them as an unsigned number.
  307. static Bits DistanceBetweenSignAndMagnitudeNumbers(const Bits &sam1,
  308. const Bits &sam2) {
  309. const Bits biased1 = SignAndMagnitudeToBiased(sam1);
  310. const Bits biased2 = SignAndMagnitudeToBiased(sam2);
  311. return (biased1 >= biased2) ? (biased1 - biased2) : (biased2 - biased1);
  312. }
  313.  
  314. FloatingPointUnion u_;
  315. };
  316.  
  317. return (diff < EPSILON) && (-diff > EPSILON);
  318.  
  319. return (diff < EPSILON) && (diff > -EPSILON);
  320.  
  321. return a==b;
  322.  
  323. #define VERYSMALL (1.0E-150)
  324. #define EPSILON (1.0E-8)
  325. bool AreSame(double a, double b)
  326. {
  327. double absDiff = fabs(a - b);
  328. if (absDiff < VERYSMALL)
  329. {
  330. return true;
  331. }
  332.  
  333. double maxAbs = max(fabs(a) - fabs(b));
  334. return (absDiff/maxAbs) < EPSILON;
  335. }
  336.  
  337. bool approximatelyEqual(float a, float b, float epsilon)
  338. {
  339. return fabs(a - b) <= ( (fabs(a) < fabs(b) ? fabs(b) : fabs(a)) * epsilon);
  340. }
  341.  
  342. bool essentiallyEqual(float a, float b, float epsilon)
  343. {
  344. return fabs(a - b) <= ( (fabs(a) > fabs(b) ? fabs(b) : fabs(a)) * epsilon);
  345. }
  346.  
  347. bool definitelyGreaterThan(float a, float b, float epsilon)
  348. {
  349. return (a - b) > ( (fabs(a) < fabs(b) ? fabs(b) : fabs(a)) * epsilon);
  350. }
  351.  
  352. bool definitelyLessThan(float a, float b, float epsilon)
  353. {
  354. return (b - a) > ( (fabs(a) < fabs(b) ? fabs(b) : fabs(a)) * epsilon);
  355. }
  356.  
  357. if (d1 == d2) return true;
  358.  
  359. if (fabs(a - b) <= DBL_EPSILON * fmax(fabs(a), fabs(b)))
  360. {
  361. // ...
  362. }
  363.  
  364. if (fabs(a - b) <= 16 * DBL_EPSILON * fmax(fabs(a), fabs(b)))
  365. {
  366. // ...
  367. }
Add Comment
Please, Sign In to add comment