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