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

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/sfrem.c $Revision: 1.1 $
  26.  *
  27.  *  Purpose:
  28.  * Single Precision Floating-point Remainder
  29.  *
  30.  *  External Interfaces:
  31.  * sgl_frem(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 Remainder
  44.  */
  45. int
  46. sgl_frem (sgl_floating_point * srcptr1, sgl_floating_point * srcptr2,
  47.   sgl_floating_point * dstptr, unsigned int *status)
  48. {
  49. register unsigned int opnd1, opnd2, result;
  50. register int opnd1_exponent, opnd2_exponent, dest_exponent, stepcount;
  51. register boolean roundup = FALSE;
  52. opnd1 = *srcptr1;
  53. opnd2 = *srcptr2;
  54. /*
  55.  * check first operand for NaN's or infinity
  56.  */
  57. if ((opnd1_exponent = Sgl_exponent(opnd1)) == SGL_INFINITY_EXPONENT) {
  58. if (Sgl_iszero_mantissa(opnd1)) {
  59. if (Sgl_isnotnan(opnd2)) {
  60. /* invalid since first operand is infinity */
  61. if (Is_invalidtrap_enabled()) 
  62.                                  return(INVALIDEXCEPTION);
  63.                                 Set_invalidflag();
  64.                                 Sgl_makequietnan(result);
  65. *dstptr = result;
  66. return(NOEXCEPTION);
  67. }
  68. }
  69. else {
  70.                  /*
  71.                    * is NaN; signaling or quiet?
  72.                    */
  73.                  if (Sgl_isone_signaling(opnd1)) {
  74.                          /* trap if INVALIDTRAP enabled */
  75.                          if (Is_invalidtrap_enabled()) 
  76.                              return(INVALIDEXCEPTION);
  77.                          /* make NaN quiet */
  78.                          Set_invalidflag();
  79.                          Sgl_set_quiet(opnd1);
  80.                  }
  81. /* 
  82.  * is second operand a signaling NaN? 
  83.  */
  84. else if (Sgl_is_signalingnan(opnd2)) {
  85.                          /* trap if INVALIDTRAP enabled */
  86.                          if (Is_invalidtrap_enabled()) 
  87.                              return(INVALIDEXCEPTION);
  88.                          /* make NaN quiet */
  89.                          Set_invalidflag();
  90.                          Sgl_set_quiet(opnd2);
  91.                  *dstptr = opnd2;
  92.                  return(NOEXCEPTION);
  93. }
  94.                  /*
  95.                    * return quiet NaN
  96.                    */
  97.                  *dstptr = opnd1;
  98.                  return(NOEXCEPTION);
  99. }
  100. /*
  101.  * check second operand for NaN's or infinity
  102.  */
  103. if ((opnd2_exponent = Sgl_exponent(opnd2)) == SGL_INFINITY_EXPONENT) {
  104. if (Sgl_iszero_mantissa(opnd2)) {
  105. /*
  106.  * return first operand
  107.  */
  108.                  *dstptr = opnd1;
  109. return(NOEXCEPTION);
  110. }
  111.                 /*
  112.                  * is NaN; signaling or quiet?
  113.                  */
  114.                 if (Sgl_isone_signaling(opnd2)) {
  115.                         /* trap if INVALIDTRAP enabled */
  116.                         if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
  117.                         /* make NaN quiet */
  118.                         Set_invalidflag();
  119.                         Sgl_set_quiet(opnd2);
  120.                 }
  121.                 /*
  122.                  * return quiet NaN
  123.                  */
  124.                 *dstptr = opnd2;
  125.                 return(NOEXCEPTION);
  126. }
  127. /*
  128.  * check second operand for zero
  129.  */
  130. if (Sgl_iszero_exponentmantissa(opnd2)) {
  131. /* invalid since second operand is zero */
  132. if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
  133.                 Set_invalidflag();
  134.                 Sgl_makequietnan(result);
  135. *dstptr = result;
  136. return(NOEXCEPTION);
  137. }
  138. /* 
  139.  * get sign of result
  140.  */
  141. result = opnd1;  
  142. /* 
  143.  * check for denormalized operands
  144.  */
  145. if (opnd1_exponent == 0) {
  146. /* check for zero */
  147. if (Sgl_iszero_mantissa(opnd1)) {
  148. *dstptr = opnd1;
  149. return(NOEXCEPTION);
  150. }
  151. /* normalize, then continue */
  152. opnd1_exponent = 1;
  153. Sgl_normalize(opnd1,opnd1_exponent);
  154. }
  155. else {
  156. Sgl_clear_signexponent_set_hidden(opnd1);
  157. }
  158. if (opnd2_exponent == 0) {
  159. /* normalize, then continue */
  160. opnd2_exponent = 1;
  161. Sgl_normalize(opnd2,opnd2_exponent);
  162. }
  163. else {
  164. Sgl_clear_signexponent_set_hidden(opnd2);
  165. }
  166. /* find result exponent and divide step loop count */
  167. dest_exponent = opnd2_exponent - 1;
  168. stepcount = opnd1_exponent - opnd2_exponent;
  169. /*
  170.  * check for opnd1/opnd2 < 1
  171.  */
  172. if (stepcount < 0) {
  173. /*
  174.  * check for opnd1/opnd2 > 1/2
  175.  *
  176.  * In this case n will round to 1, so 
  177.  *    r = opnd1 - opnd2 
  178.  */
  179. if (stepcount == -1 && Sgl_isgreaterthan(opnd1,opnd2)) {
  180. Sgl_all(result) = ~Sgl_all(result);   /* set sign */
  181. /* align opnd2 with opnd1 */
  182. Sgl_leftshiftby1(opnd2); 
  183. Sgl_subtract(opnd2,opnd1,opnd2);
  184. /* now normalize */
  185.                  while (Sgl_iszero_hidden(opnd2)) {
  186.                          Sgl_leftshiftby1(opnd2);
  187.                          dest_exponent--;
  188. }
  189. Sgl_set_exponentmantissa(result,opnd2);
  190. goto testforunderflow;
  191. }
  192. /*
  193.  * opnd1/opnd2 <= 1/2
  194.  *
  195.  * In this case n will round to zero, so 
  196.  *    r = opnd1
  197.  */
  198. Sgl_set_exponentmantissa(result,opnd1);
  199. dest_exponent = opnd1_exponent;
  200. goto testforunderflow;
  201. }
  202. /*
  203.  * Generate result
  204.  *
  205.  * Do iterative subtract until remainder is less than operand 2.
  206.  */
  207. while (stepcount-- > 0 && Sgl_all(opnd1)) {
  208. if (Sgl_isnotlessthan(opnd1,opnd2))
  209. Sgl_subtract(opnd1,opnd2,opnd1);
  210. Sgl_leftshiftby1(opnd1);
  211. }
  212. /*
  213.  * Do last subtract, then determine which way to round if remainder 
  214.  * is exactly 1/2 of opnd2 
  215.  */
  216. if (Sgl_isnotlessthan(opnd1,opnd2)) {
  217. Sgl_subtract(opnd1,opnd2,opnd1);
  218. roundup = TRUE;
  219. }
  220. if (stepcount > 0 || Sgl_iszero(opnd1)) {
  221. /* division is exact, remainder is zero */
  222. Sgl_setzero_exponentmantissa(result);
  223. *dstptr = result;
  224. return(NOEXCEPTION);
  225. }
  226. /* 
  227.  * Check for cases where opnd1/opnd2 < n 
  228.  *
  229.  * In this case the result's sign will be opposite that of
  230.  * opnd1.  The mantissa also needs some correction.
  231.  */
  232. Sgl_leftshiftby1(opnd1);
  233. if (Sgl_isgreaterthan(opnd1,opnd2)) {
  234. Sgl_invert_sign(result);
  235. Sgl_subtract((opnd2<<1),opnd1,opnd1);
  236. }
  237. /* check for remainder being exactly 1/2 of opnd2 */
  238. else if (Sgl_isequal(opnd1,opnd2) && roundup) { 
  239. Sgl_invert_sign(result);
  240. }
  241. /* normalize result's mantissa */
  242.         while (Sgl_iszero_hidden(opnd1)) {
  243.                 dest_exponent--;
  244.                 Sgl_leftshiftby1(opnd1);
  245.         }
  246. Sgl_set_exponentmantissa(result,opnd1);
  247.         /* 
  248.          * Test for underflow
  249.          */
  250.     testforunderflow:
  251. if (dest_exponent <= 0) {
  252.                 /* trap if UNDERFLOWTRAP enabled */
  253.                 if (Is_underflowtrap_enabled()) {
  254.                         /*
  255.                          * Adjust bias of result
  256.                          */
  257.                         Sgl_setwrapped_exponent(result,dest_exponent,unfl);
  258. *dstptr = result;
  259. /* frem is always exact */
  260. return(UNDERFLOWEXCEPTION);
  261.                 }
  262.                 /*
  263.                  * denormalize result or set to signed zero
  264.                  */
  265.                 if (dest_exponent >= (1 - SGL_P)) {
  266. Sgl_rightshift_exponentmantissa(result,1-dest_exponent);
  267.                 }
  268.                 else {
  269. Sgl_setzero_exponentmantissa(result);
  270. }
  271. }
  272. else Sgl_set_exponent(result,dest_exponent);
  273. *dstptr = result;
  274. return(NOEXCEPTION);
  275. }