Line data Source code
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 : * Copyright (c) 2010, Even Rouault <even dot rouault at spatialys.com>
11 : *
12 : * This code is based on the code from OpenEXR project with the following
13 : * copyright:
14 : *
15 : * Copyright (c) 2002, Industrial Light & Magic, a division of Lucas
16 : * Digital Ltd. LLC
17 : *
18 : * All rights reserved.
19 : *
20 : * Redistribution and use in source and binary forms, with or without
21 : * modification, are permitted provided that the following conditions are
22 : * met:
23 : * * Redistributions of source code must retain the above copyright
24 : * notice, this list of conditions and the following disclaimer.
25 : * * Redistributions in binary form must reproduce the above
26 : * copyright notice, this list of conditions and the following disclaimer
27 : * in the documentation and/or other materials provided with the
28 : * distribution.
29 : * * Neither the name of Industrial Light & Magic nor the names of
30 : * its contributors may be used to endorse or promote products derived
31 : * from this software without specific prior written permission.
32 : *
33 : * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
34 : * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
35 : * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
36 : * A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
37 : * OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
38 : * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
39 : * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
40 : * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
41 : * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
42 : * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
43 : * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
44 : *
45 : ****************************************************************************/
46 :
47 : #ifndef CPL_FLOAT_H_INCLUDED
48 : #define CPL_FLOAT_H_INCLUDED
49 :
50 : #include "cpl_port.h"
51 :
52 : #ifdef __cplusplus
53 : #include <algorithm>
54 : #include <cmath>
55 : #include <cstdint>
56 : #include <cstring>
57 : #include <limits>
58 : #ifdef HAVE_STD_FLOAT16_T
59 : #include <stdfloat>
60 : #endif
61 : #endif
62 :
63 : CPL_C_START
64 : GUInt32 CPL_DLL CPLHalfToFloat(GUInt16 iHalf);
65 : GUInt32 CPL_DLL CPLTripleToFloat(GUInt32 iTriple);
66 : CPL_C_END
67 :
68 : #ifdef __cplusplus
69 :
70 : GUInt16 CPL_DLL CPLFloatToHalf(GUInt32 iFloat32, bool &bHasWarned);
71 :
72 : GUInt16 CPL_DLL CPLConvertFloatToHalf(float fFloat32);
73 : float CPL_DLL CPLConvertHalfToFloat(GUInt16 nHalf);
74 :
75 : namespace cpl
76 : {
77 :
78 : // We define our own version of `std::numeric_limits` so that we can
79 : // specialize it for `cpl::Float16` if necessary. Specializing
80 : // `std::numeric_limits` doesn't always work because some libraries
81 : // use `std::numeric_limits`, and one cannot specialize a type
82 : // template after it has been used.
83 : template <typename T> struct NumericLimits : std::numeric_limits<T>
84 : {
85 : };
86 :
87 : #ifndef HAVE_STD_FLOAT16_T
88 :
89 : // Define a type `cpl::Float16`. If the compiler supports it natively
90 : // (as `_Float16`), then this class is a simple wrapper. Otherwise we
91 : // store the values in a `GUInt16` as bit pattern.
92 :
93 : //! @cond Doxygen_Suppress
94 : struct Float16
95 : {
96 : struct make_from_bits_and_value
97 : {
98 : };
99 :
100 : #ifdef HAVE__FLOAT16
101 :
102 : // How we represent a `Float16` internally
103 : using repr = _Float16;
104 :
105 : // How we compute on `Float16` values
106 : using compute = _Float16;
107 :
108 : // Create a Float16 in a constexpr manner. Since we can't convert
109 : // bits in a constexpr function, we need to take both the bit
110 : // pattern and a float value as input, and can then choose which
111 : // of the two to use.
112 : constexpr Float16(make_from_bits_and_value, CPL_UNUSED std::uint16_t bits,
113 : float fValue)
114 : : rValue(repr(fValue))
115 : {
116 : }
117 :
118 : static constexpr repr computeToRepr(compute fValue)
119 : {
120 : return fValue;
121 : }
122 :
123 : static constexpr compute reprToCompute(repr rValue)
124 : {
125 : return rValue;
126 : }
127 :
128 : template <typename T> static constexpr repr toRepr(T fValue)
129 : {
130 : return static_cast<repr>(fValue);
131 : }
132 :
133 : template <typename T> static constexpr T fromRepr(repr rValue)
134 : {
135 : return static_cast<T>(rValue);
136 : }
137 :
138 : #else // #ifndef HAVE__FLOAT16
139 :
140 : // How we represent a `Float16` internally
141 : using repr = std::uint16_t;
142 :
143 : // How we compute on `Float16` values
144 : using compute = float;
145 :
146 : // Create a Float16 in a constexpr manner. Since we can't convert
147 : // bits in a constexpr function, we need to take both the bit
148 : // pattern and a float value as input, and can then choose which
149 : // of the two to use.
150 4781 : constexpr Float16(make_from_bits_and_value, std::uint16_t bits,
151 : CPL_UNUSED float fValue)
152 4781 : : rValue(bits)
153 : {
154 4781 : }
155 :
156 2445769 : static unsigned float2unsigned(float f)
157 : {
158 : unsigned u;
159 2445769 : std::memcpy(&u, &f, 4);
160 2445769 : return u;
161 : }
162 :
163 2447791 : static float unsigned2float(unsigned u)
164 : {
165 : float f;
166 2447791 : std::memcpy(&f, &u, 4);
167 2447791 : return f;
168 : }
169 :
170 : // Copied from cpl_float.cpp so that we can inline for performance
171 2445769 : static std::uint16_t computeToRepr(float fFloat32)
172 : {
173 2445769 : std::uint32_t iFloat32 = float2unsigned(fFloat32);
174 :
175 2445769 : std::uint32_t iSign = (iFloat32 >> 31) & 0x00000001;
176 2445769 : std::uint32_t iExponent = (iFloat32 >> 23) & 0x000000ff;
177 2445769 : std::uint32_t iMantissa = iFloat32 & 0x007fffff;
178 :
179 2445769 : if (iExponent == 255)
180 : {
181 58946 : if (iMantissa == 0)
182 : {
183 : // Positive or negative infinity.
184 58941 : return static_cast<std::int16_t>((iSign << 15) | 0x7C00);
185 : }
186 :
187 : // NaN -- preserve sign and significand bits.
188 5 : if (iMantissa >> 13)
189 3 : return static_cast<std::int16_t>((iSign << 15) | 0x7C00 |
190 3 : (iMantissa >> 13));
191 2 : return static_cast<std::int16_t>((iSign << 15) | 0x7E00);
192 : }
193 :
194 2386821 : if (iExponent <= 127 - 15)
195 : {
196 : // Zero, float32 denormalized number or float32 too small normalized
197 : // number
198 114196 : if (13 + 1 + 127 - 15 - iExponent >= 32)
199 110885 : return static_cast<std::int16_t>(iSign << 15);
200 :
201 : // Return a denormalized number
202 3311 : return static_cast<std::int16_t>(
203 3311 : (iSign << 15) |
204 3311 : ((iMantissa | 0x00800000) >> (13 + 1 + 127 - 15 - iExponent)));
205 : }
206 :
207 2272626 : if (iExponent - (127 - 15) >= 31)
208 : {
209 18452 : return static_cast<std::int16_t>((iSign << 15) |
210 18452 : 0x7C00); // Infinity
211 : }
212 :
213 : // Normalized number.
214 2254174 : iExponent = iExponent - (127 - 15);
215 2254174 : iMantissa = iMantissa >> 13;
216 :
217 : // Assemble sign, exponent and mantissa.
218 : // coverity[overflow_sink]
219 2254174 : return static_cast<std::int16_t>((iSign << 15) | (iExponent << 10) |
220 2254174 : iMantissa);
221 : }
222 :
223 : // Copied from cpl_float.cpp so that we can inline for performance
224 2447791 : static float reprToCompute(std::uint16_t iHalf)
225 : {
226 2447791 : std::uint32_t iSign = (iHalf >> 15) & 0x00000001;
227 2447791 : int iExponent = (iHalf >> 10) & 0x0000001f;
228 2447791 : std::uint32_t iMantissa = iHalf & 0x000003ff;
229 :
230 2447791 : if (iExponent == 31)
231 : {
232 96789 : if (iMantissa == 0)
233 : {
234 : // Positive or negative infinity.
235 96784 : return unsigned2float((iSign << 31) | 0x7f800000);
236 : }
237 :
238 : // NaN -- preserve sign and significand bits.
239 5 : return unsigned2float((iSign << 31) | 0x7f800000 |
240 5 : (iMantissa << 13));
241 : }
242 :
243 2351004 : if (iExponent == 0)
244 : {
245 150494 : if (iMantissa == 0)
246 : {
247 : // Plus or minus zero.
248 147764 : return unsigned2float(iSign << 31);
249 : }
250 :
251 : // Denormalized number -- renormalize it.
252 18172 : while (!(iMantissa & 0x00000400))
253 : {
254 15442 : iMantissa <<= 1;
255 15442 : iExponent -= 1;
256 : }
257 :
258 2730 : iExponent += 1;
259 2730 : iMantissa &= ~0x00000400U;
260 : }
261 :
262 : // Normalized number.
263 2203233 : iExponent = iExponent + (127 - 15);
264 2203233 : iMantissa = iMantissa << 13;
265 :
266 : // Assemble sign, exponent and mantissa.
267 : /* coverity[overflow_sink] */
268 2203233 : return unsigned2float((iSign << 31) |
269 2203233 : (static_cast<std::uint32_t>(iExponent) << 23) |
270 2203233 : iMantissa);
271 : }
272 :
273 2445774 : template <typename T> static repr toRepr(T fValue)
274 : {
275 2445774 : return computeToRepr(static_cast<compute>(fValue));
276 : }
277 :
278 914038 : template <typename T> static T fromRepr(repr rValue)
279 : {
280 914038 : return static_cast<T>(reprToCompute(rValue));
281 : }
282 :
283 : #endif // #ifndef HAVE__FLOAT16
284 :
285 : private:
286 : repr rValue;
287 :
288 : public:
289 1533748 : compute get() const
290 : {
291 1533748 : return reprToCompute(rValue);
292 : }
293 :
294 : // cppcheck-suppress uninitMemberVar
295 : Float16() = default;
296 : Float16(const Float16 &) = default;
297 : Float16(Float16 &&) = default;
298 : Float16 &operator=(const Float16 &) = default;
299 : Float16 &operator=(Float16 &&) = default;
300 :
301 : // Constructors and conversion operators
302 :
303 : #ifdef HAVE__FLOAT16
304 : // cppcheck-suppress noExplicitConstructor
305 : constexpr Float16(_Float16 hfValue) : rValue(hfValue)
306 : {
307 : }
308 :
309 : constexpr operator _Float16() const
310 : {
311 : return rValue;
312 : }
313 : #endif
314 :
315 : // cppcheck-suppress-macro noExplicitConstructor
316 : #define GDAL_DEFINE_CONVERSION(TYPE) \
317 : \
318 : Float16(TYPE fValue) : rValue(toRepr(fValue)) \
319 : { \
320 : } \
321 : \
322 : operator TYPE() const \
323 : { \
324 : return fromRepr<TYPE>(rValue); \
325 : }
326 :
327 1088202 : GDAL_DEFINE_CONVERSION(float)
328 1967089 : GDAL_DEFINE_CONVERSION(double)
329 : GDAL_DEFINE_CONVERSION(char)
330 1066 : GDAL_DEFINE_CONVERSION(signed char)
331 8171 : GDAL_DEFINE_CONVERSION(short)
332 25245 : GDAL_DEFINE_CONVERSION(int)
333 11637 : GDAL_DEFINE_CONVERSION(long)
334 4330 : GDAL_DEFINE_CONVERSION(long long)
335 243147 : GDAL_DEFINE_CONVERSION(unsigned char)
336 3271 : GDAL_DEFINE_CONVERSION(unsigned short)
337 3281 : GDAL_DEFINE_CONVERSION(unsigned int)
338 2324 : GDAL_DEFINE_CONVERSION(unsigned long)
339 2049 : GDAL_DEFINE_CONVERSION(unsigned long long)
340 :
341 : #undef GDAL_DEFINE_CONVERSION
342 :
343 : // Arithmetic operators
344 :
345 201 : friend Float16 operator+(Float16 x)
346 : {
347 201 : return +x.get();
348 : }
349 :
350 963 : friend Float16 operator-(Float16 x)
351 : {
352 963 : return -x.get();
353 : }
354 :
355 : #define GDAL_DEFINE_ARITHOP(OP) \
356 : \
357 : friend Float16 operator OP(Float16 x, Float16 y) \
358 : { \
359 : return x.get() OP y.get(); \
360 : } \
361 : \
362 : friend double operator OP(double x, Float16 y) \
363 : { \
364 : return x OP y.get(); \
365 : } \
366 : \
367 : friend float operator OP(float x, Float16 y) \
368 : { \
369 : return x OP y.get(); \
370 : } \
371 : \
372 : friend Float16 operator OP(int x, Float16 y) \
373 : { \
374 : return x OP y.get(); \
375 : } \
376 : \
377 : friend double operator OP(Float16 x, double y) \
378 : { \
379 : return x.get() OP y; \
380 : } \
381 : \
382 : friend float operator OP(Float16 x, float y) \
383 : { \
384 : return x.get() OP y; \
385 : } \
386 : \
387 : friend Float16 operator OP(Float16 x, int y) \
388 : { \
389 : return x.get() OP y; \
390 : }
391 :
392 67896 : GDAL_DEFINE_ARITHOP(+)
393 40818 : GDAL_DEFINE_ARITHOP(-)
394 40401 : GDAL_DEFINE_ARITHOP(*)
395 61257 : GDAL_DEFINE_ARITHOP(/)
396 :
397 : #undef GDAL_DEFINE_ARITHOP
398 :
399 : // Comparison operators
400 :
401 : #define GDAL_DEFINE_COMPARISON(OP) \
402 : \
403 : friend bool operator OP(Float16 x, Float16 y) \
404 : { \
405 : return x.get() OP y.get(); \
406 : } \
407 : \
408 : friend bool operator OP(float x, Float16 y) \
409 : { \
410 : return x OP y.get(); \
411 : } \
412 : \
413 : friend bool operator OP(double x, Float16 y) \
414 : { \
415 : return x OP y.get(); \
416 : } \
417 : \
418 : friend bool operator OP(int x, Float16 y) \
419 : { \
420 : return x OP y.get(); \
421 : } \
422 : \
423 : friend bool operator OP(Float16 x, float y) \
424 : { \
425 : return x.get() OP y; \
426 : } \
427 : \
428 : friend bool operator OP(Float16 x, double y) \
429 : { \
430 : return x.get() OP y; \
431 : } \
432 : \
433 : friend bool operator OP(Float16 x, int y) \
434 : { \
435 : return x.get() OP y; \
436 : }
437 :
438 421768 : GDAL_DEFINE_COMPARISON(==)
439 40409 : GDAL_DEFINE_COMPARISON(!=)
440 42740 : GDAL_DEFINE_COMPARISON(<)
441 48120 : GDAL_DEFINE_COMPARISON(>)
442 41719 : GDAL_DEFINE_COMPARISON(<=)
443 41984 : GDAL_DEFINE_COMPARISON(>=)
444 :
445 : #undef GDAL_DEFINE_COMPARISON
446 :
447 : // Standard math functions
448 :
449 40602 : friend bool isfinite(Float16 x)
450 : {
451 : using std::isfinite;
452 40602 : return isfinite(float(x));
453 : }
454 :
455 4184 : friend bool isinf(Float16 x)
456 : {
457 : using std::isinf;
458 4184 : return isinf(float(x));
459 : }
460 :
461 4753 : friend bool isnan(Float16 x)
462 : {
463 : using std::isnan;
464 4753 : return isnan(float(x));
465 : }
466 :
467 : friend bool isnormal(Float16 x)
468 : {
469 : using std::isnormal;
470 : return isnormal(float(x));
471 : }
472 :
473 0 : friend bool signbit(Float16 x)
474 : {
475 : using std::signbit;
476 0 : return signbit(float(x));
477 : }
478 :
479 201 : friend Float16 abs(Float16 x)
480 : {
481 : using std::abs;
482 201 : return Float16(abs(float(x)));
483 : }
484 :
485 201 : friend Float16 cbrt(Float16 x)
486 : {
487 : using std::cbrt;
488 201 : return Float16(cbrt(float(x)));
489 : }
490 :
491 201 : friend Float16 ceil(Float16 x)
492 : {
493 : using std::ceil;
494 201 : return Float16(ceil(float(x)));
495 : }
496 :
497 : friend Float16 copysign(Float16 x, Float16 y)
498 : {
499 : using std::copysign;
500 : return Float16(copysign(float(x), float(y)));
501 : }
502 :
503 21258 : friend Float16 fabs(Float16 x)
504 : {
505 : using std::fabs;
506 21258 : return Float16(fabs(float(x)));
507 : }
508 :
509 201 : friend Float16 floor(Float16 x)
510 : {
511 : using std::floor;
512 201 : return Float16(floor(float(x)));
513 : }
514 :
515 40401 : friend Float16 fmax(Float16 x, Float16 y)
516 : {
517 : using std::fmax;
518 40401 : return Float16(fmax(float(x), float(y)));
519 : }
520 :
521 40401 : friend Float16 fmin(Float16 x, Float16 y)
522 : {
523 : using std::fmin;
524 40401 : return Float16(fmin(float(x), float(y)));
525 : }
526 :
527 40401 : friend Float16 hypot(Float16 x, Float16 y)
528 : {
529 : using std::hypot;
530 40401 : return Float16(hypot(float(x), float(y)));
531 : }
532 :
533 40401 : friend Float16 max(Float16 x, Float16 y)
534 : {
535 : using std::max;
536 40401 : return Float16(max(float(x), float(y)));
537 : }
538 :
539 40401 : friend Float16 min(Float16 x, Float16 y)
540 : {
541 : using std::min;
542 40401 : return Float16(min(float(x), float(y)));
543 : }
544 :
545 : // Adapted from the LLVM Project, under the Apache License v2.0
546 8 : friend Float16 nextafter(Float16 x, Float16 y)
547 : {
548 8 : if (isnan(x))
549 0 : return x;
550 8 : if (isnan(y))
551 0 : return y;
552 8 : if (x == y)
553 0 : return y;
554 :
555 : std::uint16_t bits;
556 8 : if (x != Float16(0))
557 : {
558 8 : std::memcpy(&bits, &x.rValue, 2);
559 8 : if ((x < y) == (x > Float16(0)))
560 2 : ++bits;
561 : else
562 6 : --bits;
563 : }
564 : else
565 : {
566 0 : bits = (signbit(y) << 15) | 0x0001;
567 : }
568 :
569 : Float16 r;
570 8 : std::memcpy(&r.rValue, &bits, 2);
571 :
572 8 : return r;
573 : }
574 :
575 40401 : friend Float16 pow(Float16 x, Float16 y)
576 : {
577 : using std::pow;
578 40401 : return Float16(pow(float(x), float(y)));
579 : }
580 :
581 40401 : friend Float16 pow(Float16 x, int n)
582 : {
583 : using std::pow;
584 40401 : return Float16(pow(float(x), n));
585 : }
586 :
587 201 : friend Float16 round(Float16 x)
588 : {
589 : using std::round;
590 201 : return Float16(round(float(x)));
591 : }
592 :
593 101 : friend Float16 sqrt(Float16 x)
594 : {
595 : using std::sqrt;
596 101 : return Float16(sqrt(float(x)));
597 : }
598 : };
599 :
600 : template <> struct NumericLimits<Float16>
601 : {
602 : static constexpr bool is_specialized = true;
603 : static constexpr bool is_signed = true;
604 : static constexpr bool is_integer = false;
605 : static constexpr bool is_exact = false;
606 : static constexpr bool has_infinity = true;
607 : static constexpr bool has_quiet_NaN = true;
608 : static constexpr bool has_signaling_NaN = true;
609 : static constexpr bool has_denorm = true;
610 : static constexpr bool is_iec559 = true;
611 :
612 : static constexpr int digits = 11;
613 : static constexpr int digits10 = 3;
614 : static constexpr int max_digits10 = 5;
615 : static constexpr int radix = 2;
616 :
617 0 : static constexpr Float16 epsilon()
618 : {
619 0 : return Float16(Float16::make_from_bits_and_value{}, 0x1400, 0.000977f);
620 : }
621 :
622 : static constexpr Float16 min()
623 : {
624 : return Float16(Float16::make_from_bits_and_value{}, 0x0001, 6.0e-8f);
625 : }
626 :
627 1626 : static constexpr Float16 lowest()
628 : {
629 1626 : return Float16(Float16::make_from_bits_and_value{}, 0xfbff, -65504.0f);
630 : }
631 :
632 3150 : static constexpr Float16 max()
633 : {
634 3150 : return Float16(Float16::make_from_bits_and_value{}, 0x7bff, +65504.0f);
635 : }
636 :
637 5 : static constexpr Float16 infinity()
638 : {
639 5 : return Float16(Float16::make_from_bits_and_value{}, 0x7c00,
640 5 : std::numeric_limits<float>::infinity());
641 : }
642 :
643 : static constexpr Float16 quiet_NaN()
644 : {
645 : return Float16(Float16::make_from_bits_and_value{}, 0x7e00,
646 : std::numeric_limits<float>::quiet_NaN());
647 : }
648 :
649 : static constexpr Float16 signaling_NaN()
650 : {
651 : return Float16(Float16::make_from_bits_and_value{}, 0xfe00,
652 : std::numeric_limits<float>::signaling_NaN());
653 : }
654 : };
655 :
656 : //! @endcond
657 :
658 : #endif // #ifndef HAVE_STD_FLOAT16_T
659 :
660 : } // namespace cpl
661 :
662 : #ifdef HAVE_STD_FLOAT16_T
663 : using GFloat16 = std::float16_t;
664 : #else
665 : using GFloat16 = cpl::Float16;
666 : #endif
667 :
668 : // Define some GDAL wrappers. Their C equivalents are defined in `cpl_port.h`.
669 : // (These wrappers are not necessary any more in C++, one can always
670 : // call `isnan` etc directly.)
671 :
672 315520571 : template <typename T> constexpr int CPLIsNan(T x)
673 : {
674 : // We need to write `using std::isnan` instead of directly using
675 : // `std::isnan` because `std::isnan` only supports the types
676 : // `float` and `double`. The `isnan` for `cpl::Float16` is found in the
677 : // `cpl` namespace via argument-dependent lookup
678 : // <https://en.cppreference.com/w/cpp/language/adl>.
679 : using std::isnan;
680 315520571 : return isnan(x);
681 : }
682 :
683 44081 : template <typename T> constexpr int CPLIsInf(T x)
684 : {
685 : using std::isinf;
686 44081 : return isinf(x);
687 : }
688 :
689 : template <typename T> constexpr int CPLIsFinite(T x)
690 : {
691 : using std::isfinite;
692 : return isfinite(x);
693 : }
694 :
695 : #endif // #ifdef __cplusplus
696 :
697 : double CPL_DLL CPLGreatestCommonDivisor(double x, double y);
698 :
699 : #endif // CPL_FLOAT_H_INCLUDED
|
