1 /*
   2  * CDDL HEADER START
   3  *
   4  * The contents of this file are subject to the terms of the
   5  * Common Development and Distribution License (the "License").
   6  * You may not use this file except in compliance with the License.
   7  *
   8  * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
   9  * or http://www.opensolaris.org/os/licensing.
  10  * See the License for the specific language governing permissions
  11  * and limitations under the License.
  12  *
  13  * When distributing Covered Code, include this CDDL HEADER in each
  14  * file and include the License file at usr/src/OPENSOLARIS.LICENSE.
  15  * If applicable, add the following below this CDDL HEADER, with the
  16  * fields enclosed by brackets "[]" replaced with your own identifying
  17  * information: Portions Copyright [yyyy] [name of copyright owner]
  18  *
  19  * CDDL HEADER END
  20  */
  21 
  22 /*
  23  * Copyright 2011 Nexenta Systems, Inc.  All rights reserved.
  24  */
  25 /*
  26  * Copyright 2006 Sun Microsystems, Inc.  All rights reserved.
  27  * Use is subject to license terms.
  28  */
  29 
  30 #pragma weak frexp = __frexp
  31 
  32 /*
  33  * frexp(x, exp) returns the normalized significand of x and sets
  34  * *exp so that x = r*2^(*exp) where r is the return value.  If x
  35  * is finite and nonzero, 1/2 <= |r| < 1.
  36  *
  37  * If x is zero, infinite or NaN, frexp returns x and sets *exp = 0.
  38  * (The relevant standards do not specify *exp when x is infinite or
  39  * NaN, but this code sets it anyway.)
  40  *
  41  * If x is a signaling NaN, this code returns x without attempting
  42  * to raise the invalid operation exception.  If x is subnormal,
  43  * this code treats it as nonzero regardless of nonstandard mode.
  44  */
  45 
  46 #include "libm.h"
  47 
  48 double
  49 __frexp(double x, int *exp) {
  50         union {
  51                 unsigned i[2];
  52                 double d;
  53         } xx, yy;
  54         double t;
  55         unsigned hx;
  56         int e;
  57 
  58         xx.d = x;
  59         hx = xx.i[HIWORD] & ~0x80000000;
  60 
  61         if (hx >= 0x7ff00000) { /* x is infinite or NaN */
  62                 *exp = 0;
  63                 return (x);
  64         }
  65 
  66         e = 0;
  67         if (hx < 0x00100000) { /* x is subnormal or zero */
  68                 if ((hx | xx.i[LOWORD]) == 0) {
  69                         *exp = 0;
  70                         return (x);
  71                 }
  72 
  73                 /*
  74                  * normalize x by regarding it as an integer
  75                  *
  76                  * Here we use 32-bit integer arithmetic to avoid trapping
  77                  * or emulating 64-bit arithmetic.  If 64-bit arithmetic is
  78                  * available (e.g., in SPARC V9), do this instead:
  79                  *
  80                  *  long lx = ((long) hx << 32) | xx.i[LOWORD];
  81                  *  xx.d = (xx.i[HIWORD] < 0)? -lx : lx;
  82                  *
  83                  * If subnormal arithmetic doesn't trap, just multiply x by
  84                  * a power of two.
  85                  */
  86                 yy.i[HIWORD] = 0x43300000 | hx;
  87                 yy.i[LOWORD] = xx.i[LOWORD];
  88                 t = yy.d;
  89                 yy.i[HIWORD] = 0x43300000;
  90                 yy.i[LOWORD] = 0;
  91                 t -= yy.d; /* t = |x| scaled */
  92                 xx.d = ((int)xx.i[HIWORD] < 0)? -t : t;
  93                 hx = xx.i[HIWORD] & ~0x80000000;
  94                 e = -1074;
  95         }
  96 
  97         /* now xx.d is normal */
  98         xx.i[HIWORD] = (xx.i[HIWORD] & ~0x7ff00000) | 0x3fe00000;
  99         *exp = e + (hx >> 20) - 0x3fe;
 100         return (xx.d);
 101 }