Home | History | Annotate | Line # | Download | only in gen
compat_modf_ieee754.c revision 1.5 
      1 /* $NetBSD: compat_modf_ieee754.c,v 1.5 2016/10/07 11:10:44 christos Exp $ */
      2 
      3 /*
      4  * Copyright (c) 1994, 1995 Carnegie-Mellon University.
      5  * All rights reserved.
      6  *
      7  * Author: Chris G. Demetriou
      8  *
      9  * Permission to use, copy, modify and distribute this software and
     10  * its documentation is hereby granted, provided that both the copyright
     11  * notice and this permission notice appear in all copies of the
     12  * software, derivative works or modified versions, and any portions
     13  * thereof, and that both notices appear in supporting documentation.
     14  *
     15  * CARNEGIE MELLON ALLOWS FREE USE OF THIS SOFTWARE IN ITS "AS IS"
     16  * CONDITION.  CARNEGIE MELLON DISCLAIMS ANY LIABILITY OF ANY KIND
     17  * FOR ANY DAMAGES WHATSOEVER RESULTING FROM THE USE OF THIS SOFTWARE.
     18  *
     19  * Carnegie Mellon requests users of this software to return to
     20  *
     21  *  Software Distribution Coordinator  or  Software.Distribution (at) CS.CMU.EDU
     22  *  School of Computer Science
     23  *  Carnegie Mellon University
     24  *  Pittsburgh PA 15213-3890
     25  *
     26  * any improvements or extensions that they make and grant Carnegie the
     27  * rights to redistribute these changes.
     28  */
     29 
     30 #include <sys/types.h>
     31 #include <machine/ieee.h>
     32 #include <errno.h>
     33 
     34 double modf(double, double *);
     35 
     36 /*
     37  * double modf(double val, double *iptr)
     38  * returns: f and i such that |f| < 1.0, (f + i) = val, and
     39  *	sign(f) == sign(i) == sign(val).
     40  *
     41  * Beware signedness when doing subtraction, and also operand size!
     42  */
     43 double
     44 modf(double val, double *iptr)
     45 {
     46 	union ieee_double_u u, v;
     47 	u_int64_t frac;
     48 
     49 	/*
     50 	 * If input is +/-Inf or NaN, return +/-0 or NaN.
     51 	 */
     52 	u.dblu_d = val;
     53 	if (u.dblu_dbl.dbl_exp == DBL_EXP_INFNAN) {
     54 		*iptr = u.dblu_d;
     55 		return (0.0 / u.dblu_d);
     56 	}
     57 
     58 	/*
     59 	 * If input can't have a fractional part, return
     60 	 * (appropriately signed) zero, and make i be the input.
     61 	 */
     62 	if ((int)u.dblu_dbl.dbl_exp - DBL_EXP_BIAS > DBL_FRACBITS - 1) {
     63 		*iptr = u.dblu_d;
     64 		v.dblu_d = 0.0;
     65 		v.dblu_dbl.dbl_sign = u.dblu_dbl.dbl_sign;
     66 		return (v.dblu_d);
     67 	}
     68 
     69 	/*
     70 	 * If |input| < 1.0, return it, and set i to the appropriately
     71 	 * signed zero.
     72 	 */
     73 	if (u.dblu_dbl.dbl_exp < DBL_EXP_BIAS) {
     74 		v.dblu_d = 0.0;
     75 		v.dblu_dbl.dbl_sign = u.dblu_dbl.dbl_sign;
     76 		*iptr = v.dblu_d;
     77 		return (u.dblu_d);
     78 	}
     79 
     80 	/*
     81 	 * There can be a fractional part of the input.
     82 	 * If you look at the math involved for a few seconds, it's
     83 	 * plain to see that the integral part is the input, with the
     84 	 * low (DBL_FRACBITS - (exponent - DBL_EXP_BIAS)) bits zeroed,
     85 	 * the fractional part is the part with the rest of the
     86 	 * bits zeroed.  Just zeroing the high bits to get the
     87 	 * fractional part would yield a fraction in need of
     88 	 * normalization.  Therefore, we take the easy way out, and
     89 	 * just use subtraction to get the fractional part.
     90 	 */
     91 	v.dblu_d = u.dblu_d;
     92 	/* Zero the low bits of the fraction, the sleazy way. */
     93 	frac = ((u_int64_t)v.dblu_dbl.dbl_frach << 32) + v.dblu_dbl.dbl_fracl;
     94 	frac >>= DBL_FRACBITS - (u.dblu_dbl.dbl_exp - DBL_EXP_BIAS);
     95 	frac <<= DBL_FRACBITS - (u.dblu_dbl.dbl_exp - DBL_EXP_BIAS);
     96 	v.dblu_dbl.dbl_fracl = (unsigned int)(frac & 0xffffffffULL);
     97 	v.dblu_dbl.dbl_frach = (unsigned int)(frac >> 32);
     98 	*iptr = v.dblu_d;
     99 
    100 	u.dblu_d -= v.dblu_d;
    101 	u.dblu_dbl.dbl_sign = v.dblu_dbl.dbl_sign;
    102 	return (u.dblu_d);
    103 }
    104