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 #if defined(ELFOBJ) 31 #pragma weak log10l = __log10l 32 #endif 33 34 /* 35 * log10l(X) 36 * 37 * Method : 38 * Let log10_2hi = leading 98(SPARC)/49(x86) bits of log10(2) and 39 * log10_2lo = log10(2) - log10_2hi, 40 * ivln10 = 1/log(10) rounded. 41 * Then 42 * n = ilogb(x), 43 * if (n<0) n = n+1; 44 * x = scalbn(x,-n); 45 * LOG10(x) := n*log10_2hi + (n*log10_2lo + ivln10*log(x)) 46 * 47 * Note1: 48 * For fear of destroying log10(10**n)=n, the rounding mode is 49 * set to Round-to-Nearest. 50 * 51 * Special cases: 52 * log10(x) is NaN with signal if x < 0; 53 * log10(+INF) is +INF with no signal; log10(0) is -INF with signal; 54 * log10(NaN) is that NaN with no signal; 55 * log10(10**N) = N for N=0,1,...,22. 56 * 57 * Constants: 58 * The hexadecimal values are the intended ones for the following constants. 59 * The decimal values may be used, provided that the compiler will convert 60 * from decimal to binary accurately enough to produce the hexadecimal values 61 * shown. 62 */ 63 64 #include "libm.h" 65 #include "longdouble.h" 66 67 #if defined(__x86) 68 #define __swapRD __swap87RD 69 #endif 70 extern enum fp_direction_type __swapRD(enum fp_direction_type); 71 72 static const long double 73 zero = 0.0L, 74 ivln10 = 4.342944819032518276511289189166050822944e-0001L, 75 one = 1.0L, 76 #if defined(__x86) 77 log10_2hi = 3.010299956639803653501985536422580480576e-01L, 78 log10_2lo = 8.298635403410822349787106337291183585413e-16L; 79 #elif defined(__sparc) 80 log10_2hi = 3.010299956639811952137388947242098603469e-01L, 81 log10_2lo = 2.831664213089468167896664371953210945664e-31L; 82 #else 83 #error Unknown Architecture! 84 #endif 85 86 long double 87 log10l(long double x) { 88 long double y, z; 89 enum fp_direction_type rd; 90 int n; 91 92 if (!finitel(x)) 93 return (x + fabsl(x)); /* x is +-INF or NaN */ 94 else if (x > zero) { 95 n = ilogbl(x); 96 if (n < 0) 97 n += 1; 98 rd = __swapRD(fp_nearest); 99 y = n; 100 x = scalbnl(x, -n); 101 z = y * log10_2lo + ivln10 * logl(x); 102 z += y * log10_2hi; 103 if (rd != fp_nearest) 104 (void) __swapRD(rd); 105 return (z); 106 } else if (x == zero) /* -INF */ 107 return (-one / zero); 108 else /* x <0, return NaN */ 109 return (zero / zero); 110 } 111