105436638acc7c010349a69c3395f1a57c642dc62Ying Wang/* Split a double into fraction and mantissa. 205436638acc7c010349a69c3395f1a57c642dc62Ying Wang Copyright (C) 2007-2012 Free Software Foundation, Inc. 305436638acc7c010349a69c3395f1a57c642dc62Ying Wang 405436638acc7c010349a69c3395f1a57c642dc62Ying Wang This program is free software: you can redistribute it and/or modify 505436638acc7c010349a69c3395f1a57c642dc62Ying Wang it under the terms of the GNU General Public License as published by 605436638acc7c010349a69c3395f1a57c642dc62Ying Wang the Free Software Foundation; either version 3 of the License, or 705436638acc7c010349a69c3395f1a57c642dc62Ying Wang (at your option) any later version. 805436638acc7c010349a69c3395f1a57c642dc62Ying Wang 905436638acc7c010349a69c3395f1a57c642dc62Ying Wang This program is distributed in the hope that it will be useful, 1005436638acc7c010349a69c3395f1a57c642dc62Ying Wang but WITHOUT ANY WARRANTY; without even the implied warranty of 1105436638acc7c010349a69c3395f1a57c642dc62Ying Wang MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 1205436638acc7c010349a69c3395f1a57c642dc62Ying Wang GNU General Public License for more details. 1305436638acc7c010349a69c3395f1a57c642dc62Ying Wang 1405436638acc7c010349a69c3395f1a57c642dc62Ying Wang You should have received a copy of the GNU General Public License 1505436638acc7c010349a69c3395f1a57c642dc62Ying Wang along with this program. If not, see <http://www.gnu.org/licenses/>. */ 1605436638acc7c010349a69c3395f1a57c642dc62Ying Wang 1705436638acc7c010349a69c3395f1a57c642dc62Ying Wang/* Written by Paolo Bonzini <bonzini@gnu.org>, 2003, and 1805436638acc7c010349a69c3395f1a57c642dc62Ying Wang Bruno Haible <bruno@clisp.org>, 2007. */ 1905436638acc7c010349a69c3395f1a57c642dc62Ying Wang 2005436638acc7c010349a69c3395f1a57c642dc62Ying Wang#if ! defined USE_LONG_DOUBLE 2105436638acc7c010349a69c3395f1a57c642dc62Ying Wang# include <config.h> 2205436638acc7c010349a69c3395f1a57c642dc62Ying Wang#endif 2305436638acc7c010349a69c3395f1a57c642dc62Ying Wang 2405436638acc7c010349a69c3395f1a57c642dc62Ying Wang/* Specification. */ 2505436638acc7c010349a69c3395f1a57c642dc62Ying Wang#include <math.h> 2605436638acc7c010349a69c3395f1a57c642dc62Ying Wang 2705436638acc7c010349a69c3395f1a57c642dc62Ying Wang#include <float.h> 2805436638acc7c010349a69c3395f1a57c642dc62Ying Wang#ifdef USE_LONG_DOUBLE 2905436638acc7c010349a69c3395f1a57c642dc62Ying Wang# include "isnanl-nolibm.h" 3005436638acc7c010349a69c3395f1a57c642dc62Ying Wang# include "fpucw.h" 3105436638acc7c010349a69c3395f1a57c642dc62Ying Wang#else 3205436638acc7c010349a69c3395f1a57c642dc62Ying Wang# include "isnand-nolibm.h" 3305436638acc7c010349a69c3395f1a57c642dc62Ying Wang#endif 3405436638acc7c010349a69c3395f1a57c642dc62Ying Wang 3505436638acc7c010349a69c3395f1a57c642dc62Ying Wang/* This file assumes FLT_RADIX = 2. If FLT_RADIX is a power of 2 greater 3605436638acc7c010349a69c3395f1a57c642dc62Ying Wang than 2, or not even a power of 2, some rounding errors can occur, so that 3705436638acc7c010349a69c3395f1a57c642dc62Ying Wang then the returned mantissa is only guaranteed to be <= 1.0, not < 1.0. */ 3805436638acc7c010349a69c3395f1a57c642dc62Ying Wang 3905436638acc7c010349a69c3395f1a57c642dc62Ying Wang#ifdef USE_LONG_DOUBLE 4005436638acc7c010349a69c3395f1a57c642dc62Ying Wang# define FUNC frexpl 4105436638acc7c010349a69c3395f1a57c642dc62Ying Wang# define DOUBLE long double 4205436638acc7c010349a69c3395f1a57c642dc62Ying Wang# define ISNAN isnanl 4305436638acc7c010349a69c3395f1a57c642dc62Ying Wang# define DECL_ROUNDING DECL_LONG_DOUBLE_ROUNDING 4405436638acc7c010349a69c3395f1a57c642dc62Ying Wang# define BEGIN_ROUNDING() BEGIN_LONG_DOUBLE_ROUNDING () 4505436638acc7c010349a69c3395f1a57c642dc62Ying Wang# define END_ROUNDING() END_LONG_DOUBLE_ROUNDING () 4605436638acc7c010349a69c3395f1a57c642dc62Ying Wang# define L_(literal) literal##L 4705436638acc7c010349a69c3395f1a57c642dc62Ying Wang#else 4805436638acc7c010349a69c3395f1a57c642dc62Ying Wang# define FUNC frexp 4905436638acc7c010349a69c3395f1a57c642dc62Ying Wang# define DOUBLE double 5005436638acc7c010349a69c3395f1a57c642dc62Ying Wang# define ISNAN isnand 5105436638acc7c010349a69c3395f1a57c642dc62Ying Wang# define DECL_ROUNDING 5205436638acc7c010349a69c3395f1a57c642dc62Ying Wang# define BEGIN_ROUNDING() 5305436638acc7c010349a69c3395f1a57c642dc62Ying Wang# define END_ROUNDING() 5405436638acc7c010349a69c3395f1a57c642dc62Ying Wang# define L_(literal) literal 5505436638acc7c010349a69c3395f1a57c642dc62Ying Wang#endif 5605436638acc7c010349a69c3395f1a57c642dc62Ying Wang 5705436638acc7c010349a69c3395f1a57c642dc62Ying WangDOUBLE 5805436638acc7c010349a69c3395f1a57c642dc62Ying WangFUNC (DOUBLE x, int *expptr) 5905436638acc7c010349a69c3395f1a57c642dc62Ying Wang{ 6005436638acc7c010349a69c3395f1a57c642dc62Ying Wang int sign; 6105436638acc7c010349a69c3395f1a57c642dc62Ying Wang int exponent; 6205436638acc7c010349a69c3395f1a57c642dc62Ying Wang DECL_ROUNDING 6305436638acc7c010349a69c3395f1a57c642dc62Ying Wang 6405436638acc7c010349a69c3395f1a57c642dc62Ying Wang /* Test for NaN, infinity, and zero. */ 6505436638acc7c010349a69c3395f1a57c642dc62Ying Wang if (ISNAN (x) || x + x == x) 6605436638acc7c010349a69c3395f1a57c642dc62Ying Wang { 6705436638acc7c010349a69c3395f1a57c642dc62Ying Wang *expptr = 0; 6805436638acc7c010349a69c3395f1a57c642dc62Ying Wang return x; 6905436638acc7c010349a69c3395f1a57c642dc62Ying Wang } 7005436638acc7c010349a69c3395f1a57c642dc62Ying Wang 7105436638acc7c010349a69c3395f1a57c642dc62Ying Wang sign = 0; 7205436638acc7c010349a69c3395f1a57c642dc62Ying Wang if (x < 0) 7305436638acc7c010349a69c3395f1a57c642dc62Ying Wang { 7405436638acc7c010349a69c3395f1a57c642dc62Ying Wang x = - x; 7505436638acc7c010349a69c3395f1a57c642dc62Ying Wang sign = -1; 7605436638acc7c010349a69c3395f1a57c642dc62Ying Wang } 7705436638acc7c010349a69c3395f1a57c642dc62Ying Wang 7805436638acc7c010349a69c3395f1a57c642dc62Ying Wang BEGIN_ROUNDING (); 7905436638acc7c010349a69c3395f1a57c642dc62Ying Wang 8005436638acc7c010349a69c3395f1a57c642dc62Ying Wang { 8105436638acc7c010349a69c3395f1a57c642dc62Ying Wang /* Since the exponent is an 'int', it fits in 64 bits. Therefore the 8205436638acc7c010349a69c3395f1a57c642dc62Ying Wang loops are executed no more than 64 times. */ 8305436638acc7c010349a69c3395f1a57c642dc62Ying Wang DOUBLE pow2[64]; /* pow2[i] = 2^2^i */ 8405436638acc7c010349a69c3395f1a57c642dc62Ying Wang DOUBLE powh[64]; /* powh[i] = 2^-2^i */ 8505436638acc7c010349a69c3395f1a57c642dc62Ying Wang int i; 8605436638acc7c010349a69c3395f1a57c642dc62Ying Wang 8705436638acc7c010349a69c3395f1a57c642dc62Ying Wang exponent = 0; 8805436638acc7c010349a69c3395f1a57c642dc62Ying Wang if (x >= L_(1.0)) 8905436638acc7c010349a69c3395f1a57c642dc62Ying Wang { 9005436638acc7c010349a69c3395f1a57c642dc62Ying Wang /* A positive exponent. */ 9105436638acc7c010349a69c3395f1a57c642dc62Ying Wang DOUBLE pow2_i; /* = pow2[i] */ 9205436638acc7c010349a69c3395f1a57c642dc62Ying Wang DOUBLE powh_i; /* = powh[i] */ 9305436638acc7c010349a69c3395f1a57c642dc62Ying Wang 9405436638acc7c010349a69c3395f1a57c642dc62Ying Wang /* Invariants: pow2_i = 2^2^i, powh_i = 2^-2^i, 9505436638acc7c010349a69c3395f1a57c642dc62Ying Wang x * 2^exponent = argument, x >= 1.0. */ 9605436638acc7c010349a69c3395f1a57c642dc62Ying Wang for (i = 0, pow2_i = L_(2.0), powh_i = L_(0.5); 9705436638acc7c010349a69c3395f1a57c642dc62Ying Wang ; 9805436638acc7c010349a69c3395f1a57c642dc62Ying Wang i++, pow2_i = pow2_i * pow2_i, powh_i = powh_i * powh_i) 9905436638acc7c010349a69c3395f1a57c642dc62Ying Wang { 10005436638acc7c010349a69c3395f1a57c642dc62Ying Wang if (x >= pow2_i) 10105436638acc7c010349a69c3395f1a57c642dc62Ying Wang { 10205436638acc7c010349a69c3395f1a57c642dc62Ying Wang exponent += (1 << i); 10305436638acc7c010349a69c3395f1a57c642dc62Ying Wang x *= powh_i; 10405436638acc7c010349a69c3395f1a57c642dc62Ying Wang } 10505436638acc7c010349a69c3395f1a57c642dc62Ying Wang else 10605436638acc7c010349a69c3395f1a57c642dc62Ying Wang break; 10705436638acc7c010349a69c3395f1a57c642dc62Ying Wang 10805436638acc7c010349a69c3395f1a57c642dc62Ying Wang pow2[i] = pow2_i; 10905436638acc7c010349a69c3395f1a57c642dc62Ying Wang powh[i] = powh_i; 11005436638acc7c010349a69c3395f1a57c642dc62Ying Wang } 11105436638acc7c010349a69c3395f1a57c642dc62Ying Wang /* Avoid making x too small, as it could become a denormalized 11205436638acc7c010349a69c3395f1a57c642dc62Ying Wang number and thus lose precision. */ 11305436638acc7c010349a69c3395f1a57c642dc62Ying Wang while (i > 0 && x < pow2[i - 1]) 11405436638acc7c010349a69c3395f1a57c642dc62Ying Wang { 11505436638acc7c010349a69c3395f1a57c642dc62Ying Wang i--; 11605436638acc7c010349a69c3395f1a57c642dc62Ying Wang powh_i = powh[i]; 11705436638acc7c010349a69c3395f1a57c642dc62Ying Wang } 11805436638acc7c010349a69c3395f1a57c642dc62Ying Wang exponent += (1 << i); 11905436638acc7c010349a69c3395f1a57c642dc62Ying Wang x *= powh_i; 12005436638acc7c010349a69c3395f1a57c642dc62Ying Wang /* Here 2^-2^i <= x < 1.0. */ 12105436638acc7c010349a69c3395f1a57c642dc62Ying Wang } 12205436638acc7c010349a69c3395f1a57c642dc62Ying Wang else 12305436638acc7c010349a69c3395f1a57c642dc62Ying Wang { 12405436638acc7c010349a69c3395f1a57c642dc62Ying Wang /* A negative or zero exponent. */ 12505436638acc7c010349a69c3395f1a57c642dc62Ying Wang DOUBLE pow2_i; /* = pow2[i] */ 12605436638acc7c010349a69c3395f1a57c642dc62Ying Wang DOUBLE powh_i; /* = powh[i] */ 12705436638acc7c010349a69c3395f1a57c642dc62Ying Wang 12805436638acc7c010349a69c3395f1a57c642dc62Ying Wang /* Invariants: pow2_i = 2^2^i, powh_i = 2^-2^i, 12905436638acc7c010349a69c3395f1a57c642dc62Ying Wang x * 2^exponent = argument, x < 1.0. */ 13005436638acc7c010349a69c3395f1a57c642dc62Ying Wang for (i = 0, pow2_i = L_(2.0), powh_i = L_(0.5); 13105436638acc7c010349a69c3395f1a57c642dc62Ying Wang ; 13205436638acc7c010349a69c3395f1a57c642dc62Ying Wang i++, pow2_i = pow2_i * pow2_i, powh_i = powh_i * powh_i) 13305436638acc7c010349a69c3395f1a57c642dc62Ying Wang { 13405436638acc7c010349a69c3395f1a57c642dc62Ying Wang if (x < powh_i) 13505436638acc7c010349a69c3395f1a57c642dc62Ying Wang { 13605436638acc7c010349a69c3395f1a57c642dc62Ying Wang exponent -= (1 << i); 13705436638acc7c010349a69c3395f1a57c642dc62Ying Wang x *= pow2_i; 13805436638acc7c010349a69c3395f1a57c642dc62Ying Wang } 13905436638acc7c010349a69c3395f1a57c642dc62Ying Wang else 14005436638acc7c010349a69c3395f1a57c642dc62Ying Wang break; 14105436638acc7c010349a69c3395f1a57c642dc62Ying Wang 14205436638acc7c010349a69c3395f1a57c642dc62Ying Wang pow2[i] = pow2_i; 14305436638acc7c010349a69c3395f1a57c642dc62Ying Wang powh[i] = powh_i; 14405436638acc7c010349a69c3395f1a57c642dc62Ying Wang } 14505436638acc7c010349a69c3395f1a57c642dc62Ying Wang /* Here 2^-2^i <= x < 1.0. */ 14605436638acc7c010349a69c3395f1a57c642dc62Ying Wang } 14705436638acc7c010349a69c3395f1a57c642dc62Ying Wang 14805436638acc7c010349a69c3395f1a57c642dc62Ying Wang /* Invariants: x * 2^exponent = argument, and 2^-2^i <= x < 1.0. */ 14905436638acc7c010349a69c3395f1a57c642dc62Ying Wang while (i > 0) 15005436638acc7c010349a69c3395f1a57c642dc62Ying Wang { 15105436638acc7c010349a69c3395f1a57c642dc62Ying Wang i--; 15205436638acc7c010349a69c3395f1a57c642dc62Ying Wang if (x < powh[i]) 15305436638acc7c010349a69c3395f1a57c642dc62Ying Wang { 15405436638acc7c010349a69c3395f1a57c642dc62Ying Wang exponent -= (1 << i); 15505436638acc7c010349a69c3395f1a57c642dc62Ying Wang x *= pow2[i]; 15605436638acc7c010349a69c3395f1a57c642dc62Ying Wang } 15705436638acc7c010349a69c3395f1a57c642dc62Ying Wang } 15805436638acc7c010349a69c3395f1a57c642dc62Ying Wang /* Here 0.5 <= x < 1.0. */ 15905436638acc7c010349a69c3395f1a57c642dc62Ying Wang } 16005436638acc7c010349a69c3395f1a57c642dc62Ying Wang 16105436638acc7c010349a69c3395f1a57c642dc62Ying Wang if (sign < 0) 16205436638acc7c010349a69c3395f1a57c642dc62Ying Wang x = - x; 16305436638acc7c010349a69c3395f1a57c642dc62Ying Wang 16405436638acc7c010349a69c3395f1a57c642dc62Ying Wang END_ROUNDING (); 16505436638acc7c010349a69c3395f1a57c642dc62Ying Wang 16605436638acc7c010349a69c3395f1a57c642dc62Ying Wang *expptr = exponent; 16705436638acc7c010349a69c3395f1a57c642dc62Ying Wang return x; 16805436638acc7c010349a69c3395f1a57c642dc62Ying Wang} 169