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