1178580Simp/* 2178580Simp * Copyright (c) 1994, 1995 Carnegie-Mellon University. 3178580Simp * All rights reserved. 4178580Simp * 5178580Simp * Author: Chris G. Demetriou 6178580Simp * 7178580Simp * Permission to use, copy, modify and distribute this software and 8178580Simp * its documentation is hereby granted, provided that both the copyright 9178580Simp * notice and this permission notice appear in all copies of the 10178580Simp * software, derivative works or modified versions, and any portions 11178580Simp * thereof, and that both notices appear in supporting documentation. 12178580Simp * 13178580Simp * CARNEGIE MELLON ALLOWS FREE USE OF THIS SOFTWARE IN ITS "AS IS" 14178580Simp * CONDITION. CARNEGIE MELLON DISCLAIMS ANY LIABILITY OF ANY KIND 15178580Simp * FOR ANY DAMAGES WHATSOEVER RESULTING FROM THE USE OF THIS SOFTWARE. 16178580Simp * 17178580Simp * Carnegie Mellon requests users of this software to return to 18178580Simp * 19178580Simp * Software Distribution Coordinator or Software.Distribution@CS.CMU.EDU 20178580Simp * School of Computer Science 21178580Simp * Carnegie Mellon University 22178580Simp * Pittsburgh PA 15213-3890 23178580Simp * 24178580Simp * any improvements or extensions that they make and grant Carnegie the 25178580Simp * rights to redistribute these changes. 26178580Simp * 27178580Simp * $NetBSD: modf.c,v 1.1 1995/02/10 17:50:25 cgd Exp $ 28178580Simp */ 29178580Simp 30178580Simp#include <sys/cdefs.h> 31178580Simp__FBSDID("$FreeBSD$"); 32178580Simp 33178580Simp#include <sys/types.h> 34178580Simp#include <errno.h> 35178580Simp#include <math.h> 36178580Simp#include <machine/ieee.h> 37178580Simp 38178580Simp/* 39178580Simp * double modf(double val, double *iptr) 40178580Simp * returns: f and i such that |f| < 1.0, (f + i) = val, and 41178580Simp * sign(f) == sign(i) == sign(val). 42178580Simp * 43178580Simp * Beware signedness when doing subtraction, and also operand size! 44178580Simp */ 45178580Simpdouble 46178580Simpmodf(val, iptr) 47178580Simp double val, *iptr; 48178580Simp{ 49178580Simp union doub { 50178580Simp double v; 51178580Simp struct ieee_double s; 52178580Simp } u, v; 53178580Simp u_int64_t frac; 54178580Simp 55178580Simp /* 56178580Simp * If input is Inf or NaN, return it and leave i alone. 57178580Simp */ 58178580Simp u.v = val; 59178580Simp if (u.s.dbl_exp == DBL_EXP_INFNAN) 60178580Simp return (u.v); 61178580Simp 62178580Simp /* 63178580Simp * If input can't have a fractional part, return 64178580Simp * (appropriately signed) zero, and make i be the input. 65178580Simp */ 66178580Simp if ((int)u.s.dbl_exp - DBL_EXP_BIAS > DBL_FRACBITS - 1) { 67178580Simp *iptr = u.v; 68178580Simp v.v = 0.0; 69178580Simp v.s.dbl_sign = u.s.dbl_sign; 70178580Simp return (v.v); 71178580Simp } 72178580Simp 73178580Simp /* 74178580Simp * If |input| < 1.0, return it, and set i to the appropriately 75178580Simp * signed zero. 76178580Simp */ 77178580Simp if (u.s.dbl_exp < DBL_EXP_BIAS) { 78178580Simp v.v = 0.0; 79178580Simp v.s.dbl_sign = u.s.dbl_sign; 80178580Simp *iptr = v.v; 81178580Simp return (u.v); 82178580Simp } 83178580Simp 84178580Simp /* 85178580Simp * There can be a fractional part of the input. 86178580Simp * If you look at the math involved for a few seconds, it's 87178580Simp * plain to see that the integral part is the input, with the 88178580Simp * low (DBL_FRACBITS - (exponent - DBL_EXP_BIAS)) bits zeroed, 89218909Sbrucec * the fractional part is the part with the rest of the 90178580Simp * bits zeroed. Just zeroing the high bits to get the 91178580Simp * fractional part would yield a fraction in need of 92178580Simp * normalization. Therefore, we take the easy way out, and 93178580Simp * just use subtraction to get the fractional part. 94178580Simp */ 95178580Simp v.v = u.v; 96178580Simp /* Zero the low bits of the fraction, the sleazy way. */ 97178580Simp frac = ((u_int64_t)v.s.dbl_frach << 32) + v.s.dbl_fracl; 98178580Simp frac >>= DBL_FRACBITS - (u.s.dbl_exp - DBL_EXP_BIAS); 99178580Simp frac <<= DBL_FRACBITS - (u.s.dbl_exp - DBL_EXP_BIAS); 100178580Simp v.s.dbl_fracl = frac & 0xffffffff; 101178580Simp v.s.dbl_frach = frac >> 32; 102178580Simp *iptr = v.v; 103178580Simp 104178580Simp u.v -= v.v; 105178580Simp u.s.dbl_sign = v.s.dbl_sign; 106178580Simp return (u.v); 107178580Simp} 108