/src/gdal/port/cpl_float.cpp
Line | Count | Source (jump to first uncovered line) |
1 | | /****************************************************************************** |
2 | | * |
3 | | * Project: CPL |
4 | | * Purpose: Floating point conversion functions. Convert 16- and 24-bit |
5 | | * floating point numbers into the 32-bit IEEE 754 compliant ones. |
6 | | * Author: Andrey Kiselev, dron@remotesensing.org |
7 | | * |
8 | | ****************************************************************************** |
9 | | * Copyright (c) 2005, Andrey Kiselev <dron@remotesensing.org> |
10 | | * |
11 | | * This code is based on the code from OpenEXR project with the following |
12 | | * copyright: |
13 | | * |
14 | | * Copyright (c) 2002, Industrial Light & Magic, a division of Lucas |
15 | | * Digital Ltd. LLC |
16 | | * |
17 | | * All rights reserved. |
18 | | * |
19 | | * Redistribution and use in source and binary forms, with or without |
20 | | * modification, are permitted provided that the following conditions are |
21 | | * met: |
22 | | * * Redistributions of source code must retain the above copyright |
23 | | * notice, this list of conditions and the following disclaimer. |
24 | | * * Redistributions in binary form must reproduce the above |
25 | | * copyright notice, this list of conditions and the following disclaimer |
26 | | * in the documentation and/or other materials provided with the |
27 | | * distribution. |
28 | | * * Neither the name of Industrial Light & Magic nor the names of |
29 | | * its contributors may be used to endorse or promote products derived |
30 | | * from this software without specific prior written permission. |
31 | | * |
32 | | * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS |
33 | | * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT |
34 | | * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR |
35 | | * A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT |
36 | | * OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, |
37 | | * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT |
38 | | * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, |
39 | | * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY |
40 | | * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT |
41 | | * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE |
42 | | * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. |
43 | | * |
44 | | ****************************************************************************/ |
45 | | |
46 | | #include "cpl_float.h" |
47 | | #include "cpl_error.h" |
48 | | |
49 | | #include <algorithm> |
50 | | #include <cmath> |
51 | | #include <cstring> |
52 | | #include <limits> |
53 | | #include <numeric> |
54 | | #include <optional> |
55 | | |
56 | | /************************************************************************/ |
57 | | /* HalfToFloat() */ |
58 | | /* */ |
59 | | /* 16-bit floating point number to 32-bit one. */ |
60 | | /************************************************************************/ |
61 | | |
62 | | GUInt32 CPLHalfToFloat(GUInt16 iHalf) |
63 | 0 | { |
64 | |
|
65 | 0 | GUInt32 iSign = (iHalf >> 15) & 0x00000001; |
66 | 0 | int iExponent = (iHalf >> 10) & 0x0000001f; |
67 | 0 | GUInt32 iMantissa = iHalf & 0x000003ff; |
68 | |
|
69 | 0 | if (iExponent == 0) |
70 | 0 | { |
71 | 0 | if (iMantissa == 0) |
72 | 0 | { |
73 | | /* -------------------------------------------------------------------- |
74 | | */ |
75 | | /* Plus or minus zero. */ |
76 | | /* -------------------------------------------------------------------- |
77 | | */ |
78 | |
|
79 | 0 | return iSign << 31; |
80 | 0 | } |
81 | 0 | else |
82 | 0 | { |
83 | | /* -------------------------------------------------------------------- |
84 | | */ |
85 | | /* Denormalized number -- renormalize it. */ |
86 | | /* -------------------------------------------------------------------- |
87 | | */ |
88 | |
|
89 | 0 | while (!(iMantissa & 0x00000400)) |
90 | 0 | { |
91 | 0 | iMantissa <<= 1; |
92 | 0 | iExponent -= 1; |
93 | 0 | } |
94 | |
|
95 | 0 | iExponent += 1; |
96 | 0 | iMantissa &= ~0x00000400U; |
97 | 0 | } |
98 | 0 | } |
99 | 0 | else if (iExponent == 31) |
100 | 0 | { |
101 | 0 | if (iMantissa == 0) |
102 | 0 | { |
103 | | /* -------------------------------------------------------------------- |
104 | | */ |
105 | | /* Positive or negative infinity. */ |
106 | | /* -------------------------------------------------------------------- |
107 | | */ |
108 | |
|
109 | 0 | return (iSign << 31) | 0x7f800000; |
110 | 0 | } |
111 | 0 | else |
112 | 0 | { |
113 | | /* -------------------------------------------------------------------- |
114 | | */ |
115 | | /* NaN -- preserve sign and significand bits. */ |
116 | | /* -------------------------------------------------------------------- |
117 | | */ |
118 | |
|
119 | 0 | return (iSign << 31) | 0x7f800000 | (iMantissa << 13); |
120 | 0 | } |
121 | 0 | } |
122 | | |
123 | | /* -------------------------------------------------------------------- */ |
124 | | /* Normalized number. */ |
125 | | /* -------------------------------------------------------------------- */ |
126 | | |
127 | 0 | iExponent = iExponent + (127 - 15); |
128 | 0 | iMantissa = iMantissa << 13; |
129 | | |
130 | | /* -------------------------------------------------------------------- */ |
131 | | /* Assemble sign, exponent and mantissa. */ |
132 | | /* -------------------------------------------------------------------- */ |
133 | | |
134 | | /* coverity[overflow_sink] */ |
135 | 0 | return (iSign << 31) | (static_cast<GUInt32>(iExponent) << 23) | iMantissa; |
136 | 0 | } |
137 | | |
138 | | /************************************************************************/ |
139 | | /* TripleToFloat() */ |
140 | | /* */ |
141 | | /* 24-bit floating point number to 32-bit one. */ |
142 | | /************************************************************************/ |
143 | | |
144 | | GUInt32 CPLTripleToFloat(GUInt32 iTriple) |
145 | 150k | { |
146 | | |
147 | 150k | GUInt32 iSign = (iTriple >> 23) & 0x00000001; |
148 | 150k | int iExponent = (iTriple >> 16) & 0x0000007f; |
149 | 150k | GUInt32 iMantissa = iTriple & 0x0000ffff; |
150 | | |
151 | 150k | if (iExponent == 0) |
152 | 22.6k | { |
153 | 22.6k | if (iMantissa == 0) |
154 | 7.35k | { |
155 | | /* -------------------------------------------------------------------- |
156 | | */ |
157 | | /* Plus or minus zero. */ |
158 | | /* -------------------------------------------------------------------- |
159 | | */ |
160 | | |
161 | 7.35k | return iSign << 31; |
162 | 7.35k | } |
163 | 15.3k | else |
164 | 15.3k | { |
165 | | /* -------------------------------------------------------------------- |
166 | | */ |
167 | | /* Denormalized number -- renormalize it. */ |
168 | | /* -------------------------------------------------------------------- |
169 | | */ |
170 | | |
171 | 159k | while (!(iMantissa & 0x00010000)) |
172 | 144k | { |
173 | 144k | iMantissa <<= 1; |
174 | 144k | iExponent -= 1; |
175 | 144k | } |
176 | | |
177 | 15.3k | iExponent += 1; |
178 | 15.3k | iMantissa &= ~0x00010000U; |
179 | 15.3k | } |
180 | 22.6k | } |
181 | 128k | else if (iExponent == 127) |
182 | 55.2k | { |
183 | 55.2k | if (iMantissa == 0) |
184 | 2.14k | { |
185 | | /* -------------------------------------------------------------------- |
186 | | */ |
187 | | /* Positive or negative infinity. */ |
188 | | /* -------------------------------------------------------------------- |
189 | | */ |
190 | | |
191 | 2.14k | return (iSign << 31) | 0x7f800000; |
192 | 2.14k | } |
193 | 53.1k | else |
194 | 53.1k | { |
195 | | /* -------------------------------------------------------------------- |
196 | | */ |
197 | | /* NaN -- preserve sign and significand bits. */ |
198 | | /* -------------------------------------------------------------------- |
199 | | */ |
200 | | |
201 | 53.1k | return (iSign << 31) | 0x7f800000 | (iMantissa << 7); |
202 | 53.1k | } |
203 | 55.2k | } |
204 | | |
205 | | /* -------------------------------------------------------------------- */ |
206 | | /* Normalized number. */ |
207 | | /* -------------------------------------------------------------------- */ |
208 | | |
209 | 88.0k | iExponent = iExponent + (127 - 63); |
210 | 88.0k | iMantissa = iMantissa << 7; |
211 | | |
212 | | /* -------------------------------------------------------------------- */ |
213 | | /* Assemble sign, exponent and mantissa. */ |
214 | | /* -------------------------------------------------------------------- */ |
215 | | |
216 | | /* coverity[overflow_sink] */ |
217 | 88.0k | return (iSign << 31) | (static_cast<GUInt32>(iExponent) << 23) | iMantissa; |
218 | 150k | } |
219 | | |
220 | | /************************************************************************/ |
221 | | /* FloatToHalf() */ |
222 | | /************************************************************************/ |
223 | | |
224 | | GUInt16 CPLFloatToHalf(GUInt32 iFloat32, bool &bHasWarned) |
225 | 0 | { |
226 | 0 | GUInt32 iSign = (iFloat32 >> 31) & 0x00000001; |
227 | 0 | GUInt32 iExponent = (iFloat32 >> 23) & 0x000000ff; |
228 | 0 | GUInt32 iMantissa = iFloat32 & 0x007fffff; |
229 | |
|
230 | 0 | if (iExponent == 255) |
231 | 0 | { |
232 | 0 | if (iMantissa == 0) |
233 | 0 | { |
234 | | /* -------------------------------------------------------------------- |
235 | | */ |
236 | | /* Positive or negative infinity. */ |
237 | | /* -------------------------------------------------------------------- |
238 | | */ |
239 | |
|
240 | 0 | return static_cast<GUInt16>((iSign << 15) | 0x7C00); |
241 | 0 | } |
242 | 0 | else |
243 | 0 | { |
244 | | /* -------------------------------------------------------------------- |
245 | | */ |
246 | | /* NaN -- preserve sign and significand bits. */ |
247 | | /* -------------------------------------------------------------------- |
248 | | */ |
249 | 0 | if (iMantissa >> 13) |
250 | 0 | return static_cast<GUInt16>((iSign << 15) | 0x7C00 | |
251 | 0 | (iMantissa >> 13)); |
252 | | |
253 | 0 | return static_cast<GUInt16>((iSign << 15) | 0x7E00); |
254 | 0 | } |
255 | 0 | } |
256 | | |
257 | 0 | if (iExponent <= 127 - 15) |
258 | 0 | { |
259 | | // Zero, float32 denormalized number or float32 too small normalized |
260 | | // number |
261 | 0 | if (13 + 1 + 127 - 15 - iExponent >= 32) |
262 | 0 | return static_cast<GUInt16>(iSign << 15); |
263 | | |
264 | | // Return a denormalized number |
265 | 0 | return static_cast<GUInt16>( |
266 | 0 | (iSign << 15) | |
267 | 0 | ((iMantissa | 0x00800000) >> (13 + 1 + 127 - 15 - iExponent))); |
268 | 0 | } |
269 | 0 | if (iExponent - (127 - 15) >= 31) |
270 | 0 | { |
271 | 0 | if (!bHasWarned) |
272 | 0 | { |
273 | 0 | bHasWarned = true; |
274 | 0 | float fVal = 0.0f; |
275 | 0 | memcpy(&fVal, &iFloat32, 4); |
276 | 0 | CPLError( |
277 | 0 | CE_Failure, CPLE_AppDefined, |
278 | 0 | "Value %.8g is beyond range of float16. Converted to %sinf", |
279 | 0 | fVal, (fVal > 0) ? "+" : "-"); |
280 | 0 | } |
281 | 0 | return static_cast<GUInt16>((iSign << 15) | 0x7C00); // Infinity |
282 | 0 | } |
283 | | |
284 | | /* -------------------------------------------------------------------- */ |
285 | | /* Normalized number. */ |
286 | | /* -------------------------------------------------------------------- */ |
287 | | |
288 | 0 | iExponent = iExponent - (127 - 15); |
289 | 0 | iMantissa = iMantissa >> 13; |
290 | | |
291 | | /* -------------------------------------------------------------------- */ |
292 | | /* Assemble sign, exponent and mantissa. */ |
293 | | /* -------------------------------------------------------------------- */ |
294 | | |
295 | | // coverity[overflow_sink] |
296 | 0 | return static_cast<GUInt16>((iSign << 15) | (iExponent << 10) | iMantissa); |
297 | 0 | } |
298 | | |
299 | | GUInt16 CPLConvertFloatToHalf(float fFloat32) |
300 | 0 | { |
301 | 0 | GUInt32 nFloat32; |
302 | 0 | std::memcpy(&nFloat32, &fFloat32, sizeof nFloat32); |
303 | 0 | bool bHasWarned = true; |
304 | 0 | return CPLFloatToHalf(nFloat32, bHasWarned); |
305 | 0 | } |
306 | | |
307 | | float CPLConvertHalfToFloat(GUInt16 nHalf) |
308 | 0 | { |
309 | 0 | GUInt32 nFloat32 = CPLHalfToFloat(nHalf); |
310 | 0 | float fFloat32; |
311 | 0 | std::memcpy(&fFloat32, &nFloat32, sizeof fFloat32); |
312 | 0 | return fFloat32; |
313 | 0 | } |
314 | | |
315 | | namespace |
316 | | { |
317 | | |
318 | | template <typename T> struct Fraction |
319 | | { |
320 | | using value_type = T; |
321 | | |
322 | | T num; |
323 | | T denom; |
324 | | }; |
325 | | |
326 | | /** Approximate a floating point number as a fraction, using the method describe |
327 | | * in Richards, Ian (1981). Continued Fractions Without Tears. Mathematics |
328 | | * Magazine, Vol. 54, No. 4. https://doi.org/10.2307/2689627 |
329 | | * |
330 | | * If the fraction cannot be approximated within the specified error tolerance |
331 | | * in a certain amount of iterations, a warning will be raised and std::nullopt |
332 | | * will be returned. |
333 | | * |
334 | | * @param x the number to approximate as a fraction |
335 | | * @param err the maximum allowable absolute error in the approximation |
336 | | * |
337 | | * @return the approximated value, or std::nullopt |
338 | | * |
339 | | */ |
340 | | std::optional<Fraction<std::uint64_t>> FloatToFraction(double x, double err) |
341 | 0 | { |
342 | 0 | using inttype = std::uint64_t; |
343 | 0 | constexpr int MAX_ITER = 1000; |
344 | |
|
345 | 0 | const double sign = std::signbit(x) ? -1 : 1; |
346 | |
|
347 | 0 | double g(std::abs(x)); |
348 | 0 | inttype a(0); |
349 | 0 | inttype b(1); |
350 | 0 | inttype c(1); |
351 | 0 | inttype d(0); |
352 | |
|
353 | 0 | Fraction<std::uint64_t> ret; |
354 | |
|
355 | 0 | for (int i = 0; i < MAX_ITER; i++) |
356 | 0 | { |
357 | 0 | if (!(g >= 0 && |
358 | 0 | g <= static_cast<double>(std::numeric_limits<inttype>::max()))) |
359 | 0 | { |
360 | 0 | break; |
361 | 0 | } |
362 | 0 | const inttype s = static_cast<inttype>(std::floor(g)); |
363 | 0 | ret.num = a + s * c; |
364 | 0 | ret.denom = b + s * d; |
365 | |
|
366 | 0 | a = c; |
367 | 0 | b = d; |
368 | 0 | c = ret.num; |
369 | 0 | d = ret.denom; |
370 | 0 | g = 1.0 / (g - static_cast<double>(s)); |
371 | |
|
372 | 0 | const double approx = sign * static_cast<double>(ret.num) / |
373 | 0 | static_cast<double>(ret.denom); |
374 | |
|
375 | 0 | if (std::abs(approx - x) < err) |
376 | 0 | { |
377 | 0 | return ret; |
378 | 0 | } |
379 | 0 | } |
380 | | |
381 | 0 | CPLError(CE_Warning, CPLE_AppDefined, |
382 | 0 | "Failed to approximate %g as a fraction with error < %g in %d " |
383 | 0 | "iterations", |
384 | 0 | x, err, MAX_ITER); |
385 | 0 | return std::nullopt; |
386 | 0 | } |
387 | | } // namespace |
388 | | |
389 | | /** Return the largest value by which two input values can be |
390 | | * divided, with the result being an integer. If no suitable |
391 | | * value can be found, zero will be returned. |
392 | | */ |
393 | | double CPLGreatestCommonDivisor(double a, double b) |
394 | 0 | { |
395 | 0 | if (a == 0 || !std::isfinite(a) || b == 0 || !std::isfinite(b)) |
396 | 0 | { |
397 | 0 | CPLError(CE_Failure, CPLE_AppDefined, |
398 | 0 | "Input values must be finite non-null values"); |
399 | 0 | return 0; |
400 | 0 | } |
401 | | |
402 | 0 | if (a == b) |
403 | 0 | { |
404 | 0 | return a; |
405 | 0 | } |
406 | | |
407 | | // Check if one resolution is an integer factor of the other. |
408 | | // This is fast and succeeds in some cases where the method below fails. |
409 | 0 | if (a > b && std::abs(std::round(a / b) - a / b) < 1e-8) |
410 | 0 | { |
411 | 0 | return b; |
412 | 0 | } |
413 | 0 | if (b > a && std::abs(std::round(b / a) - b / a) < 1e-8) |
414 | 0 | { |
415 | 0 | return a; |
416 | 0 | } |
417 | | |
418 | 0 | const auto approx_a = FloatToFraction(a, 1e-10); |
419 | 0 | if (!approx_a.has_value()) |
420 | 0 | { |
421 | 0 | CPLError(CE_Failure, CPLE_AppDefined, |
422 | 0 | "Could not approximate resolution %.18g as a fraction", a); |
423 | 0 | return 0; |
424 | 0 | } |
425 | | |
426 | 0 | const auto approx_b = FloatToFraction(b, 1e-10); |
427 | 0 | if (!approx_b.has_value()) |
428 | 0 | { |
429 | 0 | CPLError(CE_Failure, CPLE_AppDefined, |
430 | 0 | "Could not approximate resolution %.18g as a fraction", b); |
431 | 0 | return 0; |
432 | 0 | } |
433 | | |
434 | 0 | const double sign = std::signbit(a) ? -1 : 1; |
435 | |
|
436 | 0 | const auto &frac_a = approx_a.value(); |
437 | 0 | const auto &frac_b = approx_b.value(); |
438 | |
|
439 | 0 | const auto common_denom = std::lcm(frac_a.denom, frac_b.denom); |
440 | |
|
441 | 0 | const auto num_a = static_cast<std::uint64_t>( |
442 | 0 | frac_a.num * std::round(common_denom / frac_a.denom)); |
443 | 0 | const auto num_b = static_cast<std::uint64_t>( |
444 | 0 | frac_b.num * std::round(common_denom / frac_b.denom)); |
445 | |
|
446 | 0 | const auto common_num = std::gcd(num_a, num_b); |
447 | | |
448 | | // Add std::numeric_limits<double>::min() to avoid Coverity Scan warning |
449 | | // about div by zero |
450 | 0 | const auto common = sign * static_cast<double>(common_num) / |
451 | 0 | (static_cast<double>(common_denom) + |
452 | 0 | std::numeric_limits<double>::min()); |
453 | |
|
454 | 0 | const auto disaggregation_factor = std::max(a / common, b / common); |
455 | 0 | if (disaggregation_factor > 10000) |
456 | 0 | { |
457 | 0 | CPLError(CE_Failure, CPLE_AppDefined, |
458 | 0 | "Common resolution between %.18g and %.18g calculated at " |
459 | 0 | "%.18g which " |
460 | 0 | "would cause excessive disaggregation", |
461 | 0 | a, b, common); |
462 | 0 | return 0; |
463 | 0 | } |
464 | | |
465 | 0 | return common; |
466 | 0 | } |