ieee754sp.c 5.2 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244
  1. /* IEEE754 floating point arithmetic
  2. * single precision
  3. */
  4. /*
  5. * MIPS floating point support
  6. * Copyright (C) 1994-2000 Algorithmics Ltd.
  7. *
  8. * ########################################################################
  9. *
  10. * This program is free software; you can distribute it and/or modify it
  11. * under the terms of the GNU General Public License (Version 2) as
  12. * published by the Free Software Foundation.
  13. *
  14. * This program is distributed in the hope it will be useful, but WITHOUT
  15. * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  16. * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
  17. * for more details.
  18. *
  19. * You should have received a copy of the GNU General Public License along
  20. * with this program; if not, write to the Free Software Foundation, Inc.,
  21. * 59 Temple Place - Suite 330, Boston MA 02111-1307, USA.
  22. *
  23. * ########################################################################
  24. */
  25. #include "ieee754sp.h"
  26. int ieee754sp_class(ieee754sp x)
  27. {
  28. COMPXSP;
  29. EXPLODEXSP;
  30. return xc;
  31. }
  32. int ieee754sp_isnan(ieee754sp x)
  33. {
  34. return ieee754sp_class(x) >= IEEE754_CLASS_SNAN;
  35. }
  36. int ieee754sp_issnan(ieee754sp x)
  37. {
  38. assert(ieee754sp_isnan(x));
  39. return (SPMANT(x) & SP_MBIT(SP_MBITS-1));
  40. }
  41. ieee754sp ieee754sp_xcpt(ieee754sp r, const char *op, ...)
  42. {
  43. struct ieee754xctx ax;
  44. if (!TSTX())
  45. return r;
  46. ax.op = op;
  47. ax.rt = IEEE754_RT_SP;
  48. ax.rv.sp = r;
  49. va_start(ax.ap, op);
  50. ieee754_xcpt(&ax);
  51. va_end(ax.ap);
  52. return ax.rv.sp;
  53. }
  54. ieee754sp ieee754sp_nanxcpt(ieee754sp r, const char *op, ...)
  55. {
  56. struct ieee754xctx ax;
  57. assert(ieee754sp_isnan(r));
  58. if (!ieee754sp_issnan(r)) /* QNAN does not cause invalid op !! */
  59. return r;
  60. if (!SETANDTESTCX(IEEE754_INVALID_OPERATION)) {
  61. /* not enabled convert to a quiet NaN */
  62. SPMANT(r) &= (~SP_MBIT(SP_MBITS-1));
  63. if (ieee754sp_isnan(r))
  64. return r;
  65. else
  66. return ieee754sp_indef();
  67. }
  68. ax.op = op;
  69. ax.rt = 0;
  70. ax.rv.sp = r;
  71. va_start(ax.ap, op);
  72. ieee754_xcpt(&ax);
  73. va_end(ax.ap);
  74. return ax.rv.sp;
  75. }
  76. ieee754sp ieee754sp_bestnan(ieee754sp x, ieee754sp y)
  77. {
  78. assert(ieee754sp_isnan(x));
  79. assert(ieee754sp_isnan(y));
  80. if (SPMANT(x) > SPMANT(y))
  81. return x;
  82. else
  83. return y;
  84. }
  85. static unsigned get_rounding(int sn, unsigned xm)
  86. {
  87. /* inexact must round of 3 bits
  88. */
  89. if (xm & (SP_MBIT(3) - 1)) {
  90. switch (ieee754_csr.rm) {
  91. case IEEE754_RZ:
  92. break;
  93. case IEEE754_RN:
  94. xm += 0x3 + ((xm >> 3) & 1);
  95. /* xm += (xm&0x8)?0x4:0x3 */
  96. break;
  97. case IEEE754_RU: /* toward +Infinity */
  98. if (!sn) /* ?? */
  99. xm += 0x8;
  100. break;
  101. case IEEE754_RD: /* toward -Infinity */
  102. if (sn) /* ?? */
  103. xm += 0x8;
  104. break;
  105. }
  106. }
  107. return xm;
  108. }
  109. /* generate a normal/denormal number with over,under handling
  110. * sn is sign
  111. * xe is an unbiased exponent
  112. * xm is 3bit extended precision value.
  113. */
  114. ieee754sp ieee754sp_format(int sn, int xe, unsigned xm)
  115. {
  116. assert(xm); /* we don't gen exact zeros (probably should) */
  117. assert((xm >> (SP_MBITS + 1 + 3)) == 0); /* no execess */
  118. assert(xm & (SP_HIDDEN_BIT << 3));
  119. if (xe < SP_EMIN) {
  120. /* strip lower bits */
  121. int es = SP_EMIN - xe;
  122. if (ieee754_csr.nod) {
  123. SETCX(IEEE754_UNDERFLOW);
  124. SETCX(IEEE754_INEXACT);
  125. switch(ieee754_csr.rm) {
  126. case IEEE754_RN:
  127. case IEEE754_RZ:
  128. return ieee754sp_zero(sn);
  129. case IEEE754_RU: /* toward +Infinity */
  130. if(sn == 0)
  131. return ieee754sp_min(0);
  132. else
  133. return ieee754sp_zero(1);
  134. case IEEE754_RD: /* toward -Infinity */
  135. if(sn == 0)
  136. return ieee754sp_zero(0);
  137. else
  138. return ieee754sp_min(1);
  139. }
  140. }
  141. if (xe == SP_EMIN - 1
  142. && get_rounding(sn, xm) >> (SP_MBITS + 1 + 3))
  143. {
  144. /* Not tiny after rounding */
  145. SETCX(IEEE754_INEXACT);
  146. xm = get_rounding(sn, xm);
  147. xm >>= 1;
  148. /* Clear grs bits */
  149. xm &= ~(SP_MBIT(3) - 1);
  150. xe++;
  151. }
  152. else {
  153. /* sticky right shift es bits
  154. */
  155. SPXSRSXn(es);
  156. assert((xm & (SP_HIDDEN_BIT << 3)) == 0);
  157. assert(xe == SP_EMIN);
  158. }
  159. }
  160. if (xm & (SP_MBIT(3) - 1)) {
  161. SETCX(IEEE754_INEXACT);
  162. if ((xm & (SP_HIDDEN_BIT << 3)) == 0) {
  163. SETCX(IEEE754_UNDERFLOW);
  164. }
  165. /* inexact must round of 3 bits
  166. */
  167. xm = get_rounding(sn, xm);
  168. /* adjust exponent for rounding add overflowing
  169. */
  170. if (xm >> (SP_MBITS + 1 + 3)) {
  171. /* add causes mantissa overflow */
  172. xm >>= 1;
  173. xe++;
  174. }
  175. }
  176. /* strip grs bits */
  177. xm >>= 3;
  178. assert((xm >> (SP_MBITS + 1)) == 0); /* no execess */
  179. assert(xe >= SP_EMIN);
  180. if (xe > SP_EMAX) {
  181. SETCX(IEEE754_OVERFLOW);
  182. SETCX(IEEE754_INEXACT);
  183. /* -O can be table indexed by (rm,sn) */
  184. switch (ieee754_csr.rm) {
  185. case IEEE754_RN:
  186. return ieee754sp_inf(sn);
  187. case IEEE754_RZ:
  188. return ieee754sp_max(sn);
  189. case IEEE754_RU: /* toward +Infinity */
  190. if (sn == 0)
  191. return ieee754sp_inf(0);
  192. else
  193. return ieee754sp_max(1);
  194. case IEEE754_RD: /* toward -Infinity */
  195. if (sn == 0)
  196. return ieee754sp_max(0);
  197. else
  198. return ieee754sp_inf(1);
  199. }
  200. }
  201. /* gen norm/denorm/zero */
  202. if ((xm & SP_HIDDEN_BIT) == 0) {
  203. /* we underflow (tiny/zero) */
  204. assert(xe == SP_EMIN);
  205. if (ieee754_csr.mx & IEEE754_UNDERFLOW)
  206. SETCX(IEEE754_UNDERFLOW);
  207. return buildsp(sn, SP_EMIN - 1 + SP_EBIAS, xm);
  208. } else {
  209. assert((xm >> (SP_MBITS + 1)) == 0); /* no execess */
  210. assert(xm & SP_HIDDEN_BIT);
  211. return buildsp(sn, xe + SP_EBIAS, xm & ~SP_HIDDEN_BIT);
  212. }
  213. }