1165855Sdas/*- 2165855Sdas * Copyright (c) 2007 David Schultz <das@FreeBSD.ORG> 3165855Sdas * All rights reserved. 4165855Sdas * 5165855Sdas * Redistribution and use in source and binary forms, with or without 6165855Sdas * modification, are permitted provided that the following conditions 7165855Sdas * are met: 8165855Sdas * 1. Redistributions of source code must retain the above copyright 9165855Sdas * notice, this list of conditions and the following disclaimer. 10165855Sdas * 2. Redistributions in binary form must reproduce the above copyright 11165855Sdas * notice, this list of conditions and the following disclaimer in the 12165855Sdas * documentation and/or other materials provided with the distribution. 13165855Sdas * 14165855Sdas * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND 15165855Sdas * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE 16165855Sdas * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE 17165855Sdas * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE 18165855Sdas * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL 19165855Sdas * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS 20165855Sdas * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) 21165855Sdas * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT 22165855Sdas * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY 23165855Sdas * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF 24165855Sdas * SUCH DAMAGE. 25165855Sdas * 26165855Sdas * Derived from s_modf.c, which has the following Copyright: 27165855Sdas * ==================================================== 28165855Sdas * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. 29165855Sdas * 30165855Sdas * Developed at SunPro, a Sun Microsystems, Inc. business. 31165855Sdas * Permission to use, copy, modify, and distribute this 32165855Sdas * software is freely granted, provided that this notice 33165855Sdas * is preserved. 34165855Sdas * ==================================================== 35165855Sdas * 36165855Sdas * $FreeBSD: releng/10.2/lib/msun/src/s_modfl.c 165855 2007-01-07 07:54:21Z das $ 37165855Sdas */ 38165855Sdas 39165855Sdas#include <float.h> 40165855Sdas#include <math.h> 41165855Sdas#include <sys/types.h> 42165855Sdas 43165855Sdas#include "fpmath.h" 44165855Sdas 45165855Sdas#if LDBL_MANL_SIZE > 32 46165855Sdas#define MASK ((uint64_t)-1) 47165855Sdas#else 48165855Sdas#define MASK ((uint32_t)-1) 49165855Sdas#endif 50165855Sdas/* Return the last n bits of a word, representing the fractional part. */ 51165855Sdas#define GETFRAC(bits, n) ((bits) & ~(MASK << (n))) 52165855Sdas/* The number of fraction bits in manh, not counting the integer bit */ 53165855Sdas#define HIBITS (LDBL_MANT_DIG - LDBL_MANL_SIZE) 54165855Sdas 55165855Sdasstatic const long double zero[] = { 0.0L, -0.0L }; 56165855Sdas 57165855Sdaslong double 58165855Sdasmodfl(long double x, long double *iptr) 59165855Sdas{ 60165855Sdas union IEEEl2bits ux; 61165855Sdas int e; 62165855Sdas 63165855Sdas ux.e = x; 64165855Sdas e = ux.bits.exp - LDBL_MAX_EXP + 1; 65165855Sdas if (e < HIBITS) { /* Integer part is in manh. */ 66165855Sdas if (e < 0) { /* |x|<1 */ 67165855Sdas *iptr = zero[ux.bits.sign]; 68165855Sdas return (x); 69165855Sdas } else { 70165855Sdas if ((GETFRAC(ux.bits.manh, HIBITS - 1 - e) | 71165855Sdas ux.bits.manl) == 0) { /* X is an integer. */ 72165855Sdas *iptr = x; 73165855Sdas return (zero[ux.bits.sign]); 74165855Sdas } else { 75165855Sdas /* Clear all but the top e+1 bits. */ 76165855Sdas ux.bits.manh >>= HIBITS - 1 - e; 77165855Sdas ux.bits.manh <<= HIBITS - 1 - e; 78165855Sdas ux.bits.manl = 0; 79165855Sdas *iptr = ux.e; 80165855Sdas return (x - ux.e); 81165855Sdas } 82165855Sdas } 83165855Sdas } else if (e >= LDBL_MANT_DIG - 1) { /* x has no fraction part. */ 84165855Sdas *iptr = x; 85165855Sdas if (x != x) /* Handle NaNs. */ 86165855Sdas return (x); 87165855Sdas return (zero[ux.bits.sign]); 88165855Sdas } else { /* Fraction part is in manl. */ 89165855Sdas if (GETFRAC(ux.bits.manl, LDBL_MANT_DIG - 1 - e) == 0) { 90165855Sdas /* x is integral. */ 91165855Sdas *iptr = x; 92165855Sdas return (zero[ux.bits.sign]); 93165855Sdas } else { 94165855Sdas /* Clear all but the top e+1 bits. */ 95165855Sdas ux.bits.manl >>= LDBL_MANT_DIG - 1 - e; 96165855Sdas ux.bits.manl <<= LDBL_MANT_DIG - 1 - e; 97165855Sdas *iptr = ux.e; 98165855Sdas return (x - ux.e); 99165855Sdas } 100165855Sdas } 101165855Sdas} 102