1*25c28e83SPiotr Jasiukajtis /* 2*25c28e83SPiotr Jasiukajtis * CDDL HEADER START 3*25c28e83SPiotr Jasiukajtis * 4*25c28e83SPiotr Jasiukajtis * The contents of this file are subject to the terms of the 5*25c28e83SPiotr Jasiukajtis * Common Development and Distribution License (the "License"). 6*25c28e83SPiotr Jasiukajtis * You may not use this file except in compliance with the License. 7*25c28e83SPiotr Jasiukajtis * 8*25c28e83SPiotr Jasiukajtis * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE 9*25c28e83SPiotr Jasiukajtis * or http://www.opensolaris.org/os/licensing. 10*25c28e83SPiotr Jasiukajtis * See the License for the specific language governing permissions 11*25c28e83SPiotr Jasiukajtis * and limitations under the License. 12*25c28e83SPiotr Jasiukajtis * 13*25c28e83SPiotr Jasiukajtis * When distributing Covered Code, include this CDDL HEADER in each 14*25c28e83SPiotr Jasiukajtis * file and include the License file at usr/src/OPENSOLARIS.LICENSE. 15*25c28e83SPiotr Jasiukajtis * If applicable, add the following below this CDDL HEADER, with the 16*25c28e83SPiotr Jasiukajtis * fields enclosed by brackets "[]" replaced with your own identifying 17*25c28e83SPiotr Jasiukajtis * information: Portions Copyright [yyyy] [name of copyright owner] 18*25c28e83SPiotr Jasiukajtis * 19*25c28e83SPiotr Jasiukajtis * CDDL HEADER END 20*25c28e83SPiotr Jasiukajtis */ 21*25c28e83SPiotr Jasiukajtis /* 22*25c28e83SPiotr Jasiukajtis * Copyright 2011 Nexenta Systems, Inc. All rights reserved. 23*25c28e83SPiotr Jasiukajtis */ 24*25c28e83SPiotr Jasiukajtis /* 25*25c28e83SPiotr Jasiukajtis * Copyright 2006 Sun Microsystems, Inc. All rights reserved. 26*25c28e83SPiotr Jasiukajtis * Use is subject to license terms. 27*25c28e83SPiotr Jasiukajtis */ 28*25c28e83SPiotr Jasiukajtis 29*25c28e83SPiotr Jasiukajtis /* INDENT OFF */ 30*25c28e83SPiotr Jasiukajtis /* 31*25c28e83SPiotr Jasiukajtis * exp10(x) 32*25c28e83SPiotr Jasiukajtis * Code by K.C. Ng for SUN 4.0 libm. 33*25c28e83SPiotr Jasiukajtis * Method : 34*25c28e83SPiotr Jasiukajtis * n = nint(x*(log10/log2)); 35*25c28e83SPiotr Jasiukajtis * exp10(x) = 10**x = exp(x*ln(10)) = exp(n*ln2+(x*ln10-n*ln2)) 36*25c28e83SPiotr Jasiukajtis * = 2**n*exp(ln10*(x-n*log2/log10))) 37*25c28e83SPiotr Jasiukajtis * If x is an integer < 23 then use repeat multiplication. For 38*25c28e83SPiotr Jasiukajtis * 10**22 is the largest representable integer. 39*25c28e83SPiotr Jasiukajtis */ 40*25c28e83SPiotr Jasiukajtis /* INDENT ON */ 41*25c28e83SPiotr Jasiukajtis 42*25c28e83SPiotr Jasiukajtis #include "libm.h" 43*25c28e83SPiotr Jasiukajtis 44*25c28e83SPiotr Jasiukajtis static const double C[] = { 45*25c28e83SPiotr Jasiukajtis 3.3219280948736234787, /* log(10)/log(2) */ 46*25c28e83SPiotr Jasiukajtis 2.3025850929940456840, /* log(10) */ 47*25c28e83SPiotr Jasiukajtis 3.0102999565860955045E-1, /* log(2)/log(10) high */ 48*25c28e83SPiotr Jasiukajtis 5.3716447674669983622E-12, /* log(2)/log(10) low */ 49*25c28e83SPiotr Jasiukajtis 0.0, 50*25c28e83SPiotr Jasiukajtis 0.5, 51*25c28e83SPiotr Jasiukajtis 1.0, 52*25c28e83SPiotr Jasiukajtis 10.0, 53*25c28e83SPiotr Jasiukajtis 1.0e300, 54*25c28e83SPiotr Jasiukajtis 1.0e-300, 55*25c28e83SPiotr Jasiukajtis }; 56*25c28e83SPiotr Jasiukajtis 57*25c28e83SPiotr Jasiukajtis #define lg10 C[0] 58*25c28e83SPiotr Jasiukajtis #define ln10 C[1] 59*25c28e83SPiotr Jasiukajtis #define logt2hi C[2] 60*25c28e83SPiotr Jasiukajtis #define logt2lo C[3] 61*25c28e83SPiotr Jasiukajtis #define zero C[4] 62*25c28e83SPiotr Jasiukajtis #define half C[5] 63*25c28e83SPiotr Jasiukajtis #define one C[6] 64*25c28e83SPiotr Jasiukajtis #define ten C[7] 65*25c28e83SPiotr Jasiukajtis #define huge C[8] 66*25c28e83SPiotr Jasiukajtis #define tiny C[9] 67*25c28e83SPiotr Jasiukajtis 68*25c28e83SPiotr Jasiukajtis double 69*25c28e83SPiotr Jasiukajtis exp10(double x) { 70*25c28e83SPiotr Jasiukajtis double t, pt; 71*25c28e83SPiotr Jasiukajtis int ix, hx, k; 72*25c28e83SPiotr Jasiukajtis 73*25c28e83SPiotr Jasiukajtis ix = ((int *)&x)[HIWORD]; 74*25c28e83SPiotr Jasiukajtis hx = ix & ~0x80000000; 75*25c28e83SPiotr Jasiukajtis 76*25c28e83SPiotr Jasiukajtis if (hx >= 0x4074a000) { /* |x| >= 330 or x is nan */ 77*25c28e83SPiotr Jasiukajtis if (hx >= 0x7ff00000) { /* x is inf or nan */ 78*25c28e83SPiotr Jasiukajtis if (ix == 0xfff00000 && ((int *)&x)[LOWORD] == 0) 79*25c28e83SPiotr Jasiukajtis return (zero); 80*25c28e83SPiotr Jasiukajtis return (x * x); 81*25c28e83SPiotr Jasiukajtis } 82*25c28e83SPiotr Jasiukajtis t = (ix < 0)? tiny : huge; 83*25c28e83SPiotr Jasiukajtis return (t * t); 84*25c28e83SPiotr Jasiukajtis } 85*25c28e83SPiotr Jasiukajtis 86*25c28e83SPiotr Jasiukajtis if (hx < 0x3c000000) 87*25c28e83SPiotr Jasiukajtis return (one + x); 88*25c28e83SPiotr Jasiukajtis 89*25c28e83SPiotr Jasiukajtis k = (int)x; 90*25c28e83SPiotr Jasiukajtis if (0 <= k && k < 23 && (double)k == x) { 91*25c28e83SPiotr Jasiukajtis /* x is a small positive integer */ 92*25c28e83SPiotr Jasiukajtis t = one; 93*25c28e83SPiotr Jasiukajtis pt = ten; 94*25c28e83SPiotr Jasiukajtis if (k & 1) 95*25c28e83SPiotr Jasiukajtis t = ten; 96*25c28e83SPiotr Jasiukajtis k >>= 1; 97*25c28e83SPiotr Jasiukajtis while (k) { 98*25c28e83SPiotr Jasiukajtis pt *= pt; 99*25c28e83SPiotr Jasiukajtis if (k & 1) 100*25c28e83SPiotr Jasiukajtis t *= pt; 101*25c28e83SPiotr Jasiukajtis k >>= 1; 102*25c28e83SPiotr Jasiukajtis } 103*25c28e83SPiotr Jasiukajtis return (t); 104*25c28e83SPiotr Jasiukajtis } 105*25c28e83SPiotr Jasiukajtis t = x * lg10; 106*25c28e83SPiotr Jasiukajtis k = (int)((ix < 0)? t - half : t + half); 107*25c28e83SPiotr Jasiukajtis return (scalbn(exp(ln10 * ((x - k * logt2hi) - k * logt2lo)), k)); 108*25c28e83SPiotr Jasiukajtis } 109