开通博客赚积分
发布资源赚积分
分类

源码开发语言/平台
当前位置：首页 > 源码/资料 > Linux/Unix编程 > 查看源码
dfdiv.c
资源名称：linux-2.4.20.tar.gz　[点击查看]
上传用户：jlfgdled
上传日期：2013-04-10
资源大小：33168k
文件大小：12k
源码类别：

Linux/Unix编程

开发平台：

Unix_Linux

dfdiv.c：源码内容
  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/dfdiv.c $Revision: 1.1 $
  26.  *
  27.  *  Purpose:
  28.  * Double Precision Floating-point Divide
  29.  *
  30.  *  External Interfaces:
  31.  * dbl_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. */
  42. #include "float.h"
  43. #include "dbl_float.h"
  45. /*
  46.  *  Double Precision Floating-point Divide
  47.  */
  49. int
  50. dbl_fdiv (dbl_floating_point * srcptr1, dbl_floating_point * srcptr2,
  51.   dbl_floating_point * dstptr, unsigned int *status)
  52. {
  53. register unsigned int opnd1p1, opnd1p2, opnd2p1, opnd2p2;
  54. register unsigned int opnd3p1, opnd3p2, resultp1, resultp2;
  55. register int dest_exponent, count;
  56. register boolean inexact = FALSE, guardbit = FALSE, stickybit = FALSE;
  57. boolean is_tiny;
  59. Dbl_copyfromptr(srcptr1,opnd1p1,opnd1p2);
  60. Dbl_copyfromptr(srcptr2,opnd2p1,opnd2p2);
  61. /* 
  62.  * set sign bit of result 
  63.  */
  64. if (Dbl_sign(opnd1p1) ^ Dbl_sign(opnd2p1)) 
  65. Dbl_setnegativezerop1(resultp1);  
  66. else Dbl_setzerop1(resultp1);
  67. /*
  68.  * check first operand for NaN's or infinity
  69.  */
  70. if (Dbl_isinfinity_exponent(opnd1p1)) {
  71. if (Dbl_iszero_mantissa(opnd1p1,opnd1p2)) {
  72. if (Dbl_isnotnan(opnd2p1,opnd2p2)) {
  73. if (Dbl_isinfinity(opnd2p1,opnd2p2)) {
  74. /* 
  75.  * invalid since both operands 
  76.  * are infinity 
  77.  */
  78. if (Is_invalidtrap_enabled())
  79.                                  return(INVALIDEXCEPTION);
  80.                                  Set_invalidflag();
  81.                                  Dbl_makequietnan(resultp1,resultp2);
  82. Dbl_copytoptr(resultp1,resultp2,dstptr);
  83. return(NOEXCEPTION);
  84. }
  85. /*
  86.    * return infinity
  87.    */
  88. Dbl_setinfinity_exponentmantissa(resultp1,resultp2);
  89. Dbl_copytoptr(resultp1,resultp2,dstptr);
  90. return(NOEXCEPTION);
  91. }
  92. }
  93. else {
  94.                  /*
  95.                    * is NaN; signaling or quiet?
  96.                    */
  97.                  if (Dbl_isone_signaling(opnd1p1)) {
  98.                          /* trap if INVALIDTRAP enabled */
  99.                          if (Is_invalidtrap_enabled())
  100.                              return(INVALIDEXCEPTION);
  101.                          /* make NaN quiet */
  102.                          Set_invalidflag();
  103.                          Dbl_set_quiet(opnd1p1);
  104.                  }
  105. /* 
  106.  * is second operand a signaling NaN? 
  107.  */
  108. else if (Dbl_is_signalingnan(opnd2p1)) {
  109.                          /* trap if INVALIDTRAP enabled */
  110.                          if (Is_invalidtrap_enabled())
  111.                              return(INVALIDEXCEPTION);
  112.                          /* make NaN quiet */
  113.                          Set_invalidflag();
  114.                          Dbl_set_quiet(opnd2p1);
  115. Dbl_copytoptr(opnd2p1,opnd2p2,dstptr);
  116.                  return(NOEXCEPTION);
  117. }
  118.                  /*
  119.                    * return quiet NaN
  120.                    */
  121. Dbl_copytoptr(opnd1p1,opnd1p2,dstptr);
  122.                  return(NOEXCEPTION);
  123. }
  124. }
  125. /*
  126.  * check second operand for NaN's or infinity
  127.  */
  128. if (Dbl_isinfinity_exponent(opnd2p1)) {
  129. if (Dbl_iszero_mantissa(opnd2p1,opnd2p2)) {
  130. /*
  131.  * return zero
  132.  */
  133. Dbl_setzero_exponentmantissa(resultp1,resultp2);
  134. Dbl_copytoptr(resultp1,resultp2,dstptr);
  135. return(NOEXCEPTION);
  136. }
  137.                 /*
  138.                  * is NaN; signaling or quiet?
  139.                  */
  140.                 if (Dbl_isone_signaling(opnd2p1)) {
  141.                         /* trap if INVALIDTRAP enabled */
  142.                         if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
  143.                         /* make NaN quiet */
  144.                         Set_invalidflag();
  145.                         Dbl_set_quiet(opnd2p1);
  146.                 }
  147.                 /*
  148.                  * return quiet NaN
  149.                  */
  150. Dbl_copytoptr(opnd2p1,opnd2p2,dstptr);
  151.                 return(NOEXCEPTION);
  152. }
  153.         /*
  154.          * check for division by zero
  155.          */
  156.         if (Dbl_iszero_exponentmantissa(opnd2p1,opnd2p2)) {
  157.                 if (Dbl_iszero_exponentmantissa(opnd1p1,opnd1p2)) {
  158.                         /* invalid since both operands are zero */
  159.                         if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
  160.                         Set_invalidflag();
  161.                         Dbl_makequietnan(resultp1,resultp2);
  162.                         Dbl_copytoptr(resultp1,resultp2,dstptr);
  163.                         return(NOEXCEPTION);
  164.                 }
  165.                 if (Is_divisionbyzerotrap_enabled())
  166.                         return(DIVISIONBYZEROEXCEPTION);
  167.                 Set_divisionbyzeroflag();
  168.                 Dbl_setinfinity_exponentmantissa(resultp1,resultp2);
  169.                 Dbl_copytoptr(resultp1,resultp2,dstptr);
  170.                 return(NOEXCEPTION);
  171.         }
  172. /*
  173.  * Generate exponent 
  174.  */
  175. dest_exponent = Dbl_exponent(opnd1p1) - Dbl_exponent(opnd2p1) + DBL_BIAS;
  177. /*
  178.  * Generate mantissa
  179.  */
  180. if (Dbl_isnotzero_exponent(opnd1p1)) {
  181. /* set hidden bit */
  182. Dbl_clear_signexponent_set_hidden(opnd1p1);
  183. }
  184. else {
  185. /* check for zero */
  186. if (Dbl_iszero_mantissa(opnd1p1,opnd1p2)) {
  187. Dbl_setzero_exponentmantissa(resultp1,resultp2);
  188. Dbl_copytoptr(resultp1,resultp2,dstptr);
  189. return(NOEXCEPTION);
  190. }
  191.                 /* is denormalized, want to normalize */
  192.                 Dbl_clear_signexponent(opnd1p1);
  193.                 Dbl_leftshiftby1(opnd1p1,opnd1p2);
  194. Dbl_normalize(opnd1p1,opnd1p2,dest_exponent);
  195. }
  196. /* opnd2 needs to have hidden bit set with msb in hidden bit */
  197. if (Dbl_isnotzero_exponent(opnd2p1)) {
  198. Dbl_clear_signexponent_set_hidden(opnd2p1);
  199. }
  200. else {
  201.                 /* is denormalized; want to normalize */
  202.                 Dbl_clear_signexponent(opnd2p1);
  203.                 Dbl_leftshiftby1(opnd2p1,opnd2p2);
  204.                 while (Dbl_iszero_hiddenhigh7mantissa(opnd2p1)) {
  205.                         dest_exponent+=8;
  206.                         Dbl_leftshiftby8(opnd2p1,opnd2p2);
  207.                 }
  208.                 if (Dbl_iszero_hiddenhigh3mantissa(opnd2p1)) {
  209.                         dest_exponent+=4;
  210.                         Dbl_leftshiftby4(opnd2p1,opnd2p2);
  211.                 }
  212.                 while (Dbl_iszero_hidden(opnd2p1)) {
  213.                         dest_exponent++;
  214.                         Dbl_leftshiftby1(opnd2p1,opnd2p2);
  215.                 }
  216. }
  218. /* Divide the source mantissas */
  220. /* 
  221.  * A non-restoring divide algorithm is used.
  222.  */
  223. Twoword_subtract(opnd1p1,opnd1p2,opnd2p1,opnd2p2);
  224. Dbl_setzero(opnd3p1,opnd3p2);
  225. for (count=1; count <= DBL_P && (opnd1p1 || opnd1p2); count++) {
  226. Dbl_leftshiftby1(opnd1p1,opnd1p2);
  227. Dbl_leftshiftby1(opnd3p1,opnd3p2);
  228. if (Dbl_iszero_sign(opnd1p1)) {
  229. Dbl_setone_lowmantissap2(opnd3p2);
  230. Twoword_subtract(opnd1p1,opnd1p2,opnd2p1,opnd2p2);
  231. }
  232. else {
  233. Twoword_add(opnd1p1, opnd1p2, opnd2p1, opnd2p2);
  234. }
  235. }
  236. if (count <= DBL_P) {
  237. Dbl_leftshiftby1(opnd3p1,opnd3p2);
  238. Dbl_setone_lowmantissap2(opnd3p2);
  239. Dbl_leftshift(opnd3p1,opnd3p2,(DBL_P-count));
  240. if (Dbl_iszero_hidden(opnd3p1)) {
  241. Dbl_leftshiftby1(opnd3p1,opnd3p2);
  242. dest_exponent--;
  243. }
  244. }
  245. else {
  246. if (Dbl_iszero_hidden(opnd3p1)) {
  247. /* need to get one more bit of result */
  248. Dbl_leftshiftby1(opnd1p1,opnd1p2);
  249. Dbl_leftshiftby1(opnd3p1,opnd3p2);
  250. if (Dbl_iszero_sign(opnd1p1)) {
  251. Dbl_setone_lowmantissap2(opnd3p2);
  252. Twoword_subtract(opnd1p1,opnd1p2,opnd2p1,opnd2p2);
  253. }
  254. else {
  255. Twoword_add(opnd1p1,opnd1p2,opnd2p1,opnd2p2);
  256. }
  257. dest_exponent--;
  258. }
  259. if (Dbl_iszero_sign(opnd1p1)) guardbit = TRUE;
  260. stickybit = Dbl_allp1(opnd1p1) || Dbl_allp2(opnd1p2);
  261. }
  262. inexact = guardbit | stickybit;
  264. /* 
  265.  * round result 
  266.  */
  267. if (inexact && (dest_exponent > 0 || Is_underflowtrap_enabled())) {
  268. Dbl_clear_signexponent(opnd3p1);
  269. switch (Rounding_mode()) {
  270. case ROUNDPLUS: 
  271. if (Dbl_iszero_sign(resultp1)) 
  272. Dbl_increment(opnd3p1,opnd3p2);
  273. break;
  274. case ROUNDMINUS: 
  275. if (Dbl_isone_sign(resultp1)) 
  276. Dbl_increment(opnd3p1,opnd3p2);
  277. break;
  278. case ROUNDNEAREST:
  279. if (guardbit && (stickybit || 
  280.     Dbl_isone_lowmantissap2(opnd3p2))) {
  281.        Dbl_increment(opnd3p1,opnd3p2);
  282. }
  283. }
  284. if (Dbl_isone_hidden(opnd3p1)) dest_exponent++;
  285. }
  286. Dbl_set_mantissa(resultp1,resultp2,opnd3p1,opnd3p2);
  288.         /* 
  289.          * Test for overflow
  290.          */
  291. if (dest_exponent >= DBL_INFINITY_EXPONENT) {
  292.                 /* trap if OVERFLOWTRAP enabled */
  293.                 if (Is_overflowtrap_enabled()) {
  294.                         /*
  295.                          * Adjust bias of result
  296.                          */
  297.                         Dbl_setwrapped_exponent(resultp1,dest_exponent,ovfl);
  298.                         Dbl_copytoptr(resultp1,resultp2,dstptr);
  299.                         if (inexact) 
  300.                             if (Is_inexacttrap_enabled())
  301.                                 return(OVERFLOWEXCEPTION | INEXACTEXCEPTION);
  302.                             else Set_inexactflag();
  303.                         return(OVERFLOWEXCEPTION);
  304.                 }
  305. Set_overflowflag();
  306.                 /* set result to infinity or largest number */
  307. Dbl_setoverflow(resultp1,resultp2);
  308. inexact = TRUE;
  309. }
  310.         /* 
  311.          * Test for underflow
  312.          */
  313. else if (dest_exponent <= 0) {
  314.                 /* trap if UNDERFLOWTRAP enabled */
  315.                 if (Is_underflowtrap_enabled()) {
  316.                         /*
  317.                          * Adjust bias of result
  318.                          */
  319.                         Dbl_setwrapped_exponent(resultp1,dest_exponent,unfl);
  320.                         Dbl_copytoptr(resultp1,resultp2,dstptr);
  321.                         if (inexact) 
  322.                             if (Is_inexacttrap_enabled())
  323.                                 return(UNDERFLOWEXCEPTION | INEXACTEXCEPTION);
  324.                             else Set_inexactflag();
  325.                         return(UNDERFLOWEXCEPTION);
  326.                 }
  328. /* Determine if should set underflow flag */
  329. is_tiny = TRUE;
  330. if (dest_exponent == 0 && inexact) {
  331. switch (Rounding_mode()) {
  332. case ROUNDPLUS: 
  333. if (Dbl_iszero_sign(resultp1)) {
  334. Dbl_increment(opnd3p1,opnd3p2);
  335. if (Dbl_isone_hiddenoverflow(opnd3p1))
  336.                      is_tiny = FALSE;
  337. Dbl_decrement(opnd3p1,opnd3p2);
  338. }
  339. break;
  340. case ROUNDMINUS: 
  341. if (Dbl_isone_sign(resultp1)) {
  342. Dbl_increment(opnd3p1,opnd3p2);
  343. if (Dbl_isone_hiddenoverflow(opnd3p1))
  344.                      is_tiny = FALSE;
  345. Dbl_decrement(opnd3p1,opnd3p2);
  346. }
  347. break;
  348. case ROUNDNEAREST:
  349. if (guardbit && (stickybit || 
  350.     Dbl_isone_lowmantissap2(opnd3p2))) {
  351.        Dbl_increment(opnd3p1,opnd3p2);
  352. if (Dbl_isone_hiddenoverflow(opnd3p1))
  353.                      is_tiny = FALSE;
  354. Dbl_decrement(opnd3p1,opnd3p2);
  355. }
  356. break;
  357. }
  358. }
  360.                 /*
  361.                  * denormalize result or set to signed zero
  362.                  */
  363. stickybit = inexact;
  364. Dbl_denormalize(opnd3p1,opnd3p2,dest_exponent,guardbit,
  365.  stickybit,inexact);
  367. /* return rounded number */ 
  368. if (inexact) {
  369. switch (Rounding_mode()) {
  370. case ROUNDPLUS:
  371. if (Dbl_iszero_sign(resultp1)) {
  372. Dbl_increment(opnd3p1,opnd3p2);
  373. }
  374. break;
  375. case ROUNDMINUS: 
  376. if (Dbl_isone_sign(resultp1)) {
  377. Dbl_increment(opnd3p1,opnd3p2);
  378. }
  379. break;
  380. case ROUNDNEAREST:
  381. if (guardbit && (stickybit || 
  382.     Dbl_isone_lowmantissap2(opnd3p2))) {
  383.        Dbl_increment(opnd3p1,opnd3p2);
  384. }
  385. break;
  386. }
  387.                  if (is_tiny) Set_underflowflag();
  388.                 }
  389. Dbl_set_exponentmantissa(resultp1,resultp2,opnd3p1,opnd3p2);
  390. }
  391. else Dbl_set_exponent(resultp1,dest_exponent);
  392. Dbl_copytoptr(resultp1,resultp2,dstptr);
  394. /* check for inexact */
  395. if (inexact) {
  396. if (Is_inexacttrap_enabled()) return(INEXACTEXCEPTION);
  397. else Set_inexactflag();
  398. }
  399. return(NOEXCEPTION);
  400. }