PBRT
|
00001 00002 // 00003 // Copyright (c) 2002, Industrial Light & Magic, a division of Lucas 00004 // Digital Ltd. LLC 00005 // 00006 // All rights reserved. 00007 // 00008 // Redistribution and use in source and binary forms, with or without 00009 // modification, are permitted provided that the following conditions are 00010 // met: 00011 // * Redistributions of source code must retain the above copyright 00012 // notice, this list of conditions and the following disclaimer. 00013 // * Redistributions in binary form must reproduce the above 00014 // copyright notice, this list of conditions and the following disclaimer 00015 // in the documentation and/or other materials provided with the 00016 // distribution. 00017 // * Neither the name of Industrial Light & Magic nor the names of 00018 // its contributors may be used to endorse or promote products derived 00019 // from this software without specific prior written permission. 00020 // 00021 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS 00022 // "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT 00023 // LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR 00024 // A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT 00025 // OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, 00026 // SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT 00027 // LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, 00028 // DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY 00029 // THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT 00030 // (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE 00031 // OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 00032 // 00034 00035 // Primary authors: 00036 // Florian Kainz <kainz@ilm.com> 00037 // Rod Bogart <rgb@ilm.com> 00038 00039 //--------------------------------------------------------------------------- 00040 // 00041 // half -- a 16-bit floating point number class: 00042 // 00043 // Type half can represent positive and negative numbers whose 00044 // magnitude is between roughly 6.1e-5 and 6.5e+4 with a relative 00045 // error of 9.8e-4; numbers smaller than 6.1e-5 can be represented 00046 // with an absolute error of 6.0e-8. All integers from -2048 to 00047 // +2048 can be represented exactly. 00048 // 00049 // Type half behaves (almost) like the built-in C++ floating point 00050 // types. In arithmetic expressions, half, float and double can be 00051 // mixed freely. Here are a few examples: 00052 // 00053 // half a (3.5); 00054 // float b (a + sqrt (a)); 00055 // a += b; 00056 // b += a; 00057 // b = a + 7; 00058 // 00059 // Conversions from half to float are lossless; all half numbers 00060 // are exactly representable as floats. 00061 // 00062 // Conversions from float to half may not preserve a float's value 00063 // exactly. If a float is not representable as a half, then the 00064 // float value is rounded to the nearest representable half. If a 00065 // float value is exactly in the middle between the two closest 00066 // representable half values, then the float value is rounded to 00067 // the closest half whose least significant bit is zero. 00068 // 00069 // Overflows during float-to-half conversions cause arithmetic 00070 // exceptions. An overflow occurs when the float value to be 00071 // converted is too large to be represented as a half, or if the 00072 // float value is an infinity or a NAN. 00073 // 00074 // The implementation of type half makes the following assumptions 00075 // about the implementation of the built-in C++ types: 00076 // 00077 // float is an IEEE 754 single-precision number 00078 // sizeof (float) == 4 00079 // sizeof (unsigned int) == sizeof (float) 00080 // alignof (unsigned int) == alignof (float) 00081 // sizeof (unsigned short) == 2 00082 // 00083 //--------------------------------------------------------------------------- 00084 00085 #ifndef _HALF_H_ 00086 #define _HALF_H_ 00087 00088 #include <iostream> 00089 00090 #if defined(OPENEXR_DLL) 00091 #if defined(HALF_EXPORTS) 00092 #define HALF_EXPORT __declspec(dllexport) 00093 #else 00094 #define HALF_EXPORT __declspec(dllimport) 00095 #endif 00096 #define HALF_EXPORT_CONST 00097 #else 00098 #define HALF_EXPORT 00099 #define HALF_EXPORT_CONST const 00100 #endif 00101 00102 class HALF_EXPORT half 00103 { 00104 public: 00105 00106 //------------- 00107 // Constructors 00108 //------------- 00109 00110 half (); // no initialization 00111 half (float f); 00112 00113 00114 //-------------------- 00115 // Conversion to float 00116 //-------------------- 00117 00118 operator float () const; 00119 00120 00121 //------------ 00122 // Unary minus 00123 //------------ 00124 00125 half operator - () const; 00126 00127 00128 //----------- 00129 // Assignment 00130 //----------- 00131 00132 half & operator = (half h); 00133 half & operator = (float f); 00134 00135 half & operator += (half h); 00136 half & operator += (float f); 00137 00138 half & operator -= (half h); 00139 half & operator -= (float f); 00140 00141 half & operator *= (half h); 00142 half & operator *= (float f); 00143 00144 half & operator /= (half h); 00145 half & operator /= (float f); 00146 00147 00148 //--------------------------------------------------------- 00149 // Round to n-bit precision (n should be between 0 and 10). 00150 // After rounding, the significand's 10-n least significant 00151 // bits will be zero. 00152 //--------------------------------------------------------- 00153 00154 half round (unsigned int n) const; 00155 00156 00157 //-------------------------------------------------------------------- 00158 // Classification: 00159 // 00160 // h.isFinite() returns true if h is a normalized number, 00161 // a denormalized number or zero 00162 // 00163 // h.isNormalized() returns true if h is a normalized number 00164 // 00165 // h.isDenormalized() returns true if h is a denormalized number 00166 // 00167 // h.isZero() returns true if h is zero 00168 // 00169 // h.isNan() returns true if h is a NAN 00170 // 00171 // h.isInfinity() returns true if h is a positive 00172 // or a negative infinity 00173 // 00174 // h.isNegative() returns true if the sign bit of h 00175 // is set (negative) 00176 //-------------------------------------------------------------------- 00177 00178 bool isFinite () const; 00179 bool isNormalized () const; 00180 bool isDenormalized () const; 00181 bool isZero () const; 00182 bool isNan () const; 00183 bool isInfinity () const; 00184 bool isNegative () const; 00185 00186 00187 //-------------------------------------------- 00188 // Special values 00189 // 00190 // posInf() returns +infinity 00191 // 00192 // negInf() returns -infinity 00193 // 00194 // qNan() returns a NAN with the bit 00195 // pattern 0111111111111111 00196 // 00197 // sNan() returns a NAN with the bit 00198 // pattern 0111110111111111 00199 //-------------------------------------------- 00200 00201 static half posInf (); 00202 static half negInf (); 00203 static half qNan (); 00204 static half sNan (); 00205 00206 00207 //-------------------------------------- 00208 // Access to the internal representation 00209 //-------------------------------------- 00210 00211 unsigned short bits () const; 00212 void setBits (unsigned short bits); 00213 00214 00215 public: 00216 00217 union uif 00218 { 00219 unsigned int i; 00220 float f; 00221 }; 00222 00223 private: 00224 00225 static short convert (int i); 00226 static float overflow (); 00227 00228 unsigned short _h; 00229 00230 static HALF_EXPORT_CONST uif _toFloat[1 << 16]; 00231 static HALF_EXPORT_CONST unsigned short _eLut[1 << 9]; 00232 }; 00233 00234 //----------- 00235 // Stream I/O 00236 //----------- 00237 00238 HALF_EXPORT std::ostream & operator << (std::ostream &os, half h); 00239 HALF_EXPORT std::istream & operator >> (std::istream &is, half &h); 00240 00241 00242 //---------- 00243 // Debugging 00244 //---------- 00245 00246 HALF_EXPORT void printBits (std::ostream &os, half h); 00247 HALF_EXPORT void printBits (std::ostream &os, float f); 00248 HALF_EXPORT void printBits (char c[19], half h); 00249 HALF_EXPORT void printBits (char c[35], float f); 00250 00251 00252 //------------------------------------------------------------------------- 00253 // Limits 00254 // 00255 // Visual C++ will complain if HALF_MIN, HALF_NRM_MIN etc. are not float 00256 // constants, but at least one other compiler (gcc 2.96) produces incorrect 00257 // results if they are. 00258 //------------------------------------------------------------------------- 00259 00260 #if (defined _WIN32 || defined _WIN64) && defined _MSC_VER 00261 00262 #define HALF_MIN 5.96046448e-08f // Smallest positive half 00263 00264 #define HALF_NRM_MIN 6.10351562e-05f // Smallest positive normalized half 00265 00266 #define HALF_MAX 65504.0f // Largest positive half 00267 00268 #define HALF_EPSILON 0.00097656f // Smallest positive e for which 00269 // half (1.0 + e) != half (1.0) 00270 #else 00271 00272 #define HALF_MIN 5.96046448e-08 // Smallest positive half 00273 00274 #define HALF_NRM_MIN 6.10351562e-05 // Smallest positive normalized half 00275 00276 #define HALF_MAX 65504.0 // Largest positive half 00277 00278 #define HALF_EPSILON 0.00097656 // Smallest positive e for which 00279 // half (1.0 + e) != half (1.0) 00280 #endif 00281 00282 00283 #define HALF_MANT_DIG 11 // Number of digits in mantissa 00284 // (significand + hidden leading 1) 00285 00286 #define HALF_DIG 2 // Number of base 10 digits that 00287 // can be represented without change 00288 00289 #define HALF_RADIX 2 // Base of the exponent 00290 00291 #define HALF_MIN_EXP -13 // Minimum negative integer such that 00292 // HALF_RADIX raised to the power of 00293 // one less than that integer is a 00294 // normalized half 00295 00296 #define HALF_MAX_EXP 16 // Maximum positive integer such that 00297 // HALF_RADIX raised to the power of 00298 // one less than that integer is a 00299 // normalized half 00300 00301 #define HALF_MIN_10_EXP -4 // Minimum positive integer such 00302 // that 10 raised to that power is 00303 // a normalized half 00304 00305 #define HALF_MAX_10_EXP 4 // Maximum positive integer such 00306 // that 10 raised to that power is 00307 // a normalized half 00308 00309 00310 //--------------------------------------------------------------------------- 00311 // 00312 // Implementation -- 00313 // 00314 // Representation of a float: 00315 // 00316 // We assume that a float, f, is an IEEE 754 single-precision 00317 // floating point number, whose bits are arranged as follows: 00318 // 00319 // 31 (msb) 00320 // | 00321 // | 30 23 00322 // | | | 00323 // | | | 22 0 (lsb) 00324 // | | | | | 00325 // X XXXXXXXX XXXXXXXXXXXXXXXXXXXXXXX 00326 // 00327 // s e m 00328 // 00329 // S is the sign-bit, e is the exponent and m is the significand. 00330 // 00331 // If e is between 1 and 254, f is a normalized number: 00332 // 00333 // s e-127 00334 // f = (-1) * 2 * 1.m 00335 // 00336 // If e is 0, and m is not zero, f is a denormalized number: 00337 // 00338 // s -126 00339 // f = (-1) * 2 * 0.m 00340 // 00341 // If e and m are both zero, f is zero: 00342 // 00343 // f = 0.0 00344 // 00345 // If e is 255, f is an "infinity" or "not a number" (NAN), 00346 // depending on whether m is zero or not. 00347 // 00348 // Examples: 00349 // 00350 // 0 00000000 00000000000000000000000 = 0.0 00351 // 0 01111110 00000000000000000000000 = 0.5 00352 // 0 01111111 00000000000000000000000 = 1.0 00353 // 0 10000000 00000000000000000000000 = 2.0 00354 // 0 10000000 10000000000000000000000 = 3.0 00355 // 1 10000101 11110000010000000000000 = -124.0625 00356 // 0 11111111 00000000000000000000000 = +infinity 00357 // 1 11111111 00000000000000000000000 = -infinity 00358 // 0 11111111 10000000000000000000000 = NAN 00359 // 1 11111111 11111111111111111111111 = NAN 00360 // 00361 // Representation of a half: 00362 // 00363 // Here is the bit-layout for a half number, h: 00364 // 00365 // 15 (msb) 00366 // | 00367 // | 14 10 00368 // | | | 00369 // | | | 9 0 (lsb) 00370 // | | | | | 00371 // X XXXXX XXXXXXXXXX 00372 // 00373 // s e m 00374 // 00375 // S is the sign-bit, e is the exponent and m is the significand. 00376 // 00377 // If e is between 1 and 30, h is a normalized number: 00378 // 00379 // s e-15 00380 // h = (-1) * 2 * 1.m 00381 // 00382 // If e is 0, and m is not zero, h is a denormalized number: 00383 // 00384 // S -14 00385 // h = (-1) * 2 * 0.m 00386 // 00387 // If e and m are both zero, h is zero: 00388 // 00389 // h = 0.0 00390 // 00391 // If e is 31, h is an "infinity" or "not a number" (NAN), 00392 // depending on whether m is zero or not. 00393 // 00394 // Examples: 00395 // 00396 // 0 00000 0000000000 = 0.0 00397 // 0 01110 0000000000 = 0.5 00398 // 0 01111 0000000000 = 1.0 00399 // 0 10000 0000000000 = 2.0 00400 // 0 10000 1000000000 = 3.0 00401 // 1 10101 1111000001 = -124.0625 00402 // 0 11111 0000000000 = +infinity 00403 // 1 11111 0000000000 = -infinity 00404 // 0 11111 1000000000 = NAN 00405 // 1 11111 1111111111 = NAN 00406 // 00407 // Conversion: 00408 // 00409 // Converting from a float to a half requires some non-trivial bit 00410 // manipulations. In some cases, this makes conversion relatively 00411 // slow, but the most common case is accelerated via table lookups. 00412 // 00413 // Converting back from a half to a float is easier because we don't 00414 // have to do any rounding. In addition, there are only 65536 00415 // different half numbers; we can convert each of those numbers once 00416 // and store the results in a table. Later, all conversions can be 00417 // done using only simple table lookups. 00418 // 00419 //--------------------------------------------------------------------------- 00420 00421 00422 //-------------------- 00423 // Simple constructors 00424 //-------------------- 00425 00426 inline 00427 half::half () 00428 { 00429 // no initialization 00430 } 00431 00432 00433 //---------------------------- 00434 // Half-from-float constructor 00435 //---------------------------- 00436 00437 inline 00438 half::half (float f) 00439 { 00440 uif x; 00441 00442 x.f = f; 00443 00444 if (f == 0) 00445 { 00446 // 00447 // Common special case - zero. 00448 // Preserve the zero's sign bit. 00449 // 00450 00451 _h = (x.i >> 16); 00452 } 00453 else 00454 { 00455 // 00456 // We extract the combined sign and exponent, e, from our 00457 // floating-point number, f. Then we convert e to the sign 00458 // and exponent of the half number via a table lookup. 00459 // 00460 // For the most common case, where a normalized half is produced, 00461 // the table lookup returns a non-zero value; in this case, all 00462 // we have to do is round f's significand to 10 bits and combine 00463 // the result with e. 00464 // 00465 // For all other cases (overflow, zeroes, denormalized numbers 00466 // resulting from underflow, infinities and NANs), the table 00467 // lookup returns zero, and we call a longer, non-inline function 00468 // to do the float-to-half conversion. 00469 // 00470 00471 register int e = (x.i >> 23) & 0x000001ff; 00472 00473 e = _eLut[e]; 00474 00475 if (e) 00476 { 00477 // 00478 // Simple case - round the significand, m, to 10 00479 // bits and combine it with the sign and exponent. 00480 // 00481 00482 register int m = x.i & 0x007fffff; 00483 _h = e + ((m + 0x00000fff + ((m >> 13) & 1)) >> 13); 00484 } 00485 else 00486 { 00487 // 00488 // Difficult case - call a function. 00489 // 00490 00491 _h = convert (x.i); 00492 } 00493 } 00494 } 00495 00496 00497 //------------------------------------------ 00498 // Half-to-float conversion via table lookup 00499 //------------------------------------------ 00500 00501 inline 00502 half::operator float () const 00503 { 00504 return _toFloat[_h].f; 00505 } 00506 00507 00508 //------------------------- 00509 // Round to n-bit precision 00510 //------------------------- 00511 00512 inline half 00513 half::round (unsigned int n) const 00514 { 00515 // 00516 // Parameter check. 00517 // 00518 00519 if (n >= 10) 00520 return *this; 00521 00522 // 00523 // Disassemble h into the sign, s, 00524 // and the combined exponent and significand, e. 00525 // 00526 00527 unsigned short s = _h & 0x8000; 00528 unsigned short e = _h & 0x7fff; 00529 00530 // 00531 // Round the exponent and significand to the nearest value 00532 // where ones occur only in the (10-n) most significant bits. 00533 // Note that the exponent adjusts automatically if rounding 00534 // up causes the significand to overflow. 00535 // 00536 00537 e >>= 9 - n; 00538 e += e & 1; 00539 e <<= 9 - n; 00540 00541 // 00542 // Check for exponent overflow. 00543 // 00544 00545 if (e >= 0x7c00) 00546 { 00547 // 00548 // Overflow occurred -- truncate instead of rounding. 00549 // 00550 00551 e = _h; 00552 e >>= 10 - n; 00553 e <<= 10 - n; 00554 } 00555 00556 // 00557 // Put the original sign bit back. 00558 // 00559 00560 half h; 00561 h._h = s | e; 00562 00563 return h; 00564 } 00565 00566 00567 //----------------------- 00568 // Other inline functions 00569 //----------------------- 00570 00571 inline half 00572 half::operator - () const 00573 { 00574 half h; 00575 h._h = _h ^ 0x8000; 00576 return h; 00577 } 00578 00579 00580 inline half & 00581 half::operator = (half h) 00582 { 00583 _h = h._h; 00584 return *this; 00585 } 00586 00587 00588 inline half & 00589 half::operator = (float f) 00590 { 00591 *this = half (f); 00592 return *this; 00593 } 00594 00595 00596 inline half & 00597 half::operator += (half h) 00598 { 00599 *this = half (float (*this) + float (h)); 00600 return *this; 00601 } 00602 00603 00604 inline half & 00605 half::operator += (float f) 00606 { 00607 *this = half (float (*this) + f); 00608 return *this; 00609 } 00610 00611 00612 inline half & 00613 half::operator -= (half h) 00614 { 00615 *this = half (float (*this) - float (h)); 00616 return *this; 00617 } 00618 00619 00620 inline half & 00621 half::operator -= (float f) 00622 { 00623 *this = half (float (*this) - f); 00624 return *this; 00625 } 00626 00627 00628 inline half & 00629 half::operator *= (half h) 00630 { 00631 *this = half (float (*this) * float (h)); 00632 return *this; 00633 } 00634 00635 00636 inline half & 00637 half::operator *= (float f) 00638 { 00639 *this = half (float (*this) * f); 00640 return *this; 00641 } 00642 00643 00644 inline half & 00645 half::operator /= (half h) 00646 { 00647 *this = half (float (*this) / float (h)); 00648 return *this; 00649 } 00650 00651 00652 inline half & 00653 half::operator /= (float f) 00654 { 00655 *this = half (float (*this) / f); 00656 return *this; 00657 } 00658 00659 00660 inline bool 00661 half::isFinite () const 00662 { 00663 unsigned short e = (_h >> 10) & 0x001f; 00664 return e < 31; 00665 } 00666 00667 00668 inline bool 00669 half::isNormalized () const 00670 { 00671 unsigned short e = (_h >> 10) & 0x001f; 00672 return e > 0 && e < 31; 00673 } 00674 00675 00676 inline bool 00677 half::isDenormalized () const 00678 { 00679 unsigned short e = (_h >> 10) & 0x001f; 00680 unsigned short m = _h & 0x3ff; 00681 return e == 0 && m != 0; 00682 } 00683 00684 00685 inline bool 00686 half::isZero () const 00687 { 00688 return (_h & 0x7fff) == 0; 00689 } 00690 00691 00692 inline bool 00693 half::isNan () const 00694 { 00695 unsigned short e = (_h >> 10) & 0x001f; 00696 unsigned short m = _h & 0x3ff; 00697 return e == 31 && m != 0; 00698 } 00699 00700 00701 inline bool 00702 half::isInfinity () const 00703 { 00704 unsigned short e = (_h >> 10) & 0x001f; 00705 unsigned short m = _h & 0x3ff; 00706 return e == 31 && m == 0; 00707 } 00708 00709 00710 inline bool 00711 half::isNegative () const 00712 { 00713 return (_h & 0x8000) != 0; 00714 } 00715 00716 00717 inline half 00718 half::posInf () 00719 { 00720 half h; 00721 h._h = 0x7c00; 00722 return h; 00723 } 00724 00725 00726 inline half 00727 half::negInf () 00728 { 00729 half h; 00730 h._h = 0xfc00; 00731 return h; 00732 } 00733 00734 00735 inline half 00736 half::qNan () 00737 { 00738 half h; 00739 h._h = 0x7fff; 00740 return h; 00741 } 00742 00743 00744 inline half 00745 half::sNan () 00746 { 00747 half h; 00748 h._h = 0x7dff; 00749 return h; 00750 } 00751 00752 00753 inline unsigned short 00754 half::bits () const 00755 { 00756 return _h; 00757 } 00758 00759 00760 inline void 00761 half::setBits (unsigned short bits) 00762 { 00763 _h = bits; 00764 } 00765 00766 #endif