usr/src/lib/libm/common/C/exp2.c

   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  * Copyright 2011 Nexenta Systems, Inc.  All rights reserved.
  23  */
  24 /*
  25  * Copyright 2006 Sun Microsystems, Inc.  All rights reserved.
  26  * Use is subject to license terms.
  27  */
  28
  29 #pragma weak __exp2 = exp2
  30
  31 /* INDENT OFF */
  32 /*
  33  * exp2(x)
  34  * Code by K.C. Ng for SUN 4.0 libm.
  35  * Method :
  36  *      exp2(x) = 2**x = 2**((x-anint(x))+anint(x))
  37  *              = 2**anint(x)*2**(x-anint(x))
  38  *              = 2**anint(x)*exp((x-anint(x))*ln2)
  39  */
  40 /* INDENT ON */
  41
  42 #include "libm.h"
  43
  44 static const double C[] = {
  45         0.0,
  46         1.0,
  47         0.5,
  48         6.93147180559945286227e-01,
  49         1.0e300,
  50         1.0e-300,
  51 };
  52
  53 #define zero    C[0]
  54 #define one     C[1]
  55 #define half    C[2]
  56 #define ln2     C[3]
  57 #define huge    C[4]
  58 #define tiny    C[5]
  59
  60 double
  61 exp2(double x) {
  62         int     ix, hx, k;
  63         double  t;
  64
  65         ix = ((int *)&x)[HIWORD];
  66         hx = ix & ~0x80000000;
  67
  68         if (hx >= 0x4090e000) { /* |x| >= 1080 or x is nan */
  69                 if (hx >= 0x7ff00000) { /* x is inf or nan */
  70                         if (ix == 0xfff00000 && ((int *)&x)[LOWORD] == 0)
  71                                 return (zero);
  72                         return (x * x);
  73                 }
  74                 t = (ix < 0)? tiny : huge;
  75                 return (t * t);
  76         }
  77
  78         if (hx < 0x3fe00000) {  /* |x| < 0.5 */
  79                 if (hx < 0x3c000000)
  80                         return (one + x);
  81                 return (exp(ln2 * x));
  82         }
  83
  84         k = (int)x;
  85         if (x != (double)k)
  86                 k = (int)((ix < 0)? x - half : x + half);
  87         return (scalbn(exp(ln2 * (x - (double)k)), k));
  88 }