#include "SDL_internal.h"1/*2* ====================================================3* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.4*5* Developed at SunPro, a Sun Microsystems, Inc. business.6* Permission to use, copy, modify, and distribute this7* software is freely granted, provided that this notice8* is preserved.9* ====================================================10*/1112#if defined(_MSC_VER) /* Handle Microsoft VC++ compiler specifics. */13/* C4723: potential divide by zero. */14#pragma warning ( disable : 4723 )15#endif1617/* __ieee754_log10(x)18* Return the base 10 logarithm of x19*20* Method :21* Let log10_2hi = leading 40 bits of log10(2) and22* log10_2lo = log10(2) - log10_2hi,23* ivln10 = 1/log(10) rounded.24* Then25* n = ilogb(x),26* if(n<0) n = n+1;27* x = scalbn(x,-n);28* log10(x) := n*log10_2hi + (n*log10_2lo + ivln10*log(x))29*30* Note 1:31* To guarantee log10(10**n)=n, where 10**n is normal, the rounding32* mode must set to Round-to-Nearest.33* Note 2:34* [1/log(10)] rounded to 53 bits has error .198 ulps;35* log10 is monotonic at all binary break points.36*37* Special cases:38* log10(x) is NaN with signal if x < 0;39* log10(+INF) is +INF with no signal; log10(0) is -INF with signal;40* log10(NaN) is that NaN with no signal;41* log10(10**N) = N for N=0,1,...,22.42*43* Constants:44* The hexadecimal values are the intended ones for the following constants.45* The decimal values may be used, provided that the compiler will convert46* from decimal to binary accurately enough to produce the hexadecimal values47* shown.48*/4950#include "math_libm.h"51#include "math_private.h"5253static const double54two54 = 1.80143985094819840000e+16, /* 0x43500000, 0x00000000 */55ivln10 = 4.34294481903251816668e-01, /* 0x3FDBCB7B, 0x1526E50E */56log10_2hi = 3.01029995663611771306e-01, /* 0x3FD34413, 0x509F6000 */57log10_2lo = 3.69423907715893078616e-13; /* 0x3D59FEF3, 0x11F12B36 */5859static const double zero = 0.0;6061double attribute_hidden __ieee754_log10(double x)62{63double y,z;64int32_t i,k,hx;65u_int32_t lx;6667EXTRACT_WORDS(hx,lx,x);6869k=0;70if (hx < 0x00100000) { /* x < 2**-1022 */71if (((hx&0x7fffffff)|lx)==0)72return -two54/zero; /* log(+-0)=-inf */73if (hx<0) return (x-x)/zero; /* log(-#) = NaN */74k -= 54; x *= two54; /* subnormal number, scale up x */75GET_HIGH_WORD(hx,x);76}77if (hx >= 0x7ff00000) return x+x;78k += (hx>>20)-1023;79i = ((u_int32_t)k&0x80000000)>>31;80hx = (hx&0x000fffff)|((0x3ff-i)<<20);81y = (double)(k+i);82SET_HIGH_WORD(x,hx);83z = y*log10_2lo + ivln10*__ieee754_log(x);84return z+y*log10_2hi;85}8687/*88* wrapper log10(X)89*/90#ifndef _IEEE_LIBM91double log10(double x)92{93double z = __ieee754_log10(x);94if (_LIB_VERSION == _IEEE_ || isnan(x))95return z;96if (x <= 0.0) {97if(x == 0.0)98return __kernel_standard(x, x, 18); /* log10(0) */99return __kernel_standard(x, x, 19); /* log10(x<0) */100}101return z;102}103#else104strong_alias(__ieee754_log10, log10)105#endif106libm_hidden_def(log10)107108109