|
| 1 | +// Copyright 2013 The Flutter Authors. All rights reserved. |
| 2 | +// Use of this source code is governed by a BSD-style license that can be |
| 3 | +// found in the LICENSE file. |
| 4 | + |
| 5 | +part of engine; |
| 6 | + |
| 7 | +// This is a small library to handle stability for floating point operations. |
| 8 | +// |
| 9 | +// Since we are representing an infinite number of real numbers in finite |
| 10 | +// number of bits, when we perform comparisons of coordinates for paths for |
| 11 | +// example, we want to make sure that line and curve sections that are too |
| 12 | +// close to each other (number of floating point numbers |
| 13 | +// representable in bits between two numbers) are handled correctly and |
| 14 | +// don't cause algorithms to fail when we perform operations such as |
| 15 | +// subtraction or between checks. |
| 16 | +// |
| 17 | +// Small introduction into floating point comparison: |
| 18 | +// |
| 19 | +// For some good articles on the topic, see |
| 20 | +// https://randomascii.wordpress.com/category/floating-point/page/2/ |
| 21 | +// Port based on: |
| 22 | +// https://github.com/google/skia/blob/master/include/private/SkFloatBits.h |
| 23 | +// |
| 24 | +// Here is the 32 bit IEEE representation: |
| 25 | +// uint32_t mantissa : 23; |
| 26 | +// uint32_t exponent : 8; |
| 27 | +// uint32_t sign : 1; |
| 28 | +// As you can see it was carefully designed to be reinterpreted as an integer. |
| 29 | +// |
| 30 | +// Ulps stands for unit in the last place. ulp(x) is the gap between two |
| 31 | +// floating point numbers nearest x. |
| 32 | + |
| 33 | +/// Converts a sign-bit int (float interpreted as int) into a 2s complement |
| 34 | +/// int. Also converts 0x80000000 to 0. Allows result to be compared using |
| 35 | +/// int comparison. |
| 36 | +int signBitTo2sCompliment(int x) => |
| 37 | + (x & 0x80000000) != 0 ? (-(x & 0x7fffffff)) : x; |
| 38 | + |
| 39 | +/// Convert a 2s complement int to a sign-bit (i.e. int interpreted as float). |
| 40 | +int twosComplimentToSignBit(int x) { |
| 41 | + if ((x & 0x80000000) == 0) { |
| 42 | + return x; |
| 43 | + } |
| 44 | + x = ~x + 1; |
| 45 | + x |= 0x80000000; |
| 46 | + return x; |
| 47 | +} |
| 48 | + |
| 49 | +class _FloatBitConverter { |
| 50 | + final Float32List float32List; |
| 51 | + final Int32List int32List; |
| 52 | + _FloatBitConverter._(this.float32List, this.int32List); |
| 53 | + |
| 54 | + factory _FloatBitConverter() { |
| 55 | + Float32List float32List = Float32List(1); |
| 56 | + return _FloatBitConverter._( |
| 57 | + float32List, float32List.buffer.asInt32List(0, 1)); |
| 58 | + } |
| 59 | + |
| 60 | + int toInt(Float32List source, int index) { |
| 61 | + float32List[0] = source[index]; |
| 62 | + return int32List[0]; |
| 63 | + } |
| 64 | + |
| 65 | + int toBits(double x) { |
| 66 | + float32List[0] = x; |
| 67 | + return int32List[0]; |
| 68 | + } |
| 69 | + |
| 70 | + double toDouble(int bits) { |
| 71 | + int32List[0] = bits; |
| 72 | + return float32List[0]; |
| 73 | + } |
| 74 | +} |
| 75 | + |
| 76 | +// Singleton bit converter to prevent typed array allocations. |
| 77 | +final _FloatBitConverter _floatBitConverter = _FloatBitConverter(); |
| 78 | + |
| 79 | +// Converts float to bits. |
| 80 | +int float2Bits(Float32List source, int index) { |
| 81 | + return _floatBitConverter.toInt(source, index); |
| 82 | +} |
| 83 | + |
| 84 | +// Converts bits to float. |
| 85 | +double bitsToFloat(int bits) { |
| 86 | + return _floatBitConverter.toDouble(bits); |
| 87 | +} |
| 88 | + |
| 89 | +const int floatBitsExponentMask = 0x7F800000; |
| 90 | +const int floatBitsMatissaMask = 0x007FFFFF; |
| 91 | + |
| 92 | +/// Returns a float as 2s complement int to be able to compare floats to each |
| 93 | +/// other. |
| 94 | +int floatFromListAs2sCompliment(Float32List source, int index) => |
| 95 | + signBitTo2sCompliment(float2Bits(source, index)); |
| 96 | + |
| 97 | +int floatAs2sCompliment(double x) => |
| 98 | + signBitTo2sCompliment(_floatBitConverter.toBits(x)); |
| 99 | + |
| 100 | +double twosComplimentAsFloat(int x) => bitsToFloat(twosComplimentToSignBit(x)); |
| 101 | + |
| 102 | +bool _argumentsDenormalized(double a, double b, int epsilon) { |
| 103 | + double denormalizedCheck = kFltEpsilon * epsilon / 2; |
| 104 | + return a.abs() <= denormalizedCheck && b.abs() <= denormalizedCheck; |
| 105 | +} |
| 106 | + |
| 107 | +bool equalUlps(double a, double b, int epsilon, int depsilon) { |
| 108 | + if (_argumentsDenormalized(a, b, depsilon)) { |
| 109 | + return true; |
| 110 | + } |
| 111 | + int aBits = floatAs2sCompliment(a); |
| 112 | + int bBits = floatAs2sCompliment(b); |
| 113 | + // Find the difference in ULPs. |
| 114 | + return aBits < bBits + epsilon && bBits < aBits + epsilon; |
| 115 | +} |
| 116 | + |
| 117 | +/// General equality check that covers between, product and division by using |
| 118 | +/// ulps epsilon 16. |
| 119 | +bool almostEqualUlps(double a, double b) { |
| 120 | + const int kUlpsEpsilon = 16; |
| 121 | + return equalUlps(a, b, kUlpsEpsilon, kUlpsEpsilon); |
| 122 | +} |
| 123 | + |
| 124 | +/// Equality using the same error term for between comparison. |
| 125 | +bool almostBequalUlps(double a, double b) { |
| 126 | + const int kUlpsEpsilon = 2; |
| 127 | + return equalUlps(a, b, kUlpsEpsilon, kUlpsEpsilon); |
| 128 | +} |
| 129 | + |
| 130 | +/// Equality check for product. |
| 131 | +bool almostPequalUlps(double a, double b) { |
| 132 | + const int kUlpsEpsilon = 8; |
| 133 | + return equalUlps(a, b, kUlpsEpsilon, kUlpsEpsilon); |
| 134 | +} |
| 135 | + |
| 136 | +/// Equality check for division. |
| 137 | +bool almostDequalUlps(double a, double b) { |
| 138 | + const int kUlpsEpsilon = 16; |
| 139 | + return equalUlps(a, b, kUlpsEpsilon, kUlpsEpsilon); |
| 140 | +} |
| 141 | + |
| 142 | +/// Checks if 2 points are roughly equal (ulp 256) to each other. |
| 143 | +bool approximatelyEqual(double ax, double ay, double bx, double by) { |
| 144 | + if (approximatelyEqualT(ax, bx) && approximatelyEqualT(ay, by)) { |
| 145 | + return true; |
| 146 | + } |
| 147 | + if (!roughlyEqualUlps(ax, bx) || !roughlyEqualUlps(ay, by)) { |
| 148 | + return false; |
| 149 | + } |
| 150 | + final double dx = (ax - bx); |
| 151 | + final double dy = (ay - by); |
| 152 | + double dist = math.sqrt(dx * dx + dy * dy); |
| 153 | + double tiniest = math.min(math.min(math.min(ax, bx), ay), by); |
| 154 | + double largest = math.max(math.max(math.max(ax, bx), ay), by); |
| 155 | + largest = math.max(largest, -tiniest); |
| 156 | + return almostDequalUlps(largest, largest + dist); |
| 157 | +} |
| 158 | + |
| 159 | +/// Equality check for comparing curve T values in the range of 0 to 1. |
| 160 | +/// |
| 161 | +/// For general numbers (larger and smaller) use |
| 162 | +/// AlmostEqualUlps instead. |
| 163 | +bool approximatelyEqualT(double t1, double t2) { |
| 164 | + return approximatelyZero(t1 - t2); |
| 165 | +} |
| 166 | + |
| 167 | +bool approximatelyZero(double value) => value.abs() < kFltEpsilon; |
| 168 | + |
| 169 | +bool roughlyEqualUlps(double a, double b) { |
| 170 | + const int kUlpsEpsilon = 256; |
| 171 | + const int kDUlpsEpsilon = 1024; |
| 172 | + return equalUlps(a, b, kUlpsEpsilon, kDUlpsEpsilon); |
| 173 | +} |
| 174 | + |
| 175 | +bool dEqualUlpsEpsilon(double a, double b, int epsilon) { |
| 176 | + int aBits = floatAs2sCompliment(a); |
| 177 | + int bBits = floatAs2sCompliment(b); |
| 178 | + // Find the difference in ULPs. |
| 179 | + return aBits < bBits + epsilon && bBits < aBits + epsilon; |
| 180 | +} |
| 181 | + |
| 182 | +// Checks equality for division. |
| 183 | +bool almostDequalUlpsDouble(double a, double b) { |
| 184 | + final double absA = a.abs(); |
| 185 | + final double absB = b.abs(); |
| 186 | + if (absA < kScalarMax && absB < kScalarMax) { |
| 187 | + return almostDequalUlps(a, b); |
| 188 | + } |
| 189 | + return (a - b).abs() / math.max(absA, absB) < kDblEpsilonSubdivideErr; |
| 190 | +} |
| 191 | + |
| 192 | +const double kFltEpsilon = 1.19209290E-07; // == 1 / (2 ^ 23) |
| 193 | +const double kDblEpsilon = 2.22045e-16; |
| 194 | +const double kFltEpsilonCubed = kFltEpsilon * kFltEpsilon * kFltEpsilon; |
| 195 | +const double kFltEpsilonHalf = kFltEpsilon / 2; |
| 196 | +const double kFltEpsilonDouble = kFltEpsilon * 2; |
| 197 | +// Epsilon to use when ordering vectors. |
| 198 | +const double kFltEpsilonOrderableErr = kFltEpsilon * 16; |
| 199 | +const double kFltEpsilonSquared = kFltEpsilon * kFltEpsilon; |
| 200 | +// Use a compile-time constant for FLT_EPSILON_SQRT to avoid initializers. |
| 201 | +// A 17 digit constant guarantees exact results. |
| 202 | +const double kFltEpsilonSqrt = 0.00034526697709225118; // sqrt(kFltEpsilon); |
| 203 | +const double kFltEpsilonInverse = 1 / kFltEpsilon; |
| 204 | +const double kDblEpsilonErr = kDblEpsilon * 4; |
| 205 | +const double kDblEpsilonSubdivideErr = kDblEpsilon * 16; |
| 206 | +const double kRoughEpsilon = kFltEpsilon * 64; |
| 207 | +const double kMoreRoughEpsilon = kFltEpsilon * 256; |
| 208 | +const double kWayRoughEpsilon = kFltEpsilon * 2048; |
| 209 | +const double kBumpEpsilon = kFltEpsilon * 4096; |
| 210 | + |
| 211 | +// Scalar max is based on 32 bit float since [PathRef] stores values in |
| 212 | +// Float32List. |
| 213 | +const double kScalarMax = 3.402823466e+38; |
| 214 | +const double kScalarMin = -kScalarMax; |
0 commit comments