ieee754sp.c 4.6 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196
  1. /* IEEE754 floating point arithmetic
  2. * single precision
  3. */
  4. /*
  5. * MIPS floating point support
  6. * Copyright (C) 1994-2000 Algorithmics Ltd.
  7. *
  8. * This program is free software; you can distribute it and/or modify it
  9. * under the terms of the GNU General Public License (Version 2) as
  10. * published by the Free Software Foundation.
  11. *
  12. * This program is distributed in the hope it will be useful, but WITHOUT
  13. * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  14. * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
  15. * for more details.
  16. *
  17. * You should have received a copy of the GNU General Public License along
  18. * with this program; if not, write to the Free Software Foundation, Inc.,
  19. * 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA.
  20. */
  21. #include <linux/compiler.h>
  22. #include "ieee754sp.h"
  23. int ieee754sp_class(union ieee754sp x)
  24. {
  25. COMPXSP;
  26. EXPLODEXSP;
  27. return xc;
  28. }
  29. static inline int ieee754sp_isnan(union ieee754sp x)
  30. {
  31. return ieee754_class_nan(ieee754sp_class(x));
  32. }
  33. static inline int ieee754sp_issnan(union ieee754sp x)
  34. {
  35. assert(ieee754sp_isnan(x));
  36. return SPMANT(x) & SP_MBIT(SP_FBITS - 1);
  37. }
  38. /*
  39. * Raise the Invalid Operation IEEE 754 exception
  40. * and convert the signaling NaN supplied to a quiet NaN.
  41. */
  42. union ieee754sp __cold ieee754sp_nanxcpt(union ieee754sp r)
  43. {
  44. assert(ieee754sp_issnan(r));
  45. ieee754_setcx(IEEE754_INVALID_OPERATION);
  46. return ieee754sp_indef();
  47. }
  48. static unsigned ieee754sp_get_rounding(int sn, unsigned xm)
  49. {
  50. /* inexact must round of 3 bits
  51. */
  52. if (xm & (SP_MBIT(3) - 1)) {
  53. switch (ieee754_csr.rm) {
  54. case FPU_CSR_RZ:
  55. break;
  56. case FPU_CSR_RN:
  57. xm += 0x3 + ((xm >> 3) & 1);
  58. /* xm += (xm&0x8)?0x4:0x3 */
  59. break;
  60. case FPU_CSR_RU: /* toward +Infinity */
  61. if (!sn) /* ?? */
  62. xm += 0x8;
  63. break;
  64. case FPU_CSR_RD: /* toward -Infinity */
  65. if (sn) /* ?? */
  66. xm += 0x8;
  67. break;
  68. }
  69. }
  70. return xm;
  71. }
  72. /* generate a normal/denormal number with over,under handling
  73. * sn is sign
  74. * xe is an unbiased exponent
  75. * xm is 3bit extended precision value.
  76. */
  77. union ieee754sp ieee754sp_format(int sn, int xe, unsigned xm)
  78. {
  79. assert(xm); /* we don't gen exact zeros (probably should) */
  80. assert((xm >> (SP_FBITS + 1 + 3)) == 0); /* no execess */
  81. assert(xm & (SP_HIDDEN_BIT << 3));
  82. if (xe < SP_EMIN) {
  83. /* strip lower bits */
  84. int es = SP_EMIN - xe;
  85. if (ieee754_csr.nod) {
  86. ieee754_setcx(IEEE754_UNDERFLOW);
  87. ieee754_setcx(IEEE754_INEXACT);
  88. switch(ieee754_csr.rm) {
  89. case FPU_CSR_RN:
  90. case FPU_CSR_RZ:
  91. return ieee754sp_zero(sn);
  92. case FPU_CSR_RU: /* toward +Infinity */
  93. if (sn == 0)
  94. return ieee754sp_min(0);
  95. else
  96. return ieee754sp_zero(1);
  97. case FPU_CSR_RD: /* toward -Infinity */
  98. if (sn == 0)
  99. return ieee754sp_zero(0);
  100. else
  101. return ieee754sp_min(1);
  102. }
  103. }
  104. if (xe == SP_EMIN - 1 &&
  105. ieee754sp_get_rounding(sn, xm) >> (SP_FBITS + 1 + 3))
  106. {
  107. /* Not tiny after rounding */
  108. ieee754_setcx(IEEE754_INEXACT);
  109. xm = ieee754sp_get_rounding(sn, xm);
  110. xm >>= 1;
  111. /* Clear grs bits */
  112. xm &= ~(SP_MBIT(3) - 1);
  113. xe++;
  114. } else {
  115. /* sticky right shift es bits
  116. */
  117. SPXSRSXn(es);
  118. assert((xm & (SP_HIDDEN_BIT << 3)) == 0);
  119. assert(xe == SP_EMIN);
  120. }
  121. }
  122. if (xm & (SP_MBIT(3) - 1)) {
  123. ieee754_setcx(IEEE754_INEXACT);
  124. if ((xm & (SP_HIDDEN_BIT << 3)) == 0) {
  125. ieee754_setcx(IEEE754_UNDERFLOW);
  126. }
  127. /* inexact must round of 3 bits
  128. */
  129. xm = ieee754sp_get_rounding(sn, xm);
  130. /* adjust exponent for rounding add overflowing
  131. */
  132. if (xm >> (SP_FBITS + 1 + 3)) {
  133. /* add causes mantissa overflow */
  134. xm >>= 1;
  135. xe++;
  136. }
  137. }
  138. /* strip grs bits */
  139. xm >>= 3;
  140. assert((xm >> (SP_FBITS + 1)) == 0); /* no execess */
  141. assert(xe >= SP_EMIN);
  142. if (xe > SP_EMAX) {
  143. ieee754_setcx(IEEE754_OVERFLOW);
  144. ieee754_setcx(IEEE754_INEXACT);
  145. /* -O can be table indexed by (rm,sn) */
  146. switch (ieee754_csr.rm) {
  147. case FPU_CSR_RN:
  148. return ieee754sp_inf(sn);
  149. case FPU_CSR_RZ:
  150. return ieee754sp_max(sn);
  151. case FPU_CSR_RU: /* toward +Infinity */
  152. if (sn == 0)
  153. return ieee754sp_inf(0);
  154. else
  155. return ieee754sp_max(1);
  156. case FPU_CSR_RD: /* toward -Infinity */
  157. if (sn == 0)
  158. return ieee754sp_max(0);
  159. else
  160. return ieee754sp_inf(1);
  161. }
  162. }
  163. /* gen norm/denorm/zero */
  164. if ((xm & SP_HIDDEN_BIT) == 0) {
  165. /* we underflow (tiny/zero) */
  166. assert(xe == SP_EMIN);
  167. if (ieee754_csr.mx & IEEE754_UNDERFLOW)
  168. ieee754_setcx(IEEE754_UNDERFLOW);
  169. return buildsp(sn, SP_EMIN - 1 + SP_EBIAS, xm);
  170. } else {
  171. assert((xm >> (SP_FBITS + 1)) == 0); /* no execess */
  172. assert(xm & SP_HIDDEN_BIT);
  173. return buildsp(sn, xe + SP_EBIAS, xm & ~SP_HIDDEN_BIT);
  174. }
  175. }