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

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/sfdiv.c $Revision: 1.1 $
  26.  *
  27.  *  Purpose:
  28.  * Single Precision Floating-point Divide
  29.  *
  30.  *  External Interfaces:
  31.  * sgl_fdiv(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 Divide
  44.  */
  45. int
  46. sgl_fdiv (sgl_floating_point * srcptr1, sgl_floating_point * srcptr2,
  47.   sgl_floating_point * dstptr, unsigned int *status)
  48. {
  49. register unsigned int opnd1, opnd2, opnd3, result;
  50. register int dest_exponent, count;
  51. register boolean inexact = FALSE, guardbit = FALSE, stickybit = FALSE;
  52. boolean is_tiny;
  53. opnd1 = *srcptr1;
  54. opnd2 = *srcptr2;
  55. /* 
  56.  * set sign bit of result 
  57.  */
  58. if (Sgl_sign(opnd1) ^ Sgl_sign(opnd2)) Sgl_setnegativezero(result);  
  59. else Sgl_setzero(result);
  60. /*
  61.  * check first operand for NaN's or infinity
  62.  */
  63. if (Sgl_isinfinity_exponent(opnd1)) {
  64. if (Sgl_iszero_mantissa(opnd1)) {
  65. if (Sgl_isnotnan(opnd2)) {
  66. if (Sgl_isinfinity(opnd2)) {
  67. /* 
  68.  * invalid since both operands 
  69.  * are infinity 
  70.  */
  71. if (Is_invalidtrap_enabled()) 
  72.                                  return(INVALIDEXCEPTION);
  73.                                  Set_invalidflag();
  74.                                  Sgl_makequietnan(result);
  75. *dstptr = result;
  76. return(NOEXCEPTION);
  77. }
  78. /*
  79.    * return infinity
  80.    */
  81. Sgl_setinfinity_exponentmantissa(result);
  82. *dstptr = result;
  83. return(NOEXCEPTION);
  84. }
  85. }
  86. else {
  87.                  /*
  88.                    * is NaN; signaling or quiet?
  89.                    */
  90.                  if (Sgl_isone_signaling(opnd1)) {
  91.                          /* trap if INVALIDTRAP enabled */
  92.                          if (Is_invalidtrap_enabled()) 
  93.                              return(INVALIDEXCEPTION);
  94.                          /* make NaN quiet */
  95.                          Set_invalidflag();
  96.                          Sgl_set_quiet(opnd1);
  97.                  }
  98. /* 
  99.  * is second operand a signaling NaN? 
  100.  */
  101. else if (Sgl_is_signalingnan(opnd2)) {
  102.                          /* trap if INVALIDTRAP enabled */
  103.                          if (Is_invalidtrap_enabled())
  104.                              return(INVALIDEXCEPTION);
  105.                          /* make NaN quiet */
  106.                          Set_invalidflag();
  107.                          Sgl_set_quiet(opnd2);
  108.                  *dstptr = opnd2;
  109.                  return(NOEXCEPTION);
  110. }
  111.                  /*
  112.                    * return quiet NaN
  113.                    */
  114.                  *dstptr = opnd1;
  115.                  return(NOEXCEPTION);
  116. }
  117. }
  118. /*
  119.  * check second operand for NaN's or infinity
  120.  */
  121. if (Sgl_isinfinity_exponent(opnd2)) {
  122. if (Sgl_iszero_mantissa(opnd2)) {
  123. /*
  124.  * return zero
  125.  */
  126. Sgl_setzero_exponentmantissa(result);
  127. *dstptr = result;
  128. return(NOEXCEPTION);
  129. }
  130.                 /*
  131.                  * is NaN; signaling or quiet?
  132.                  */
  133.                 if (Sgl_isone_signaling(opnd2)) {
  134.                         /* trap if INVALIDTRAP enabled */
  135.                         if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
  136.                         /* make NaN quiet */
  137.                         Set_invalidflag();
  138.                         Sgl_set_quiet(opnd2);
  139.                 }
  140.                 /*
  141.                  * return quiet NaN
  142.                  */
  143.                 *dstptr = opnd2;
  144.                 return(NOEXCEPTION);
  145. }
  146. /*
  147.  * check for division by zero
  148.  */
  149. if (Sgl_iszero_exponentmantissa(opnd2)) {
  150. if (Sgl_iszero_exponentmantissa(opnd1)) {
  151. /* invalid since both operands are zero */
  152. if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
  153.                         Set_invalidflag();
  154.                         Sgl_makequietnan(result);
  155. *dstptr = result;
  156. return(NOEXCEPTION);
  157. }
  158. if (Is_divisionbyzerotrap_enabled())
  159.                         return(DIVISIONBYZEROEXCEPTION);
  160.                 Set_divisionbyzeroflag();
  161.                 Sgl_setinfinity_exponentmantissa(result);
  162. *dstptr = result;
  163. return(NOEXCEPTION);
  164. }
  165. /*
  166.  * Generate exponent 
  167.  */
  168. dest_exponent = Sgl_exponent(opnd1) - Sgl_exponent(opnd2) + SGL_BIAS;
  169. /*
  170.  * Generate mantissa
  171.  */
  172. if (Sgl_isnotzero_exponent(opnd1)) {
  173. /* set hidden bit */
  174. Sgl_clear_signexponent_set_hidden(opnd1);
  175. }
  176. else {
  177. /* check for zero */
  178. if (Sgl_iszero_mantissa(opnd1)) {
  179. Sgl_setzero_exponentmantissa(result);
  180. *dstptr = result;
  181. return(NOEXCEPTION);
  182. }
  183.                 /* is denormalized; want to normalize */
  184.                 Sgl_clear_signexponent(opnd1);
  185.                 Sgl_leftshiftby1(opnd1);
  186. Sgl_normalize(opnd1,dest_exponent);
  187. }
  188. /* opnd2 needs to have hidden bit set with msb in hidden bit */
  189. if (Sgl_isnotzero_exponent(opnd2)) {
  190. Sgl_clear_signexponent_set_hidden(opnd2);
  191. }
  192. else {
  193.                 /* is denormalized; want to normalize */
  194.                 Sgl_clear_signexponent(opnd2);
  195.                 Sgl_leftshiftby1(opnd2);
  196. while(Sgl_iszero_hiddenhigh7mantissa(opnd2)) {
  197. Sgl_leftshiftby8(opnd2);
  198. dest_exponent += 8;
  199. }
  200. if(Sgl_iszero_hiddenhigh3mantissa(opnd2)) {
  201. Sgl_leftshiftby4(opnd2);
  202. dest_exponent += 4;
  203. }
  204. while(Sgl_iszero_hidden(opnd2)) {
  205. Sgl_leftshiftby1(opnd2);
  206. dest_exponent += 1;
  207. }
  208. }
  209. /* Divide the source mantissas */
  210. /*
  211.  * A non_restoring divide algorithm is used.
  212.  */
  213. Sgl_subtract(opnd1,opnd2,opnd1);
  214. Sgl_setzero(opnd3);
  215. for (count=1;count<=SGL_P && Sgl_all(opnd1);count++) {
  216. Sgl_leftshiftby1(opnd1);
  217. Sgl_leftshiftby1(opnd3);
  218. if (Sgl_iszero_sign(opnd1)) {
  219. Sgl_setone_lowmantissa(opnd3);
  220. Sgl_subtract(opnd1,opnd2,opnd1);
  221. }
  222. else Sgl_addition(opnd1,opnd2,opnd1);
  223. }
  224. if (count <= SGL_P) {
  225. Sgl_leftshiftby1(opnd3);
  226. Sgl_setone_lowmantissa(opnd3);
  227. Sgl_leftshift(opnd3,SGL_P-count);
  228. if (Sgl_iszero_hidden(opnd3)) {
  229. Sgl_leftshiftby1(opnd3);
  230. dest_exponent--;
  231. }
  232. }
  233. else {
  234. if (Sgl_iszero_hidden(opnd3)) {
  235. /* need to get one more bit of result */
  236. Sgl_leftshiftby1(opnd1);
  237. Sgl_leftshiftby1(opnd3);
  238. if (Sgl_iszero_sign(opnd1)) {
  239. Sgl_setone_lowmantissa(opnd3);
  240. Sgl_subtract(opnd1,opnd2,opnd1);
  241. }
  242. else Sgl_addition(opnd1,opnd2,opnd1);
  243. dest_exponent--;
  244. }
  245. if (Sgl_iszero_sign(opnd1)) guardbit = TRUE;
  246. stickybit = Sgl_all(opnd1);
  247. }
  248. inexact = guardbit | stickybit;
  249. /* 
  250.  * round result 
  251.  */
  252. if (inexact && (dest_exponent > 0 || Is_underflowtrap_enabled())) {
  253. Sgl_clear_signexponent(opnd3);
  254. switch (Rounding_mode()) {
  255. case ROUNDPLUS: 
  256. if (Sgl_iszero_sign(result)) 
  257. Sgl_increment_mantissa(opnd3);
  258. break;
  259. case ROUNDMINUS: 
  260. if (Sgl_isone_sign(result)) 
  261. Sgl_increment_mantissa(opnd3);
  262. break;
  263. case ROUNDNEAREST:
  264. if (guardbit) {
  265.     if (stickybit || Sgl_isone_lowmantissa(opnd3))
  266.            Sgl_increment_mantissa(opnd3);
  267. }
  268. }
  269. if (Sgl_isone_hidden(opnd3)) dest_exponent++;
  270. }
  271. Sgl_set_mantissa(result,opnd3);
  272.         /* 
  273.          * Test for overflow
  274.          */
  275. if (dest_exponent >= SGL_INFINITY_EXPONENT) {
  276.                 /* trap if OVERFLOWTRAP enabled */
  277.                 if (Is_overflowtrap_enabled()) {
  278.                         /*
  279.                          * Adjust bias of result
  280.                          */
  281.                         Sgl_setwrapped_exponent(result,dest_exponent,ovfl);
  282.                         *dstptr = result;
  283.                         if (inexact) 
  284.                             if (Is_inexacttrap_enabled())
  285.                                 return(OVERFLOWEXCEPTION | INEXACTEXCEPTION);
  286.                             else Set_inexactflag();
  287.                         return(OVERFLOWEXCEPTION);
  288.                 }
  289. Set_overflowflag();
  290.                 /* set result to infinity or largest number */
  291. Sgl_setoverflow(result);
  292. inexact = TRUE;
  293. }
  294.         /* 
  295.          * Test for underflow
  296.          */
  297. else if (dest_exponent <= 0) {
  298.                 /* trap if UNDERFLOWTRAP enabled */
  299.                 if (Is_underflowtrap_enabled()) {
  300.                         /*
  301.                          * Adjust bias of result
  302.                          */
  303.                         Sgl_setwrapped_exponent(result,dest_exponent,unfl);
  304.                         *dstptr = result;
  305.                         if (inexact) 
  306.                             if (Is_inexacttrap_enabled())
  307.                                 return(UNDERFLOWEXCEPTION | INEXACTEXCEPTION);
  308.                             else Set_inexactflag();
  309.                         return(UNDERFLOWEXCEPTION);
  310.                 }
  311. /* Determine if should set underflow flag */
  312. is_tiny = TRUE;
  313. if (dest_exponent == 0 && inexact) {
  314. switch (Rounding_mode()) {
  315. case ROUNDPLUS: 
  316. if (Sgl_iszero_sign(result)) {
  317. Sgl_increment(opnd3);
  318. if (Sgl_isone_hiddenoverflow(opnd3))
  319.                      is_tiny = FALSE;
  320. Sgl_decrement(opnd3);
  321. }
  322. break;
  323. case ROUNDMINUS: 
  324. if (Sgl_isone_sign(result)) {
  325. Sgl_increment(opnd3);
  326. if (Sgl_isone_hiddenoverflow(opnd3))
  327.                      is_tiny = FALSE;
  328. Sgl_decrement(opnd3);
  329. }
  330. break;
  331. case ROUNDNEAREST:
  332. if (guardbit && (stickybit || 
  333.     Sgl_isone_lowmantissa(opnd3))) {
  334.        Sgl_increment(opnd3);
  335. if (Sgl_isone_hiddenoverflow(opnd3))
  336.                      is_tiny = FALSE;
  337. Sgl_decrement(opnd3);
  338. }
  339. break;
  340. }
  341. }
  342.                 /*
  343.                  * denormalize result or set to signed zero
  344.                  */
  345. stickybit = inexact;
  346. Sgl_denormalize(opnd3,dest_exponent,guardbit,stickybit,inexact);
  347. /* return rounded number */ 
  348. if (inexact) {
  349. switch (Rounding_mode()) {
  350. case ROUNDPLUS:
  351. if (Sgl_iszero_sign(result)) {
  352. Sgl_increment(opnd3);
  353. }
  354. break;
  355. case ROUNDMINUS: 
  356. if (Sgl_isone_sign(result))  {
  357. Sgl_increment(opnd3);
  358. }
  359. break;
  360. case ROUNDNEAREST:
  361. if (guardbit && (stickybit || 
  362.     Sgl_isone_lowmantissa(opnd3))) {
  363.        Sgl_increment(opnd3);
  364. }
  365. break;
  366. }
  367.                  if (is_tiny) Set_underflowflag();
  368.                 }
  369. Sgl_set_exponentmantissa(result,opnd3);
  370. }
  371. else Sgl_set_exponent(result,dest_exponent);
  372. *dstptr = result;
  373. /* check for inexact */
  374. if (inexact) {
  375. if (Is_inexacttrap_enabled()) return(INEXACTEXCEPTION);
  376. else  Set_inexactflag();
  377. }
  378. return(NOEXCEPTION);
  379. }