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 #include "libm.h"
  31 #include "longdouble.h"
  32 
  33 /*
  34  * exp10l(x)
  35  *      n = nint(x*(log10/log2)) ;
  36  *      exp10(x) = 10**x = exp(x*ln(10)) = exp(n*ln2+(x*ln10-n*ln2))
  37  *               = 2**n*exp(ln10*(x-n*log2/log10)))
  38  *      If x is an integer <= M then use repeat multiplication. For
  39  *      10**M is the largest representable integer, where
  40  *              M = 10          single precision (24 bits)
  41  *              M = 22          double precision (53 bits)
  42  *              M = 48          quadruple precision (113 bits)
  43  */
  44 
  45 #define TINY    1.0e-20L        /* single: 1e-5, double: 1e-10, quad: 1e-20 */
  46 #define LG10OVT 4933.L          /* single:  39, double:  309, quad:  4933 */
  47 #define LG10UFT -4966.L         /* single: -45, double: -323, quad: -4966 */
  48 #define M       48
  49                         /* logt2hi : last 32 bits is zero for quad prec */
  50 #define LOGT2HI 0.30102999566398119521373889472420986034688L
  51 #define LOGT2LO 2.831664213089468167896664371953e-31L
  52 
  53 static const long double
  54         zero    = 0.0L,
  55         tiny    = TINY * TINY,
  56         one     = 1.0L,
  57         lg10    = 3.321928094887362347870319429489390175865e+0000L,
  58         ln10    = 2.302585092994045684017991454684364207601e+0000L,
  59         logt2hi = LOGT2HI,
  60         logt2lo = LOGT2LO,
  61         lg10ovt = LG10OVT,
  62         lg10uft = LG10UFT;
  63 
  64 long double
  65 exp10l(long double x) {
  66         long double t, tenp;
  67         int k;
  68 
  69         if (!finitel(x)) {
  70                 if (isnanl(x) || x > zero)
  71                         return (x + x);
  72                 else
  73                         return (zero);
  74         }
  75         if (fabsl(x) < tiny)
  76                 return (one + x);
  77         if (x <= lg10ovt)
  78                 if (x >= lg10uft) {
  79                         k = (int) x;
  80                         tenp = 10.0L;
  81                                         /* x is a small +integer */
  82                         if (0 <= k && k <= M && (long double) k == x) {
  83                                 t = one;
  84                                 if (k & 1)
  85                                         t *= tenp;
  86                                 k >>= 1;
  87                                 while (k) {
  88                                         tenp *= tenp;
  89                                         if (k & 1)
  90                                                 t *= tenp;
  91                                         k >>= 1;
  92                                 }
  93                                 return (t);
  94                         }
  95                         t = anintl(x * lg10);
  96                         return (scalbnl(expl(ln10 * ((x - t * logt2hi) -
  97                                 t * logt2lo)), (int) t));
  98                 } else
  99                         return (scalbnl(one, -50000));  /* underflow */
 100         else
 101                         return (scalbnl(one, 50000));   /* overflow  */
 102 }