ieee754dp.c 4.2 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197
  1. // SPDX-License-Identifier: GPL-2.0-only
  2. /* IEEE754 floating point arithmetic
  3. * double precision: common utilities
  4. */
  5. /*
  6. * MIPS floating point support
  7. * Copyright (C) 1994-2000 Algorithmics Ltd.
  8. */
  9. #include <linux/compiler.h>
  10. #include "ieee754dp.h"
  11. int ieee754dp_class(union ieee754dp x)
  12. {
  13. COMPXDP;
  14. EXPLODEXDP;
  15. return xc;
  16. }
  17. static inline int ieee754dp_isnan(union ieee754dp x)
  18. {
  19. return ieee754_class_nan(ieee754dp_class(x));
  20. }
  21. static inline int ieee754dp_issnan(union ieee754dp x)
  22. {
  23. int qbit;
  24. assert(ieee754dp_isnan(x));
  25. qbit = (DPMANT(x) & DP_MBIT(DP_FBITS - 1)) == DP_MBIT(DP_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 ieee754dp __cold ieee754dp_nanxcpt(union ieee754dp r)
  33. {
  34. assert(ieee754dp_issnan(r));
  35. ieee754_setcx(IEEE754_INVALID_OPERATION);
  36. if (ieee754_csr.nan2008) {
  37. DPMANT(r) |= DP_MBIT(DP_FBITS - 1);
  38. } else {
  39. DPMANT(r) &= ~DP_MBIT(DP_FBITS - 1);
  40. if (!ieee754dp_isnan(r))
  41. DPMANT(r) |= DP_MBIT(DP_FBITS - 2);
  42. }
  43. return r;
  44. }
  45. static u64 ieee754dp_get_rounding(int sn, u64 xm)
  46. {
  47. /* inexact must round of 3 bits
  48. */
  49. if (xm & (DP_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 ieee754dp ieee754dp_format(int sn, int xe, u64 xm)
  75. {
  76. assert(xm); /* we don't gen exact zeros (probably should) */
  77. assert((xm >> (DP_FBITS + 1 + 3)) == 0); /* no excess */
  78. assert(xm & (DP_HIDDEN_BIT << 3));
  79. if (xe < DP_EMIN) {
  80. /* strip lower bits */
  81. int es = DP_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 ieee754dp_zero(sn);
  89. case FPU_CSR_RU: /* toward +Infinity */
  90. if (sn == 0)
  91. return ieee754dp_min(0);
  92. else
  93. return ieee754dp_zero(1);
  94. case FPU_CSR_RD: /* toward -Infinity */
  95. if (sn == 0)
  96. return ieee754dp_zero(0);
  97. else
  98. return ieee754dp_min(1);
  99. }
  100. }
  101. if (xe == DP_EMIN - 1 &&
  102. ieee754dp_get_rounding(sn, xm) >> (DP_FBITS + 1 + 3))
  103. {
  104. /* Not tiny after rounding */
  105. ieee754_setcx(IEEE754_INEXACT);
  106. xm = ieee754dp_get_rounding(sn, xm);
  107. xm >>= 1;
  108. /* Clear grs bits */
  109. xm &= ~(DP_MBIT(3) - 1);
  110. xe++;
  111. }
  112. else {
  113. /* sticky right shift es bits
  114. */
  115. xm = XDPSRS(xm, es);
  116. xe += es;
  117. assert((xm & (DP_HIDDEN_BIT << 3)) == 0);
  118. assert(xe == DP_EMIN);
  119. }
  120. }
  121. if (xm & (DP_MBIT(3) - 1)) {
  122. ieee754_setcx(IEEE754_INEXACT);
  123. if ((xm & (DP_HIDDEN_BIT << 3)) == 0) {
  124. ieee754_setcx(IEEE754_UNDERFLOW);
  125. }
  126. /* inexact must round of 3 bits
  127. */
  128. xm = ieee754dp_get_rounding(sn, xm);
  129. /* adjust exponent for rounding add overflowing
  130. */
  131. if (xm >> (DP_FBITS + 3 + 1)) {
  132. /* add causes mantissa overflow */
  133. xm >>= 1;
  134. xe++;
  135. }
  136. }
  137. /* strip grs bits */
  138. xm >>= 3;
  139. assert((xm >> (DP_FBITS + 1)) == 0); /* no excess */
  140. assert(xe >= DP_EMIN);
  141. if (xe > DP_EMAX) {
  142. ieee754_setcx(IEEE754_OVERFLOW);
  143. ieee754_setcx(IEEE754_INEXACT);
  144. /* -O can be table indexed by (rm,sn) */
  145. switch (ieee754_csr.rm) {
  146. case FPU_CSR_RN:
  147. return ieee754dp_inf(sn);
  148. case FPU_CSR_RZ:
  149. return ieee754dp_max(sn);
  150. case FPU_CSR_RU: /* toward +Infinity */
  151. if (sn == 0)
  152. return ieee754dp_inf(0);
  153. else
  154. return ieee754dp_max(1);
  155. case FPU_CSR_RD: /* toward -Infinity */
  156. if (sn == 0)
  157. return ieee754dp_max(0);
  158. else
  159. return ieee754dp_inf(1);
  160. }
  161. }
  162. /* gen norm/denorm/zero */
  163. if ((xm & DP_HIDDEN_BIT) == 0) {
  164. /* we underflow (tiny/zero) */
  165. assert(xe == DP_EMIN);
  166. if (ieee754_csr.mx & IEEE754_UNDERFLOW)
  167. ieee754_setcx(IEEE754_UNDERFLOW);
  168. return builddp(sn, DP_EMIN - 1 + DP_EBIAS, xm);
  169. } else {
  170. assert((xm >> (DP_FBITS + 1)) == 0); /* no excess */
  171. assert(xm & DP_HIDDEN_BIT);
  172. return builddp(sn, xe + DP_EBIAS, xm & ~DP_HIDDEN_BIT);
  173. }
  174. }