sfmpy.c
上传用户:jlfgdled
上传日期:2013-04-10
资源大小:33168k
文件大小:10k
源码类别:

Linux/Unix编程

开发平台:

Unix_Linux

  1. /*
  2.  * Linux/PA-RISC Project (http://www.parisc-linux.org/)
  3.  *
  4.  * Floating-point emulation code
  5.  *  Copyright (C) 2001 Hewlett-Packard (Paul Bame) <bame@debian.org>
  6.  *
  7.  *    This program is free software; you can redistribute it and/or modify
  8.  *    it under the terms of the GNU General Public License as published by
  9.  *    the Free Software Foundation; either version 2, or (at your option)
  10.  *    any later version.
  11.  *
  12.  *    This program is distributed in the hope that it will be useful,
  13.  *    but WITHOUT ANY WARRANTY; without even the implied warranty of
  14.  *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  15.  *    GNU General Public License for more details.
  16.  *
  17.  *    You should have received a copy of the GNU General Public License
  18.  *    along with this program; if not, write to the Free Software
  19.  *    Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  20.  */
  21. /*
  22.  * BEGIN_DESC
  23.  *
  24.  *  File:
  25.  * @(#) pa/spmath/sfmpy.c $Revision: 1.1 $
  26.  *
  27.  *  Purpose:
  28.  * Single Precision Floating-point Multiply
  29.  *
  30.  *  External Interfaces:
  31.  * sgl_fmpy(srcptr1,srcptr2,dstptr,status)
  32.  *
  33.  *  Internal Interfaces:
  34.  *
  35.  *  Theory:
  36.  * <<please update with a overview of the operation of this file>>
  37.  *
  38.  * END_DESC
  39. */
  40. #include "float.h"
  41. #include "sgl_float.h"
  42. /*
  43.  *  Single Precision Floating-point Multiply
  44.  */
  45. int
  46. sgl_fmpy(
  47.     sgl_floating_point *srcptr1,
  48.     sgl_floating_point *srcptr2,
  49.     sgl_floating_point *dstptr,
  50.     unsigned int *status)
  51. {
  52. register unsigned int opnd1, opnd2, opnd3, result;
  53. register int dest_exponent, count;
  54. register boolean inexact = FALSE, guardbit = FALSE, stickybit = FALSE;
  55. boolean is_tiny;
  56. opnd1 = *srcptr1;
  57. opnd2 = *srcptr2;
  58. /* 
  59.  * set sign bit of result 
  60.  */
  61. if (Sgl_sign(opnd1) ^ Sgl_sign(opnd2)) Sgl_setnegativezero(result);  
  62. else Sgl_setzero(result);
  63. /*
  64.  * check first operand for NaN's or infinity
  65.  */
  66. if (Sgl_isinfinity_exponent(opnd1)) {
  67. if (Sgl_iszero_mantissa(opnd1)) {
  68. if (Sgl_isnotnan(opnd2)) {
  69. if (Sgl_iszero_exponentmantissa(opnd2)) {
  70. /* 
  71.  * invalid since operands are infinity 
  72.  * and zero 
  73.  */
  74. if (Is_invalidtrap_enabled()) 
  75.                                  return(INVALIDEXCEPTION);
  76.                                  Set_invalidflag();
  77.                                  Sgl_makequietnan(result);
  78. *dstptr = result;
  79. return(NOEXCEPTION);
  80. }
  81. /*
  82.    * return infinity
  83.    */
  84. Sgl_setinfinity_exponentmantissa(result);
  85. *dstptr = result;
  86. return(NOEXCEPTION);
  87. }
  88. }
  89. else {
  90.                  /*
  91.                    * is NaN; signaling or quiet?
  92.                    */
  93.                  if (Sgl_isone_signaling(opnd1)) {
  94.                          /* trap if INVALIDTRAP enabled */
  95.                          if (Is_invalidtrap_enabled()) 
  96.                              return(INVALIDEXCEPTION);
  97.                          /* make NaN quiet */
  98.                          Set_invalidflag();
  99.                          Sgl_set_quiet(opnd1);
  100.                  }
  101. /* 
  102.  * is second operand a signaling NaN? 
  103.  */
  104. else if (Sgl_is_signalingnan(opnd2)) {
  105.                          /* trap if INVALIDTRAP enabled */
  106.                          if (Is_invalidtrap_enabled()) 
  107.                              return(INVALIDEXCEPTION);
  108.                          /* make NaN quiet */
  109.                          Set_invalidflag();
  110.                          Sgl_set_quiet(opnd2);
  111.                  *dstptr = opnd2;
  112.                  return(NOEXCEPTION);
  113. }
  114.                  /*
  115.                    * return quiet NaN
  116.                    */
  117.                  *dstptr = opnd1;
  118.                  return(NOEXCEPTION);
  119. }
  120. }
  121. /*
  122.  * check second operand for NaN's or infinity
  123.  */
  124. if (Sgl_isinfinity_exponent(opnd2)) {
  125. if (Sgl_iszero_mantissa(opnd2)) {
  126. if (Sgl_iszero_exponentmantissa(opnd1)) {
  127. /* invalid since operands are zero & infinity */
  128. if (Is_invalidtrap_enabled()) 
  129.                                  return(INVALIDEXCEPTION);
  130.                                 Set_invalidflag();
  131.                                 Sgl_makequietnan(opnd2);
  132. *dstptr = opnd2;
  133. return(NOEXCEPTION);
  134. }
  135. /*
  136.  * return infinity
  137.  */
  138. Sgl_setinfinity_exponentmantissa(result);
  139. *dstptr = result;
  140. return(NOEXCEPTION);
  141. }
  142.                 /*
  143.                  * is NaN; signaling or quiet?
  144.                  */
  145.                 if (Sgl_isone_signaling(opnd2)) {
  146.                         /* trap if INVALIDTRAP enabled */
  147.                         if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
  148.                         /* make NaN quiet */
  149.                         Set_invalidflag();
  150.                         Sgl_set_quiet(opnd2);
  151.                 }
  152.                 /*
  153.                  * return quiet NaN
  154.                  */
  155.                 *dstptr = opnd2;
  156.                 return(NOEXCEPTION);
  157. }
  158. /*
  159.  * Generate exponent 
  160.  */
  161. dest_exponent = Sgl_exponent(opnd1) + Sgl_exponent(opnd2) - SGL_BIAS;
  162. /*
  163.  * Generate mantissa
  164.  */
  165. if (Sgl_isnotzero_exponent(opnd1)) {
  166. /* set hidden bit */
  167. Sgl_clear_signexponent_set_hidden(opnd1);
  168. }
  169. else {
  170. /* check for zero */
  171. if (Sgl_iszero_mantissa(opnd1)) {
  172. Sgl_setzero_exponentmantissa(result);
  173. *dstptr = result;
  174. return(NOEXCEPTION);
  175. }
  176.                 /* is denormalized, adjust exponent */
  177.                 Sgl_clear_signexponent(opnd1);
  178. Sgl_leftshiftby1(opnd1);
  179. Sgl_normalize(opnd1,dest_exponent);
  180. }
  181. /* opnd2 needs to have hidden bit set with msb in hidden bit */
  182. if (Sgl_isnotzero_exponent(opnd2)) {
  183. Sgl_clear_signexponent_set_hidden(opnd2);
  184. }
  185. else {
  186. /* check for zero */
  187. if (Sgl_iszero_mantissa(opnd2)) {
  188. Sgl_setzero_exponentmantissa(result);
  189. *dstptr = result;
  190. return(NOEXCEPTION);
  191. }
  192.                 /* is denormalized; want to normalize */
  193.                 Sgl_clear_signexponent(opnd2);
  194.                 Sgl_leftshiftby1(opnd2);
  195. Sgl_normalize(opnd2,dest_exponent);
  196. }
  197. /* Multiply two source mantissas together */
  198. Sgl_leftshiftby4(opnd2);     /* make room for guard bits */
  199. Sgl_setzero(opnd3);
  200. /*
  201.  * Four bits at a time are inspected in each loop, and a
  202.  * simple shift and add multiply algorithm is used.
  203.  */
  204. for (count=1;count<SGL_P;count+=4) {
  205. stickybit |= Slow4(opnd3);
  206. Sgl_rightshiftby4(opnd3);
  207. if (Sbit28(opnd1)) Sall(opnd3) += (Sall(opnd2) << 3);
  208. if (Sbit29(opnd1)) Sall(opnd3) += (Sall(opnd2) << 2);
  209. if (Sbit30(opnd1)) Sall(opnd3) += (Sall(opnd2) << 1);
  210. if (Sbit31(opnd1)) Sall(opnd3) += Sall(opnd2);
  211. Sgl_rightshiftby4(opnd1);
  212. }
  213. /* make sure result is left-justified */
  214. if (Sgl_iszero_sign(opnd3)) {
  215. Sgl_leftshiftby1(opnd3);
  216. }
  217. else {
  218. /* result mantissa >= 2. */
  219. dest_exponent++;
  220. }
  221. /* check for denormalized result */
  222. while (Sgl_iszero_sign(opnd3)) {
  223. Sgl_leftshiftby1(opnd3);
  224. dest_exponent--;
  225. }
  226. /*
  227.  * check for guard, sticky and inexact bits
  228.  */
  229. stickybit |= Sgl_all(opnd3) << (SGL_BITLENGTH - SGL_EXP_LENGTH + 1);
  230. guardbit = Sbit24(opnd3);
  231. inexact = guardbit | stickybit;
  232. /* re-align mantissa */
  233. Sgl_rightshiftby8(opnd3);
  234. /* 
  235.  * round result 
  236.  */
  237. if (inexact && (dest_exponent>0 || Is_underflowtrap_enabled())) {
  238. Sgl_clear_signexponent(opnd3);
  239. switch (Rounding_mode()) {
  240. case ROUNDPLUS: 
  241. if (Sgl_iszero_sign(result)) 
  242. Sgl_increment(opnd3);
  243. break;
  244. case ROUNDMINUS: 
  245. if (Sgl_isone_sign(result)) 
  246. Sgl_increment(opnd3);
  247. break;
  248. case ROUNDNEAREST:
  249. if (guardbit) {
  250.     if (stickybit || Sgl_isone_lowmantissa(opnd3))
  251.        Sgl_increment(opnd3);
  252. }
  253. }
  254. if (Sgl_isone_hidden(opnd3)) dest_exponent++;
  255. }
  256. Sgl_set_mantissa(result,opnd3);
  257.         /* 
  258.          * Test for overflow
  259.          */
  260. if (dest_exponent >= SGL_INFINITY_EXPONENT) {
  261.                 /* trap if OVERFLOWTRAP enabled */
  262.                 if (Is_overflowtrap_enabled()) {
  263.                         /*
  264.                          * Adjust bias of result
  265.                          */
  266. Sgl_setwrapped_exponent(result,dest_exponent,ovfl);
  267. *dstptr = result;
  268. if (inexact) 
  269.     if (Is_inexacttrap_enabled())
  270. return(OVERFLOWEXCEPTION | INEXACTEXCEPTION);
  271.     else Set_inexactflag();
  272. return(OVERFLOWEXCEPTION);
  273.                 }
  274. inexact = TRUE;
  275. Set_overflowflag();
  276.                 /* set result to infinity or largest number */
  277. Sgl_setoverflow(result);
  278. }
  279.         /* 
  280.          * Test for underflow
  281.          */
  282. else if (dest_exponent <= 0) {
  283.                 /* trap if UNDERFLOWTRAP enabled */
  284.                 if (Is_underflowtrap_enabled()) {
  285.                         /*
  286.                          * Adjust bias of result
  287.                          */
  288. Sgl_setwrapped_exponent(result,dest_exponent,unfl);
  289. *dstptr = result;
  290. if (inexact) 
  291.     if (Is_inexacttrap_enabled())
  292. return(UNDERFLOWEXCEPTION | INEXACTEXCEPTION);
  293.     else Set_inexactflag();
  294. return(UNDERFLOWEXCEPTION);
  295.                 }
  296. /* Determine if should set underflow flag */
  297. is_tiny = TRUE;
  298. if (dest_exponent == 0 && inexact) {
  299. switch (Rounding_mode()) {
  300. case ROUNDPLUS: 
  301. if (Sgl_iszero_sign(result)) {
  302. Sgl_increment(opnd3);
  303. if (Sgl_isone_hiddenoverflow(opnd3))
  304.                      is_tiny = FALSE;
  305. Sgl_decrement(opnd3);
  306. }
  307. break;
  308. case ROUNDMINUS: 
  309. if (Sgl_isone_sign(result)) {
  310. Sgl_increment(opnd3);
  311. if (Sgl_isone_hiddenoverflow(opnd3))
  312.                      is_tiny = FALSE;
  313. Sgl_decrement(opnd3);
  314. }
  315. break;
  316. case ROUNDNEAREST:
  317. if (guardbit && (stickybit || 
  318.     Sgl_isone_lowmantissa(opnd3))) {
  319.        Sgl_increment(opnd3);
  320. if (Sgl_isone_hiddenoverflow(opnd3))
  321.                      is_tiny = FALSE;
  322. Sgl_decrement(opnd3);
  323. }
  324. break;
  325. }
  326. }
  327.                 /*
  328.                  * denormalize result or set to signed zero
  329.                  */
  330. stickybit = inexact;
  331. Sgl_denormalize(opnd3,dest_exponent,guardbit,stickybit,inexact);
  332. /* return zero or smallest number */
  333. if (inexact) {
  334. switch (Rounding_mode()) {
  335. case ROUNDPLUS: 
  336. if (Sgl_iszero_sign(result)) {
  337. Sgl_increment(opnd3);
  338. }
  339. break;
  340. case ROUNDMINUS: 
  341. if (Sgl_isone_sign(result)) {
  342. Sgl_increment(opnd3);
  343. }
  344. break;
  345. case ROUNDNEAREST:
  346. if (guardbit && (stickybit || 
  347.     Sgl_isone_lowmantissa(opnd3))) {
  348.        Sgl_increment(opnd3);
  349. }
  350. break;
  351. }
  352.                 if (is_tiny) Set_underflowflag();
  353. }
  354. Sgl_set_exponentmantissa(result,opnd3);
  355. }
  356. else Sgl_set_exponent(result,dest_exponent);
  357. *dstptr = result;
  358. /* check for inexact */
  359. if (inexact) {
  360. if (Is_inexacttrap_enabled()) return(INEXACTEXCEPTION);
  361. else Set_inexactflag();
  362. }
  363. return(NOEXCEPTION);
  364. }