sfdiv.c 11 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392
  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., 59 Temple Place, Suite 330, Boston, MA 02111-1307 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. }