ieee754sp.c
上传用户:lgb322
上传日期:2013-02-24
资源大小:30529k
文件大小:4k
源码类别:

嵌入式Linux

开发平台:

Unix_Linux

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