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