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 /*
  27  * Copyright 2006 Sun Microsystems, Inc.  All rights reserved.
  28  * Use is subject to license terms.
  29  */
  30 
  31 #pragma weak frexp = __frexp
  32 
  33 /*
  34  * frexp(x, exp) returns the normalized significand of x and sets
  35  * *exp so that x = r*2^(*exp) where r is the return value.  If x
  36  * is finite and nonzero, 1/2 <= |r| < 1.
  37  *
  38  * If x is zero, infinite or NaN, frexp returns x and sets *exp = 0.
  39  * (The relevant standards do not specify *exp when x is infinite or
  40  * NaN, but this code sets it anyway.)
  41  *
  42  * If x is a signaling NaN, this code returns x without attempting
  43  * to raise the invalid operation exception.  If x is subnormal,
  44  * this code treats it as nonzero regardless of nonstandard mode.
  45  */
  46 
  47 #include "libm.h"
  48 
  49 double
  50 __frexp(double x, int *exp)
  51 {
  52         union {
  53                 unsigned i[2];
  54                 double d;
  55         } xx, yy;
  56 
  57         double t;
  58         unsigned hx;
  59         int e;
  60 
  61         xx.d = x;
  62         hx = xx.i[HIWORD] & ~0x80000000;
  63 
  64         if (hx >= 0x7ff00000) {              /* x is infinite or NaN */
  65                 *exp = 0;
  66                 return (x);
  67         }
  68 
  69         e = 0;
  70 
  71         if (hx < 0x00100000) {               /* x is subnormal or zero */
  72                 if ((hx | xx.i[LOWORD]) == 0) {
  73                         *exp = 0;
  74                         return (x);
  75                 }
  76 
  77                 /*
  78                  * normalize x by regarding it as an integer
  79                  *
  80                  * Here we use 32-bit integer arithmetic to avoid trapping
  81                  * or emulating 64-bit arithmetic.  If 64-bit arithmetic is
  82                  * available (e.g., in SPARC V9), do this instead:
  83                  *
  84                  *  long lx = ((long) hx << 32) | xx.i[LOWORD];
  85                  *  xx.d = (xx.i[HIWORD] < 0)? -lx : lx;
  86                  *
  87                  * If subnormal arithmetic doesn't trap, just multiply x by
  88                  * a power of two.
  89                  */
  90                 yy.i[HIWORD] = 0x43300000 | hx;
  91                 yy.i[LOWORD] = xx.i[LOWORD];
  92                 t = yy.d;
  93                 yy.i[HIWORD] = 0x43300000;
  94                 yy.i[LOWORD] = 0;
  95                 t -= yy.d;              /* t = |x| scaled */
  96                 xx.d = ((int)xx.i[HIWORD] < 0) ? -t : t;
  97                 hx = xx.i[HIWORD] & ~0x80000000;
  98                 e = -1074;
  99         }
 100 
 101         /* now xx.d is normal */
 102         xx.i[HIWORD] = (xx.i[HIWORD] & ~0x7ff00000) | 0x3fe00000;
 103         *exp = e + (hx >> 20) - 0x3fe;
 104         return (xx.d);
 105 }