1251881Speter/* Copyright (C) 1995,1996,1997,1998,1999,2002,2003 2251881Speter Free Software Foundation, Inc. 3251881Speter This file is part of the GNU C Library. 4251881Speter 5251881Speter The GNU C Library is free software; you can redistribute it and/or 6251881Speter modify it under the terms of the GNU Lesser General Public 7251881Speter License as published by the Free Software Foundation; either 8251881Speter version 2.1 of the License, or (at your option) any later version. 9251881Speter 10251881Speter The GNU C Library is distributed in the hope that it will be useful, 11251881Speter but WITHOUT ANY WARRANTY; without even the implied warranty of 12251881Speter MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 13251881Speter Lesser General Public License for more details. 14251881Speter 15251881Speter You should have received a copy of the GNU Lesser General Public 16251881Speter License along with the GNU C Library; if not, write to the Free 17251881Speter Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 18251881Speter 02111-1307 USA. */ 19251881Speter 20251881Speter#include <float.h> 21251881Speter#include <math.h> 22251881Speter#include <stdlib.h> 23251881Speter#include "gmp-impl.h" 24251881Speter 25251881Speter/* Convert a `__float128' in IEEE854 quad-precision format to a 26251881Speter multi-precision integer representing the significand scaled up by its 27251881Speter number of bits (113 for long double) and an integral power of two 28251881Speter (MPN frexpl). */ 29251881Speter 30251881Spetermp_size_t 31251881Spetermpn_extract_flt128 (mp_ptr res_ptr, mp_size_t size, 32251881Speter int *expt, int *is_neg, 33251881Speter __float128 value) 34251881Speter{ 35251881Speter ieee854_float128 u; 36251881Speter u.value = value; 37251881Speter 38251881Speter *is_neg = u.ieee.negative; 39251881Speter *expt = (int) u.ieee.exponent - IEEE854_FLOAT128_BIAS; 40251881Speter 41251881Speter#if BITS_PER_MP_LIMB == 32 42251881Speter res_ptr[0] = u.ieee.mant_low; /* Low-order 32 bits of fraction. */ 43251881Speter res_ptr[1] = (u.ieee.mant_low >> 32); 44251881Speter res_ptr[2] = u.ieee.mant_high; 45251881Speter res_ptr[3] = (u.ieee.mant_high >> 32); /* High-order 32 bits. */ 46251881Speter #define N 4 47251881Speter#elif BITS_PER_MP_LIMB == 64 48251881Speter res_ptr[0] = u.ieee.mant_low; 49251881Speter res_ptr[1] = u.ieee.mant_high; 50251881Speter #define N 2 51251881Speter#else 52251881Speter #error "mp_limb size " BITS_PER_MP_LIMB "not accounted for" 53251881Speter#endif 54251881Speter/* The format does not fill the last limb. There are some zeros. */ 55251881Speter#define NUM_LEADING_ZEROS (BITS_PER_MP_LIMB \ 56251881Speter - (FLT128_MANT_DIG - ((N - 1) * BITS_PER_MP_LIMB))) 57251881Speter 58251881Speter if (u.ieee.exponent == 0) 59251881Speter { 60251881Speter /* A biased exponent of zero is a special case. 61251881Speter Either it is a zero or it is a denormal number. */ 62251881Speter if (res_ptr[0] == 0 && res_ptr[1] == 0 63251881Speter && res_ptr[N - 2] == 0 && res_ptr[N - 1] == 0) /* Assumes N<=4. */ 64251881Speter /* It's zero. */ 65251881Speter *expt = 0; 66251881Speter else 67251881Speter { 68251881Speter /* It is a denormal number, meaning it has no implicit leading 69251881Speter one bit, and its exponent is in fact the format minimum. */ 70251881Speter int cnt; 71251881Speter 72251881Speter#if N == 2 73251881Speter if (res_ptr[N - 1] != 0) 74251881Speter { 75251881Speter count_leading_zeros (cnt, res_ptr[N - 1]); 76251881Speter cnt -= NUM_LEADING_ZEROS; 77251881Speter res_ptr[N - 1] = res_ptr[N - 1] << cnt 78251881Speter | (res_ptr[0] >> (BITS_PER_MP_LIMB - cnt)); 79251881Speter res_ptr[0] <<= cnt; 80251881Speter *expt = FLT128_MIN_EXP - 1 - cnt; 81251881Speter } 82251881Speter else 83251881Speter { 84251881Speter count_leading_zeros (cnt, res_ptr[0]); 85251881Speter if (cnt >= NUM_LEADING_ZEROS) 86251881Speter { 87251881Speter res_ptr[N - 1] = res_ptr[0] << (cnt - NUM_LEADING_ZEROS); 88251881Speter res_ptr[0] = 0; 89251881Speter } 90251881Speter else 91251881Speter { 92251881Speter res_ptr[N - 1] = res_ptr[0] >> (NUM_LEADING_ZEROS - cnt); 93251881Speter res_ptr[0] <<= BITS_PER_MP_LIMB - (NUM_LEADING_ZEROS - cnt); 94251881Speter } 95251881Speter *expt = FLT128_MIN_EXP - 1 96251881Speter - (BITS_PER_MP_LIMB - NUM_LEADING_ZEROS) - cnt; 97251881Speter } 98251881Speter#else 99251881Speter int j, k, l; 100251881Speter 101251881Speter for (j = N - 1; j > 0; j--) 102251881Speter if (res_ptr[j] != 0) 103251881Speter break; 104251881Speter 105251881Speter count_leading_zeros (cnt, res_ptr[j]); 106251881Speter cnt -= NUM_LEADING_ZEROS; 107251881Speter l = N - 1 - j; 108251881Speter if (cnt < 0) 109251881Speter { 110251881Speter cnt += BITS_PER_MP_LIMB; 111251881Speter l--; 112251881Speter } 113251881Speter if (!cnt) 114251881Speter for (k = N - 1; k >= l; k--) 115251881Speter res_ptr[k] = res_ptr[k-l]; 116251881Speter else 117251881Speter { 118251881Speter for (k = N - 1; k > l; k--) 119251881Speter res_ptr[k] = res_ptr[k-l] << cnt 120251881Speter | res_ptr[k-l-1] >> (BITS_PER_MP_LIMB - cnt); 121251881Speter res_ptr[k--] = res_ptr[0] << cnt; 122251881Speter } 123251881Speter 124251881Speter for (; k >= 0; k--) 125251881Speter res_ptr[k] = 0; 126251881Speter *expt = FLT128_MIN_EXP - 1 - l * BITS_PER_MP_LIMB - cnt; 127251881Speter#endif 128251881Speter } 129251881Speter } 130251881Speter else 131251881Speter /* Add the implicit leading one bit for a normalized number. */ 132251881Speter res_ptr[N - 1] |= (mp_limb_t) 1 << (FLT128_MANT_DIG - 1 133251881Speter - ((N - 1) * BITS_PER_MP_LIMB)); 134251881Speter 135251881Speter return N; 136251881Speter} 137251881Speter