log10l.c revision ddc0e0b53c661f6e439e3b7072b3ef353eadb4af
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis/*
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * CDDL HEADER START
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis *
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * The contents of this file are subject to the terms of the
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * Common Development and Distribution License (the "License").
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * You may not use this file except in compliance with the License.
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis *
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * or http://www.opensolaris.org/os/licensing.
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * See the License for the specific language governing permissions
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * and limitations under the License.
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis *
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * When distributing Covered Code, include this CDDL HEADER in each
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * file and include the License file at usr/src/OPENSOLARIS.LICENSE.
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * If applicable, add the following below this CDDL HEADER, with the
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * fields enclosed by brackets "[]" replaced with your own identifying
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * information: Portions Copyright [yyyy] [name of copyright owner]
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis *
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * CDDL HEADER END
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis/*
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * Copyright 2011 Nexenta Systems, Inc. All rights reserved.
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis/*
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * Copyright 2006 Sun Microsystems, Inc. All rights reserved.
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * Use is subject to license terms.
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis
ddc0e0b53c661f6e439e3b7072b3ef353eadb4afRichard Lowe#pragma weak __log10l = log10l
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis/*
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * log10l(X)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis *
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * Method :
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * Let log10_2hi = leading 98(SPARC)/49(x86) bits of log10(2) and
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * log10_2lo = log10(2) - log10_2hi,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * ivln10 = 1/log(10) rounded.
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * Then
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * n = ilogb(x),
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * if (n<0) n = n+1;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * x = scalbn(x,-n);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * LOG10(x) := n*log10_2hi + (n*log10_2lo + ivln10*log(x))
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis *
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * Note1:
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * For fear of destroying log10(10**n)=n, the rounding mode is
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * set to Round-to-Nearest.
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis *
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * Special cases:
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * log10(x) is NaN with signal if x < 0;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * log10(+INF) is +INF with no signal; log10(0) is -INF with signal;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * log10(NaN) is that NaN with no signal;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * log10(10**N) = N for N=0,1,...,22.
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis *
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * Constants:
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * The hexadecimal values are the intended ones for the following constants.
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * The decimal values may be used, provided that the compiler will convert
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * from decimal to binary accurately enough to produce the hexadecimal values
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * shown.
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis#include "libm.h"
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis#include "longdouble.h"
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis#if defined(__x86)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis#define __swapRD __swap87RD
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis#endif
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtisextern enum fp_direction_type __swapRD(enum fp_direction_type);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtisstatic const long double
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis zero = 0.0L,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis ivln10 = 4.342944819032518276511289189166050822944e-0001L,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis one = 1.0L,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis#if defined(__x86)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis log10_2hi = 3.010299956639803653501985536422580480576e-01L,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis log10_2lo = 8.298635403410822349787106337291183585413e-16L;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis#elif defined(__sparc)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis log10_2hi = 3.010299956639811952137388947242098603469e-01L,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis log10_2lo = 2.831664213089468167896664371953210945664e-31L;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis#else
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis#error Unknown Architecture!
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis#endif
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtislong double
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtislog10l(long double x) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis long double y, z;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis enum fp_direction_type rd;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis int n;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis if (!finitel(x))
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis return (x + fabsl(x)); /* x is +-INF or NaN */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis else if (x > zero) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis n = ilogbl(x);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis if (n < 0)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis n += 1;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis rd = __swapRD(fp_nearest);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis y = n;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis x = scalbnl(x, -n);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis z = y * log10_2lo + ivln10 * logl(x);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis z += y * log10_2hi;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis if (rd != fp_nearest)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis (void) __swapRD(rd);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis return (z);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis } else if (x == zero) /* -INF */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis return (-one / zero);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis else /* x <0, return NaN */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis return (zero / zero);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis}