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 2022 : constexpr Float16(make_from_bits_and_value, std::uint16_t bits,
151 : CPL_UNUSED float fValue)
152 2022 : : rValue(bits)
153 : {
154 2022 : }
155 :
156 2182288 : static unsigned float2unsigned(float f)
157 : {
158 : unsigned u;
159 2182288 : std::memcpy(&u, &f, 4);
160 2182288 : return u;
161 : }
162 :
163 2249173 : static float unsigned2float(unsigned u)
164 : {
165 : float f;
166 2249173 : std::memcpy(&f, &u, 4);
167 2249173 : return f;
168 : }
169 :
170 : // Copied from cpl_float.cpp so that we can inline for performance
171 2182288 : static std::uint16_t computeToRepr(float fFloat32)
172 : {
173 2182288 : std::uint32_t iFloat32 = float2unsigned(fFloat32);
174 :
175 2182288 : std::uint32_t iSign = (iFloat32 >> 31) & 0x00000001;
176 2182288 : std::uint32_t iExponent = (iFloat32 >> 23) & 0x000000ff;
177 2182288 : std::uint32_t iMantissa = iFloat32 & 0x007fffff;
178 :
179 2182288 : if (iExponent == 255)
180 : {
181 58941 : if (iMantissa == 0)
182 : {
183 : // Positive or negative infinity.
184 58940 : return static_cast<std::int16_t>((iSign << 15) | 0x7C00);
185 : }
186 :
187 : // NaN -- preserve sign and significand bits.
188 1 : if (iMantissa >> 13)
189 1 : return static_cast<std::int16_t>((iSign << 15) | 0x7C00 |
190 1 : (iMantissa >> 13));
191 0 : return static_cast<std::int16_t>((iSign << 15) | 0x7E00);
192 : }
193 :
194 2123345 : if (iExponent <= 127 - 15)
195 : {
196 : // Zero, float32 denormalized number or float32 too small normalized
197 : // number
198 110683 : if (13 + 1 + 127 - 15 - iExponent >= 32)
199 107373 : return static_cast<std::int16_t>(iSign << 15);
200 :
201 : // Return a denormalized number
202 3310 : return static_cast<std::int16_t>(
203 3310 : (iSign << 15) |
204 3310 : ((iMantissa | 0x00800000) >> (13 + 1 + 127 - 15 - iExponent)));
205 : }
206 :
207 2012663 : 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 1994211 : iExponent = iExponent - (127 - 15);
215 1994211 : iMantissa = iMantissa >> 13;
216 :
217 : // Assemble sign, exponent and mantissa.
218 : // coverity[overflow_sink]
219 1994211 : return static_cast<std::int16_t>((iSign << 15) | (iExponent << 10) |
220 1994211 : iMantissa);
221 : }
222 :
223 : // Copied from cpl_float.cpp so that we can inline for performance
224 2249173 : static float reprToCompute(std::uint16_t iHalf)
225 : {
226 2249173 : std::uint32_t iSign = (iHalf >> 15) & 0x00000001;
227 2249173 : int iExponent = (iHalf >> 10) & 0x0000001f;
228 2249173 : std::uint32_t iMantissa = iHalf & 0x000003ff;
229 :
230 2249173 : if (iExponent == 31)
231 : {
232 96788 : if (iMantissa == 0)
233 : {
234 : // Positive or negative infinity.
235 96783 : 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 2152387 : if (iExponent == 0)
244 : {
245 146700 : if (iMantissa == 0)
246 : {
247 : // Plus or minus zero.
248 143970 : 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 2008410 : iExponent = iExponent + (127 - 15);
264 2008410 : iMantissa = iMantissa << 13;
265 :
266 : // Assemble sign, exponent and mantissa.
267 : /* coverity[overflow_sink] */
268 2008410 : return unsigned2float((iSign << 31) |
269 2008410 : (static_cast<std::uint32_t>(iExponent) << 23) |
270 2008410 : iMantissa);
271 : }
272 :
273 2182293 : template <typename T> static repr toRepr(T fValue)
274 : {
275 2182293 : return computeToRepr(static_cast<compute>(fValue));
276 : }
277 :
278 739490 : template <typename T> static T fromRepr(repr rValue)
279 : {
280 739490 : return static_cast<T>(reprToCompute(rValue));
281 : }
282 :
283 : #endif // #ifndef HAVE__FLOAT16
284 :
285 : private:
286 : repr rValue;
287 :
288 : public:
289 1509678 : compute get() const
290 : {
291 1509678 : return reprToCompute(rValue);
292 : }
293 :
294 : Float16() = default;
295 : Float16(const Float16 &) = default;
296 : Float16(Float16 &&) = default;
297 : Float16 &operator=(const Float16 &) = default;
298 : Float16 &operator=(Float16 &&) = default;
299 :
300 : // Constructors and conversion operators
301 :
302 : #ifdef HAVE__FLOAT16
303 : // cppcheck-suppress noExplicitConstructor
304 : constexpr Float16(_Float16 hfValue) : rValue(hfValue)
305 : {
306 : }
307 :
308 : constexpr operator _Float16() const
309 : {
310 : return rValue;
311 : }
312 : #endif
313 :
314 : // cppcheck-suppress-macro noExplicitConstructor
315 : #define GDAL_DEFINE_CONVERSION(TYPE) \
316 : \
317 : Float16(TYPE fValue) : rValue(toRepr(fValue)) \
318 : { \
319 : } \
320 : \
321 : operator TYPE() const \
322 : { \
323 : return fromRepr<TYPE>(rValue); \
324 : }
325 :
326 1068304 : GDAL_DEFINE_CONVERSION(float)
327 1804929 : GDAL_DEFINE_CONVERSION(double)
328 : GDAL_DEFINE_CONVERSION(char)
329 1067 : GDAL_DEFINE_CONVERSION(signed char)
330 5315 : GDAL_DEFINE_CONVERSION(short)
331 23117 : GDAL_DEFINE_CONVERSION(int)
332 4429 : GDAL_DEFINE_CONVERSION(long)
333 4330 : GDAL_DEFINE_CONVERSION(long long)
334 1243 : GDAL_DEFINE_CONVERSION(unsigned char)
335 2333 : GDAL_DEFINE_CONVERSION(unsigned short)
336 2343 : GDAL_DEFINE_CONVERSION(unsigned int)
337 2324 : GDAL_DEFINE_CONVERSION(unsigned long)
338 2049 : GDAL_DEFINE_CONVERSION(unsigned long long)
339 :
340 : #undef GDAL_DEFINE_CONVERSION
341 :
342 : // Arithmetic operators
343 :
344 201 : friend Float16 operator+(Float16 x)
345 : {
346 201 : return +x.get();
347 : }
348 :
349 537 : friend Float16 operator-(Float16 x)
350 : {
351 537 : return -x.get();
352 : }
353 :
354 : #define GDAL_DEFINE_ARITHOP(OP) \
355 : \
356 : friend Float16 operator OP(Float16 x, Float16 y) \
357 : { \
358 : return x.get() OP y.get(); \
359 : } \
360 : \
361 : friend double operator OP(double x, Float16 y) \
362 : { \
363 : return x OP y.get(); \
364 : } \
365 : \
366 : friend float operator OP(float x, Float16 y) \
367 : { \
368 : return x OP y.get(); \
369 : } \
370 : \
371 : friend Float16 operator OP(int x, Float16 y) \
372 : { \
373 : return x OP y.get(); \
374 : } \
375 : \
376 : friend double operator OP(Float16 x, double y) \
377 : { \
378 : return x.get() OP y; \
379 : } \
380 : \
381 : friend float operator OP(Float16 x, float y) \
382 : { \
383 : return x.get() OP y; \
384 : } \
385 : \
386 : friend Float16 operator OP(Float16 x, int y) \
387 : { \
388 : return x.get() OP y; \
389 : }
390 :
391 63192 : GDAL_DEFINE_ARITHOP(+)
392 40790 : GDAL_DEFINE_ARITHOP(-)
393 40401 : GDAL_DEFINE_ARITHOP(*)
394 61257 : GDAL_DEFINE_ARITHOP(/)
395 :
396 : #undef GDAL_DEFINE_ARITHOP
397 :
398 : // Comparison operators
399 :
400 : #define GDAL_DEFINE_COMPARISON(OP) \
401 : \
402 : friend bool operator OP(Float16 x, Float16 y) \
403 : { \
404 : return x.get() OP y.get(); \
405 : } \
406 : \
407 : friend bool operator OP(float x, Float16 y) \
408 : { \
409 : return x OP y.get(); \
410 : } \
411 : \
412 : friend bool operator OP(double x, Float16 y) \
413 : { \
414 : return x OP y.get(); \
415 : } \
416 : \
417 : friend bool operator OP(int x, Float16 y) \
418 : { \
419 : return x OP y.get(); \
420 : } \
421 : \
422 : friend bool operator OP(Float16 x, float y) \
423 : { \
424 : return x.get() OP y; \
425 : } \
426 : \
427 : friend bool operator OP(Float16 x, double y) \
428 : { \
429 : return x.get() OP y; \
430 : } \
431 : \
432 : friend bool operator OP(Float16 x, int y) \
433 : { \
434 : return x.get() OP y; \
435 : }
436 :
437 421768 : GDAL_DEFINE_COMPARISON(==)
438 40409 : GDAL_DEFINE_COMPARISON(!=)
439 41364 : GDAL_DEFINE_COMPARISON(<)
440 42934 : GDAL_DEFINE_COMPARISON(>)
441 40766 : GDAL_DEFINE_COMPARISON(<=)
442 41031 : GDAL_DEFINE_COMPARISON(>=)
443 :
444 : #undef GDAL_DEFINE_COMPARISON
445 :
446 : // Standard math functions
447 :
448 40602 : friend bool isfinite(Float16 x)
449 : {
450 : using std::isfinite;
451 40602 : return isfinite(float(x));
452 : }
453 :
454 1356 : friend bool isinf(Float16 x)
455 : {
456 : using std::isinf;
457 1356 : return isinf(float(x));
458 : }
459 :
460 1897 : friend bool isnan(Float16 x)
461 : {
462 : using std::isnan;
463 1897 : return isnan(float(x));
464 : }
465 :
466 : friend bool isnormal(Float16 x)
467 : {
468 : using std::isnormal;
469 : return isnormal(float(x));
470 : }
471 :
472 0 : friend bool signbit(Float16 x)
473 : {
474 : using std::signbit;
475 0 : return signbit(float(x));
476 : }
477 :
478 201 : friend Float16 abs(Float16 x)
479 : {
480 : using std::abs;
481 201 : return Float16(abs(float(x)));
482 : }
483 :
484 201 : friend Float16 cbrt(Float16 x)
485 : {
486 : using std::cbrt;
487 201 : return Float16(cbrt(float(x)));
488 : }
489 :
490 201 : friend Float16 ceil(Float16 x)
491 : {
492 : using std::ceil;
493 201 : return Float16(ceil(float(x)));
494 : }
495 :
496 : friend Float16 copysign(Float16 x, Float16 y)
497 : {
498 : using std::copysign;
499 : return Float16(copysign(float(x), float(y)));
500 : }
501 :
502 21258 : friend Float16 fabs(Float16 x)
503 : {
504 : using std::fabs;
505 21258 : return Float16(fabs(float(x)));
506 : }
507 :
508 201 : friend Float16 floor(Float16 x)
509 : {
510 : using std::floor;
511 201 : return Float16(floor(float(x)));
512 : }
513 :
514 40401 : friend Float16 fmax(Float16 x, Float16 y)
515 : {
516 : using std::fmax;
517 40401 : return Float16(fmax(float(x), float(y)));
518 : }
519 :
520 40401 : friend Float16 fmin(Float16 x, Float16 y)
521 : {
522 : using std::fmin;
523 40401 : return Float16(fmin(float(x), float(y)));
524 : }
525 :
526 40401 : friend Float16 hypot(Float16 x, Float16 y)
527 : {
528 : using std::hypot;
529 40401 : return Float16(hypot(float(x), float(y)));
530 : }
531 :
532 40401 : friend Float16 max(Float16 x, Float16 y)
533 : {
534 : using std::max;
535 40401 : return Float16(max(float(x), float(y)));
536 : }
537 :
538 40401 : friend Float16 min(Float16 x, Float16 y)
539 : {
540 : using std::min;
541 40401 : return Float16(min(float(x), float(y)));
542 : }
543 :
544 : // Adapted from the LLVM Project, under the Apache License v2.0
545 8 : friend Float16 nextafter(Float16 x, Float16 y)
546 : {
547 8 : if (isnan(x))
548 0 : return x;
549 8 : if (isnan(y))
550 0 : return y;
551 8 : if (x == y)
552 0 : return y;
553 :
554 : std::uint16_t bits;
555 8 : if (x != Float16(0))
556 : {
557 8 : std::memcpy(&bits, &x.rValue, 2);
558 8 : if ((x < y) == (x > Float16(0)))
559 2 : ++bits;
560 : else
561 6 : --bits;
562 : }
563 : else
564 : {
565 0 : bits = (signbit(y) << 15) | 0x0001;
566 : }
567 :
568 : Float16 r;
569 8 : std::memcpy(&r.rValue, &bits, 2);
570 :
571 8 : return r;
572 : }
573 :
574 40401 : friend Float16 pow(Float16 x, Float16 y)
575 : {
576 : using std::pow;
577 40401 : return Float16(pow(float(x), float(y)));
578 : }
579 :
580 40401 : friend Float16 pow(Float16 x, int n)
581 : {
582 : using std::pow;
583 40401 : return Float16(pow(float(x), n));
584 : }
585 :
586 201 : friend Float16 round(Float16 x)
587 : {
588 : using std::round;
589 201 : return Float16(round(float(x)));
590 : }
591 :
592 101 : friend Float16 sqrt(Float16 x)
593 : {
594 : using std::sqrt;
595 101 : return Float16(sqrt(float(x)));
596 : }
597 : };
598 :
599 : template <> struct NumericLimits<Float16>
600 : {
601 : static constexpr bool is_specialized = true;
602 : static constexpr bool is_signed = true;
603 : static constexpr bool is_integer = false;
604 : static constexpr bool is_exact = false;
605 : static constexpr bool has_infinity = true;
606 : static constexpr bool has_quiet_NaN = true;
607 : static constexpr bool has_signaling_NaN = true;
608 : static constexpr bool has_denorm = true;
609 : static constexpr bool is_iec559 = true;
610 :
611 : static constexpr int digits = 11;
612 : static constexpr int digits10 = 3;
613 : static constexpr int max_digits10 = 5;
614 : static constexpr int radix = 2;
615 :
616 0 : static constexpr Float16 epsilon()
617 : {
618 0 : return Float16(Float16::make_from_bits_and_value{}, 0x1400, 0.000977f);
619 : }
620 :
621 : static constexpr Float16 min()
622 : {
623 : return Float16(Float16::make_from_bits_and_value{}, 0x0001, 6.0e-8f);
624 : }
625 :
626 674 : static constexpr Float16 lowest()
627 : {
628 674 : return Float16(Float16::make_from_bits_and_value{}, 0xfbff, -65504.0f);
629 : }
630 :
631 1346 : static constexpr Float16 max()
632 : {
633 1346 : return Float16(Float16::make_from_bits_and_value{}, 0x7bff, +65504.0f);
634 : }
635 :
636 2 : static constexpr Float16 infinity()
637 : {
638 2 : return Float16(Float16::make_from_bits_and_value{}, 0x7c00,
639 2 : std::numeric_limits<float>::infinity());
640 : }
641 :
642 : static constexpr Float16 quiet_NaN()
643 : {
644 : return Float16(Float16::make_from_bits_and_value{}, 0x7e00,
645 : std::numeric_limits<float>::quiet_NaN());
646 : }
647 :
648 : static constexpr Float16 signaling_NaN()
649 : {
650 : return Float16(Float16::make_from_bits_and_value{}, 0xfe00,
651 : std::numeric_limits<float>::signaling_NaN());
652 : }
653 : };
654 :
655 : //! @endcond
656 :
657 : #endif // #ifndef HAVE_STD_FLOAT16_T
658 :
659 : } // namespace cpl
660 :
661 : #ifdef HAVE_STD_FLOAT16_T
662 : using GFloat16 = std::float16_t;
663 : #else
664 : using GFloat16 = cpl::Float16;
665 : #endif
666 :
667 : // Define some GDAL wrappers. Their C equivalents are defined in `cpl_port.h`.
668 : // (These wrappers are not necessary any more in C++, one can always
669 : // call `isnan` etc directly.)
670 :
671 336937715 : template <typename T> constexpr int CPLIsNan(T x)
672 : {
673 : // We need to write `using std::isnan` instead of directly using
674 : // `std::isnan` because `std::isnan` only supports the types
675 : // `float` and `double`. The `isnan` for `cpl::Float16` is found in the
676 : // `cpl` namespace via argument-dependent lookup
677 : // <https://en.cppreference.com/w/cpp/language/adl>.
678 : using std::isnan;
679 336937715 : return isnan(x);
680 : }
681 :
682 43078 : template <typename T> constexpr int CPLIsInf(T x)
683 : {
684 : using std::isinf;
685 43078 : return isinf(x);
686 : }
687 :
688 : template <typename T> constexpr int CPLIsFinite(T x)
689 : {
690 : using std::isfinite;
691 : return isfinite(x);
692 : }
693 :
694 : #endif // #ifdef __cplusplus
695 :
696 : double CPL_DLL CPLGreatestCommonDivisor(double x, double y);
697 :
698 : #endif // CPL_FLOAT_H_INCLUDED
|
