GCC Code Coverage Report


Directory: ./
File: openvdb/openvdb/math/Half.h
Date: 2022-07-25 17:40:05
Exec Total Coverage
Lines: 11 12 91.7%
Functions: 1 1 100.0%
Branches: 50 92 54.3%

Line Branch Exec Source
1 ///////////////////////////////////////////////////////////////////////////
2 //
3 // Copyright (c) 2002, Industrial Light & Magic, a division of Lucas
4 // Digital Ltd. LLC
5 //
6 // All rights reserved.
7 //
8 // Redistribution and use in source and binary forms, with or without
9 // modification, are permitted provided that the following conditions are
10 // met:
11 // * Redistributions of source code must retain the above copyright
12 // notice, this list of conditions and the following disclaimer.
13 // * Redistributions in binary form must reproduce the above
14 // copyright notice, this list of conditions and the following disclaimer
15 // in the documentation and/or other materials provided with the
16 // distribution.
17 // * Neither the name of Industrial Light & Magic nor the names of
18 // its contributors may be used to endorse or promote products derived
19 // from this software without specific prior written permission.
20 //
21 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
22 // "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
23 // LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
24 // A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
25 // OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
26 // SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
27 // LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
28 // DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
29 // THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
30 // (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
31 // OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
32 //
33 ///////////////////////////////////////////////////////////////////////////
34
35 // Primary authors:
36 // Florian Kainz <kainz@ilm.com>
37 // Rod Bogart <rgb@ilm.com>
38
39 //---------------------------------------------------------------------------
40 //
41 // half -- a 16-bit floating point number class:
42 //
43 // Type half can represent positive and negative numbers whose
44 // magnitude is between roughly 6.1e-5 and 6.5e+4 with a relative
45 // error of 9.8e-4; numbers smaller than 6.1e-5 can be represented
46 // with an absolute error of 6.0e-8. All integers from -2048 to
47 // +2048 can be represented exactly.
48 //
49 // Type half behaves (almost) like the built-in C++ floating point
50 // types. In arithmetic expressions, half, float and double can be
51 // mixed freely. Here are a few examples:
52 //
53 // half a (3.5);
54 // float b (a + sqrt (a));
55 // a += b;
56 // b += a;
57 // b = a + 7;
58 //
59 // Conversions from half to float are lossless; all half numbers
60 // are exactly representable as floats.
61 //
62 // Conversions from float to half may not preserve a float's value
63 // exactly. If a float is not representable as a half, then the
64 // float value is rounded to the nearest representable half. If a
65 // float value is exactly in the middle between the two closest
66 // representable half values, then the float value is rounded to
67 // the closest half whose least significant bit is zero.
68 //
69 // Overflows during float-to-half conversions cause arithmetic
70 // exceptions. An overflow occurs when the float value to be
71 // converted is too large to be represented as a half, or if the
72 // float value is an infinity or a NAN.
73 //
74 // The implementation of type half makes the following assumptions
75 // about the implementation of the built-in C++ types:
76 //
77 // float is an IEEE 754 single-precision number
78 // sizeof (float) == 4
79 // sizeof (unsigned int) == sizeof (float)
80 // alignof (unsigned int) == alignof (float)
81 // sizeof (unsigned short) == 2
82 //
83 //---------------------------------------------------------------------------
84
85 #ifndef OPENVDB_MATH_HALF_HAS_BEEN_INCLUDED
86 #define OPENVDB_MATH_HALF_HAS_BEEN_INCLUDED
87
88 #include <openvdb/Platform.h>
89 #include <openvdb/version.h>
90 #include <iostream>
91
92 namespace openvdb {
93 OPENVDB_USE_VERSION_NAMESPACE
94 namespace OPENVDB_VERSION_NAME {
95 namespace math {
96 namespace internal {
97
98 class OPENVDB_API half
99 {
100 public:
101
102 //-------------
103 // Constructors
104 //-------------
105
106 half () = default; // no initialization
107 half (float f);
108 // rule of 5
109 ~half () = default;
110 half (const half &) = default;
111 half (half &&) noexcept = default;
112
113 //--------------------
114 // Conversion to float
115 //--------------------
116
117 operator float () const;
118
119
120 //------------
121 // Unary minus
122 //------------
123
124 half operator - () const;
125
126
127 //-----------
128 // Assignment
129 //-----------
130
131 half & operator = (const half &h) = default;
132 half & operator = (half &&h) noexcept = default;
133 half & operator = (float f);
134
135 half & operator += (half h);
136 half & operator += (float f);
137
138 half & operator -= (half h);
139 half & operator -= (float f);
140
141 half & operator *= (half h);
142 half & operator *= (float f);
143
144 half & operator /= (half h);
145 half & operator /= (float f);
146
147
148 //---------------------------------------------------------
149 // Round to n-bit precision (n should be between 0 and 10).
150 // After rounding, the significand's 10-n least significant
151 // bits will be zero.
152 //---------------------------------------------------------
153
154 half round (unsigned int n) const;
155
156
157 //--------------------------------------------------------------------
158 // Classification:
159 //
160 // h.isFinite() returns true if h is a normalized number,
161 // a denormalized number or zero
162 //
163 // h.isNormalized() returns true if h is a normalized number
164 //
165 // h.isDenormalized() returns true if h is a denormalized number
166 //
167 // h.isZero() returns true if h is zero
168 //
169 // h.isNan() returns true if h is a NAN
170 //
171 // h.isInfinity() returns true if h is a positive
172 // or a negative infinity
173 //
174 // h.isNegative() returns true if the sign bit of h
175 // is set (negative)
176 //--------------------------------------------------------------------
177
178 bool isFinite () const;
179 bool isNormalized () const;
180 bool isDenormalized () const;
181 bool isZero () const;
182 bool isNan () const;
183 bool isInfinity () const;
184 bool isNegative () const;
185
186
187 //--------------------------------------------
188 // Special values
189 //
190 // posInf() returns +infinity
191 //
192 // negInf() returns -infinity
193 //
194 // qNan() returns a NAN with the bit
195 // pattern 0111111111111111
196 //
197 // sNan() returns a NAN with the bit
198 // pattern 0111110111111111
199 //--------------------------------------------
200
201 static half posInf ();
202 static half negInf ();
203 static half qNan ();
204 static half sNan ();
205
206
207 //--------------------------------------
208 // Access to the internal representation
209 //--------------------------------------
210
211 unsigned short bits () const;
212 void setBits (unsigned short bits);
213
214
215 public:
216
217 union uif
218 {
219 unsigned int i;
220 float f;
221 };
222
223 private:
224
225 static short convert (int i);
226 static float overflow ();
227
228 unsigned short _h;
229
230 static const uif _toFloat[1 << 16];
231 static const unsigned short _eLut[1 << 9];
232 };
233
234
235
236 //-----------
237 // Stream I/O
238 //-----------
239
240 OPENVDB_API std::ostream & operator << (std::ostream &os, half h);
241 OPENVDB_API std::istream & operator >> (std::istream &is, half &h);
242
243
244 //----------
245 // Debugging
246 //----------
247
248 OPENVDB_API void printBits (std::ostream &os, half h);
249 OPENVDB_API void printBits (std::ostream &os, float f);
250 OPENVDB_API void printBits (char c[19], half h);
251 OPENVDB_API void printBits (char c[35], float f);
252
253
254 //-------------------------------------------------------------------------
255 // Limits
256 //
257 // Visual C++ will complain if HALF_MIN, HALF_NRM_MIN etc. are not float
258 // constants, but at least one other compiler (gcc 2.96) produces incorrect
259 // results if they are.
260 //-------------------------------------------------------------------------
261
262 #if (defined _WIN32 || defined _WIN64) && defined _MSC_VER
263
264 #define VDB_HALF_MIN 5.96046448e-08f // Smallest positive half
265
266 #define VDB_HALF_NRM_MIN 6.10351562e-05f // Smallest positive normalized half
267
268 #define VDB_HALF_MAX 65504.0f // Largest positive half
269
270 #define VDB_HALF_EPSILON 0.00097656f // Smallest positive e for which
271 // half (1.0 + e) != half (1.0)
272 #else
273
274 #define VDB_HALF_MIN 5.96046448e-08 // Smallest positive half
275
276 #define VDB_HALF_NRM_MIN 6.10351562e-05 // Smallest positive normalized half
277
278 #define VDB_HALF_MAX 65504.0 // Largest positive half
279
280 #define VDB_HALF_EPSILON 0.00097656 // Smallest positive e for which
281 // half (1.0 + e) != half (1.0)
282 #endif
283
284
285 #define VDB_HALF_MANT_DIG 11 // Number of digits in mantissa
286 // (significand + hidden leading 1)
287
288 //
289 // floor( (VDB_HALF_MANT_DIG - 1) * log10(2) ) => 3.01... -> 3
290 #define VDB_HALF_DIG 3 // Number of base 10 digits that
291 // can be represented without change
292
293 // ceil(VDB_HALF_MANT_DIG * log10(2) + 1) => 4.31... -> 5
294 #define VDB_HALF_DECIMAL_DIG 5 // Number of base-10 digits that are
295 // necessary to uniquely represent all
296 // distinct values
297
298 #define VDB_HALF_RADIX 2 // Base of the exponent
299
300 #define VDB_HALF_MIN_EXP -13 // Minimum negative integer such that
301 // HALF_RADIX raised to the power of
302 // one less than that integer is a
303 // normalized half
304
305 #define VDB_HALF_MAX_EXP 16 // Maximum positive integer such that
306 // HALF_RADIX raised to the power of
307 // one less than that integer is a
308 // normalized half
309
310 #define VDB_HALF_MIN_10_EXP -4 // Minimum positive integer such
311 // that 10 raised to that power is
312 // a normalized half
313
314 #define VDB_HALF_MAX_10_EXP 4 // Maximum positive integer such
315 // that 10 raised to that power is
316 // a normalized half
317
318
319 //---------------------------------------------------------------------------
320 //
321 // Implementation --
322 //
323 // Representation of a float:
324 //
325 // We assume that a float, f, is an IEEE 754 single-precision
326 // floating point number, whose bits are arranged as follows:
327 //
328 // 31 (msb)
329 // |
330 // | 30 23
331 // | | |
332 // | | | 22 0 (lsb)
333 // | | | | |
334 // X XXXXXXXX XXXXXXXXXXXXXXXXXXXXXXX
335 //
336 // s e m
337 //
338 // S is the sign-bit, e is the exponent and m is the significand.
339 //
340 // If e is between 1 and 254, f is a normalized number:
341 //
342 // s e-127
343 // f = (-1) * 2 * 1.m
344 //
345 // If e is 0, and m is not zero, f is a denormalized number:
346 //
347 // s -126
348 // f = (-1) * 2 * 0.m
349 //
350 // If e and m are both zero, f is zero:
351 //
352 // f = 0.0
353 //
354 // If e is 255, f is an "infinity" or "not a number" (NAN),
355 // depending on whether m is zero or not.
356 //
357 // Examples:
358 //
359 // 0 00000000 00000000000000000000000 = 0.0
360 // 0 01111110 00000000000000000000000 = 0.5
361 // 0 01111111 00000000000000000000000 = 1.0
362 // 0 10000000 00000000000000000000000 = 2.0
363 // 0 10000000 10000000000000000000000 = 3.0
364 // 1 10000101 11110000010000000000000 = -124.0625
365 // 0 11111111 00000000000000000000000 = +infinity
366 // 1 11111111 00000000000000000000000 = -infinity
367 // 0 11111111 10000000000000000000000 = NAN
368 // 1 11111111 11111111111111111111111 = NAN
369 //
370 // Representation of a half:
371 //
372 // Here is the bit-layout for a half number, h:
373 //
374 // 15 (msb)
375 // |
376 // | 14 10
377 // | | |
378 // | | | 9 0 (lsb)
379 // | | | | |
380 // X XXXXX XXXXXXXXXX
381 //
382 // s e m
383 //
384 // S is the sign-bit, e is the exponent and m is the significand.
385 //
386 // If e is between 1 and 30, h is a normalized number:
387 //
388 // s e-15
389 // h = (-1) * 2 * 1.m
390 //
391 // If e is 0, and m is not zero, h is a denormalized number:
392 //
393 // S -14
394 // h = (-1) * 2 * 0.m
395 //
396 // If e and m are both zero, h is zero:
397 //
398 // h = 0.0
399 //
400 // If e is 31, h is an "infinity" or "not a number" (NAN),
401 // depending on whether m is zero or not.
402 //
403 // Examples:
404 //
405 // 0 00000 0000000000 = 0.0
406 // 0 01110 0000000000 = 0.5
407 // 0 01111 0000000000 = 1.0
408 // 0 10000 0000000000 = 2.0
409 // 0 10000 1000000000 = 3.0
410 // 1 10101 1111000001 = -124.0625
411 // 0 11111 0000000000 = +infinity
412 // 1 11111 0000000000 = -infinity
413 // 0 11111 1000000000 = NAN
414 // 1 11111 1111111111 = NAN
415 //
416 // Conversion:
417 //
418 // Converting from a float to a half requires some non-trivial bit
419 // manipulations. In some cases, this makes conversion relatively
420 // slow, but the most common case is accelerated via table lookups.
421 //
422 // Converting back from a half to a float is easier because we don't
423 // have to do any rounding. In addition, there are only 65536
424 // different half numbers; we can convert each of those numbers once
425 // and store the results in a table. Later, all conversions can be
426 // done using only simple table lookups.
427 //
428 //---------------------------------------------------------------------------
429
430
431 //----------------------------
432 // Half-from-float constructor
433 //----------------------------
434
435 inline
436 6676871 half::half (float f)
437 {
438 uif x;
439
440 x.f = f;
441
442
2/2
✓ Branch 0 taken 11629 times.
✓ Branch 1 taken 6665217 times.
6676846 if (f == 0)
443 {
444 //
445 // Common special case - zero.
446 // Preserve the zero's sign bit.
447 //
448
449
3/6
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✓ Branch 6 taken 1 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 1 times.
✗ Branch 10 not taken.
11632 _h = (unsigned short)(x.i >> 16);
450 }
451 else
452 {
453 //
454 // We extract the combined sign and exponent, e, from our
455 // floating-point number, f. Then we convert e to the sign
456 // and exponent of the half number via a table lookup.
457 //
458 // For the most common case, where a normalized half is produced,
459 // the table lookup returns a non-zero value; in this case, all
460 // we have to do is round f's significand to 10 bits and combine
461 // the result with e.
462 //
463 // For all other cases (overflow, zeroes, denormalized numbers
464 // resulting from underflow, infinities and NANs), the table
465 // lookup returns zero, and we call a longer, non-inline function
466 // to do the float-to-half conversion.
467 //
468
469 6665217 int e = (x.i >> 23) & 0x000001ff;
470
471 6665236 e = _eLut[e];
472
473
18/26
✓ Branch 0 taken 6665196 times.
✓ Branch 1 taken 1 times.
✓ Branch 2 taken 1 times.
✓ Branch 3 taken 1 times.
✓ Branch 4 taken 1 times.
✓ Branch 5 taken 1 times.
✓ Branch 6 taken 2 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 2 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 2 times.
✗ Branch 12 not taken.
✓ Branch 13 taken 2 times.
✓ Branch 14 taken 1 times.
✓ Branch 15 taken 1 times.
✓ Branch 16 taken 13 times.
✓ Branch 17 taken 10 times.
✓ Branch 18 taken 1 times.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
✓ Branch 21 taken 1 times.
✓ Branch 22 taken 1 times.
✗ Branch 23 not taken.
✗ Branch 24 not taken.
✓ Branch 25 taken 1 times.
6665238 if (e)
474 {
475 //
476 // Simple case - round the significand, m, to 10
477 // bits and combine it with the sign and exponent.
478 //
479
480 6665207 int m = x.i & 0x007fffff;
481 6665218 _h = (unsigned short)(e + ((m + 0x00000fff + ((m >> 13) & 1)) >> 13));
482 }
483 else
484 {
485 //
486 // Difficult case - call a function.
487 //
488
489
8/26
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1 times.
✗ Branch 8 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✓ Branch 16 taken 2 times.
✗ Branch 17 not taken.
✓ Branch 19 taken 2 times.
✗ Branch 20 not taken.
✓ Branch 22 taken 1 times.
✗ Branch 23 not taken.
✗ Branch 25 not taken.
✗ Branch 26 not taken.
✗ Branch 28 not taken.
✗ Branch 29 not taken.
✓ Branch 31 taken 1 times.
✗ Branch 32 not taken.
✗ Branch 34 not taken.
✗ Branch 35 not taken.
✓ Branch 37 taken 1 times.
✗ Branch 38 not taken.
20 _h = convert (x.i);
490 }
491 }
492 6676846 }
493
494
495 //------------------------------------------
496 // Half-to-float conversion via table lookup
497 //------------------------------------------
498
499 inline
500 half::operator float () const
501 {
502
19/32
✓ Branch 0 taken 8 times.
✓ Branch 1 taken 1 times.
✓ Branch 2 taken 8 times.
✓ Branch 3 taken 1 times.
✓ Branch 4 taken 1 times.
✓ Branch 5 taken 8 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 25 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 12 times.
✓ Branch 10 taken 1 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 1 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 1 times.
✗ Branch 17 not taken.
✓ Branch 19 taken 1 times.
✗ Branch 20 not taken.
✓ Branch 22 taken 1 times.
✗ Branch 23 not taken.
✓ Branch 25 taken 1 times.
✗ Branch 26 not taken.
✓ Branch 28 taken 1 times.
✗ Branch 29 not taken.
✓ Branch 31 taken 1 times.
✗ Branch 32 not taken.
✓ Branch 34 taken 1 times.
✗ Branch 35 not taken.
✓ Branch 37 taken 1 times.
✗ Branch 38 not taken.
✓ Branch 40 taken 1 times.
✗ Branch 41 not taken.
524347 return _toFloat[_h].f;
503 }
504
505
506 //-------------------------
507 // Round to n-bit precision
508 //-------------------------
509
510 inline half
511 half::round (unsigned int n) const
512 {
513 //
514 // Parameter check.
515 //
516
517 if (n >= 10)
518 return *this;
519
520 //
521 // Disassemble h into the sign, s,
522 // and the combined exponent and significand, e.
523 //
524
525 unsigned short s = _h & 0x8000;
526 unsigned short e = _h & 0x7fff;
527
528 //
529 // Round the exponent and significand to the nearest value
530 // where ones occur only in the (10-n) most significant bits.
531 // Note that the exponent adjusts automatically if rounding
532 // up causes the significand to overflow.
533 //
534
535 e = (unsigned short)(e >> (9 - n));
536 e = (unsigned short)(e + (e & 1));
537 e = (unsigned short)(e << (9 - n));
538
539 //
540 // Check for exponent overflow.
541 //
542
543 if (e >= 0x7c00)
544 {
545 //
546 // Overflow occurred -- truncate instead of rounding.
547 //
548
549 e = _h;
550 e = (unsigned short)(e >> (10 - n));
551 e = (unsigned short)(e << (10 - n));
552 }
553
554 //
555 // Put the original sign bit back.
556 //
557
558 half h;
559 h._h = (unsigned short)(s | e);
560
561 return h;
562 }
563
564
565 //-----------------------
566 // Other inline functions
567 //-----------------------
568
569 inline half
570 half::operator - () const
571 {
572 half h;
573 h._h = _h ^ 0x8000;
574 return h;
575 }
576
577
578 inline half &
579 half::operator = (float f)
580 {
581 *this = half (f);
582 return *this;
583 }
584
585
586 inline half &
587 half::operator += (half h)
588 {
589 *this = half (float (*this) + float (h));
590 return *this;
591 }
592
593
594 inline half &
595 half::operator += (float f)
596 {
597 *this = half (float (*this) + f);
598 return *this;
599 }
600
601
602 inline half &
603 half::operator -= (half h)
604 {
605 *this = half (float (*this) - float (h));
606 return *this;
607 }
608
609
610 inline half &
611 half::operator -= (float f)
612 {
613 *this = half (float (*this) - f);
614 return *this;
615 }
616
617
618 inline half &
619 half::operator *= (half h)
620 {
621 *this = half (float (*this) * float (h));
622 return *this;
623 }
624
625
626 inline half &
627 half::operator *= (float f)
628 {
629 *this = half (float (*this) * f);
630 return *this;
631 }
632
633
634 inline half &
635 half::operator /= (half h)
636 {
637 *this = half (float (*this) / float (h));
638 return *this;
639 }
640
641
642 inline half &
643 half::operator /= (float f)
644 {
645 *this = half (float (*this) / f);
646 return *this;
647 }
648
649
650 inline bool
651 half::isFinite () const
652 {
653 unsigned short e = (_h >> 10) & 0x001f;
654 return e < 31;
655 }
656
657
658 inline bool
659 half::isNormalized () const
660 {
661 unsigned short e = (_h >> 10) & 0x001f;
662 return e > 0 && e < 31;
663 }
664
665
666 inline bool
667 half::isDenormalized () const
668 {
669 unsigned short e = (_h >> 10) & 0x001f;
670 unsigned short m = _h & 0x3ff;
671 return e == 0 && m != 0;
672 }
673
674
675 inline bool
676 half::isZero () const
677 {
678 return (_h & 0x7fff) == 0;
679 }
680
681
682 inline bool
683 half::isNan () const
684 {
685 unsigned short e = (_h >> 10) & 0x001f;
686 unsigned short m = _h & 0x3ff;
687 return e == 31 && m != 0;
688 }
689
690
691 inline bool
692 half::isInfinity () const
693 {
694 unsigned short e = (_h >> 10) & 0x001f;
695 unsigned short m = _h & 0x3ff;
696 return e == 31 && m == 0;
697 }
698
699
700 inline bool
701 half::isNegative () const
702 {
703 return (_h & 0x8000) != 0;
704 }
705
706
707 inline half
708 half::posInf ()
709 {
710 half h;
711 h._h = 0x7c00;
712 return h;
713 }
714
715
716 inline half
717 half::negInf ()
718 {
719 half h;
720 h._h = 0xfc00;
721 return h;
722 }
723
724
725 inline half
726 half::qNan ()
727 {
728 half h;
729 h._h = 0x7fff;
730 return h;
731 }
732
733
734 inline half
735 half::sNan ()
736 {
737 half h;
738 h._h = 0x7dff;
739 return h;
740 }
741
742
743 inline unsigned short
744 half::bits () const
745 {
746 return _h;
747 }
748
749
750 inline void
751 half::setBits (unsigned short bits)
752 {
753 _h = bits;
754 }
755
756 } // namespace internal
757 } // namespace math
758 } // namespace OPENVDB_VERSION_NAME
759 } // namespace openvdb
760
761 #endif // OPENVDB_MATH_HALF_HAS_BEEN_INCLUDED
762