ieee754sp.c 4.2 KB

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