chromium / native_client / nacl-gcc / f80d6b9ee7f94755c697ffb7194fb01dd0c537dd / . / mpfr-2.4.1 / subnormal.c

/* mpfr_subnormalize -- Subnormalize a floating point number | |

emulating sub-normal numbers. | |

Copyright 2005, 2006, 2007, 2008, 2009 Free Software Foundation, Inc. | |

Contributed by the Arenaire and Cacao projects, INRIA. | |

This file is part of the GNU MPFR Library. | |

The GNU MPFR Library is free software; you can redistribute it and/or modify | |

it under the terms of the GNU Lesser General Public License as published by | |

the Free Software Foundation; either version 2.1 of the License, or (at your | |

option) any later version. | |

The GNU MPFR Library is distributed in the hope that it will be useful, but | |

WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY | |

or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public | |

License for more details. | |

You should have received a copy of the GNU Lesser General Public License | |

along with the GNU MPFR Library; see the file COPYING.LIB. If not, write to | |

the Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, | |

MA 02110-1301, USA. */ | |

#include "mpfr-impl.h" | |

/* For GMP_RNDN, we can have a problem of double rounding. | |

In such a case, this table helps to conclude what to do (y positive): | |

Rounding Bit | Sticky Bit | inexact | Action | new inexact | |

0 | ? | ? | Trunc | sticky | |

1 | 0 | 1 | Trunc | | |

1 | 0 | 0 | Trunc if even | | |

1 | 0 | -1 | AddOneUlp | | |

1 | 1 | ? | AddOneUlp | | |

For other rounding mode, there isn't such a problem. | |

Just round it again and merge the inexact flags. | |

*/ | |

int | |

mpfr_subnormalize (mpfr_ptr y, int old_inexact, mp_rnd_t rnd) | |

{ | |

int inexact = 0; | |

/* The subnormal exponent range are from: | |

mpfr_emin to mpfr_emin + MPFR_PREC(y) - 1 */ | |

if (MPFR_LIKELY (MPFR_IS_SINGULAR (y) | |

|| (MPFR_GET_EXP (y) >= | |

__gmpfr_emin + (mp_exp_t) MPFR_PREC (y) - 1))) | |

inexact = old_inexact; | |

/* We have to emulate one bit rounding if EXP(y) = emin */ | |

else if (MPFR_GET_EXP (y) == __gmpfr_emin) | |

{ | |

/* If this is a power of 2, we don't need rounding. | |

It handles cases when rouding away and y=0.1*2^emin */ | |

if (mpfr_powerof2_raw (y)) | |

inexact = old_inexact; | |

/* We keep the same sign for y. | |

Assuming Y is the real value and y the approximation | |

and since y is not a power of 2: 0.5*2^emin < Y < 1*2^emin | |

We also know the direction of the error thanks to inexact flag */ | |

else if (rnd == GMP_RNDN) | |

{ | |

mp_limb_t *mant, rb ,sb; | |

mp_size_t s; | |

/* We need the rounding bit and the sticky bit. Read them | |

and use the previous table to conclude. */ | |

s = MPFR_LIMB_SIZE (y) - 1; | |

mant = MPFR_MANT (y) + s; | |

rb = *mant & (MPFR_LIMB_HIGHBIT >> 1); | |

if (rb == 0) | |

goto set_min; | |

sb = *mant & ((MPFR_LIMB_HIGHBIT >> 1) - 1); | |

while (sb == 0 && s-- != 0) | |

sb = *--mant; | |

if (sb != 0) | |

goto set_min_p1; | |

/* Rounding bit is 1 and sticky bit is 0. | |

We need to examine old inexact flag to conclude. */ | |

if ((old_inexact > 0 && MPFR_SIGN (y) > 0) || | |

(old_inexact < 0 && MPFR_SIGN (y) < 0)) | |

goto set_min; | |

/* If inexact != 0, return 0.1*2^(emin+1). | |

Otherwise, rounding bit = 1, sticky bit = 0 and inexact = 0 | |

So we have 0.1100000000000000000000000*2^emin exactly. | |

We return 0.1*2^(emin+1) according to the even-rounding | |

rule on subnormals. */ | |

goto set_min_p1; | |

} | |

else if (MPFR_IS_LIKE_RNDZ (rnd, MPFR_IS_NEG (y))) | |

{ | |

set_min: | |

mpfr_setmin (y, __gmpfr_emin); | |

inexact = -MPFR_SIGN (y); | |

} | |

else | |

{ | |

set_min_p1: | |

/* Note: mpfr_setmin will abort if __gmpfr_emax == __gmpfr_emin. */ | |

mpfr_setmin (y, __gmpfr_emin + 1); | |

inexact = MPFR_SIGN (y); | |

} | |

} | |

else /* Hard case: It is more or less the same problem than mpfr_cache */ | |

{ | |

mpfr_t dest; | |

mp_prec_t q; | |

int sign; | |

/* Compute the intermediary precision */ | |

q = (mpfr_uexp_t) MPFR_GET_EXP (y) - __gmpfr_emin + 1; | |

mpfr_init2 (dest, q); | |

/* Round y in dest */ | |

sign = MPFR_SIGN (y); | |

MPFR_SET_EXP (dest, MPFR_GET_EXP (y)); | |

MPFR_SET_SIGN (dest, sign); | |

MPFR_RNDRAW_EVEN (inexact, dest, | |

MPFR_MANT (y), MPFR_PREC (y), rnd, sign, | |

MPFR_SET_EXP (dest, MPFR_GET_EXP (dest)+1)); | |

if (MPFR_LIKELY (old_inexact != 0)) | |

{ | |

if (MPFR_UNLIKELY(rnd == GMP_RNDN && (inexact == MPFR_EVEN_INEX | |

|| inexact == -MPFR_EVEN_INEX))) | |

{ | |

/* if both roundings are in the same direction, we have to go | |

back in the other direction */ | |

if (SAME_SIGN (inexact, old_inexact)) | |

{ | |

if (SAME_SIGN (inexact, MPFR_INT_SIGN (y))) | |

mpfr_nexttozero (dest); | |

else | |

mpfr_nexttoinf (dest); | |

inexact = -inexact; | |

} | |

} | |

else if (MPFR_UNLIKELY (inexact == 0)) | |

inexact = old_inexact; | |

} | |

old_inexact = mpfr_set (y, dest, rnd); | |

MPFR_ASSERTN (old_inexact == 0); | |

MPFR_ASSERTN (MPFR_IS_PURE_FP (y)); | |

mpfr_clear (dest); | |

} | |

return inexact; | |

} |