gonzui

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