mirror of
https://github.com/DrBeef/Raze.git
synced 2024-11-17 01:41:30 +00:00
275 lines
6.4 KiB
C
275 lines
6.4 KiB
C
|
/* log10.c
|
|||
|
|
*
|
|||
|
|
* Common logarithm
|
|||
|
|
*
|
|||
|
|
*
|
|||
|
|
*
|
|||
|
|
* SYNOPSIS:
|
|||
|
|
*
|
|||
|
|
* double x, y, log10();
|
|||
|
|
*
|
|||
|
|
* y = log10( x );
|
|||
|
|
*
|
|||
|
|
*
|
|||
|
|
*
|
|||
|
|
* DESCRIPTION:
|
|||
|
|
*
|
|||
|
|
* Returns logarithm to the base 10 of x.
|
|||
|
|
*
|
|||
|
|
* The argument is separated into its exponent and fractional
|
|||
|
|
* parts. The logarithm of the fraction is approximated by
|
|||
|
|
*
|
|||
|
|
* log(1+x) = x - 0.5 x**2 + x**3 P(x)/Q(x).
|
|||
|
|
*
|
|||
|
|
*
|
|||
|
|
*
|
|||
|
|
* ACCURACY:
|
|||
|
|
*
|
|||
|
|
* Relative error:
|
|||
|
|
* arithmetic domain # trials peak rms
|
|||
|
|
* IEEE 0.5, 2.0 30000 1.5e-16 5.0e-17
|
|||
|
|
* IEEE 0, MAXNUM 30000 1.4e-16 4.8e-17
|
|||
|
|
* DEC 1, MAXNUM 50000 2.5e-17 6.0e-18
|
|||
|
|
*
|
|||
|
|
* In the tests over the interval [1, MAXNUM], the logarithms
|
|||
|
|
* of the random arguments were uniformly distributed over
|
|||
|
|
* [0, MAXLOG].
|
|||
|
|
*
|
|||
|
|
* ERROR MESSAGES:
|
|||
|
|
*
|
|||
|
|
* log10 singularity: x = 0; returns -INFINITY
|
|||
|
|
* log10 domain: x < 0; returns NAN
|
|||
|
|
*/
|
|||
|
|
|
|||
|
|
/*
|
|||
|
|
Cephes Math Library Release 2.8: June, 2000
|
|||
|
|
Copyright 1984, 1995, 2000 by Stephen L. Moshier
|
|||
|
|
|
|||
|
|
Redistribution and use in source and binary forms, with or without
|
|||
|
|
modification, are permitted provided that the following conditions are met:
|
|||
|
|
|
|||
|
|
1. Redistributions of source code must retain the above copyright notice,
|
|||
|
|
this list of conditions and the following disclaimer.
|
|||
|
|
2. Redistributions in binary form must reproduce the above copyright
|
|||
|
|
notice, this list of conditions and the following disclaimer in the
|
|||
|
|
documentation and/or other materials provided with the distribution.
|
|||
|
|
3. Neither the name of the <ORGANIZATION> nor the names of its
|
|||
|
|
contributors may be used to endorse or promote products derived from
|
|||
|
|
this software without specific prior written permission.
|
|||
|
|
|
|||
|
|
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
|
|||
|
|
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
|
|||
|
|
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
|
|||
|
|
ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
|
|||
|
|
LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
|
|||
|
|
CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
|
|||
|
|
SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
|
|||
|
|
INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
|
|||
|
|
CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
|
|||
|
|
ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
|
|||
|
|
POSSIBILITY OF SUCH DAMAGE.
|
|||
|
|
*/
|
|||
|
|
|
|||
|
|
#include "mconf.h"
|
|||
|
|
static char fname[] = {"log10"};
|
|||
|
|
|
|||
|
|
/* Coefficients for log(1+x) = x - x**2/2 + x**3 P(x)/Q(x)
|
|||
|
|
* 1/sqrt(2) <= x < sqrt(2)
|
|||
|
|
*/
|
|||
|
|
#ifdef UNK
|
|||
|
|
static double P[] = {
|
|||
|
|
4.58482948458143443514E-5,
|
|||
|
|
4.98531067254050724270E-1,
|
|||
|
|
6.56312093769992875930E0,
|
|||
|
|
2.97877425097986925891E1,
|
|||
|
|
6.06127134467767258030E1,
|
|||
|
|
5.67349287391754285487E1,
|
|||
|
|
1.98892446572874072159E1
|
|||
|
|
};
|
|||
|
|
static double Q[] = {
|
|||
|
|
/* 1.00000000000000000000E0, */
|
|||
|
|
1.50314182634250003249E1,
|
|||
|
|
8.27410449222435217021E1,
|
|||
|
|
2.20664384982121929218E2,
|
|||
|
|
3.07254189979530058263E2,
|
|||
|
|
2.14955586696422947765E2,
|
|||
|
|
5.96677339718622216300E1
|
|||
|
|
};
|
|||
|
|
#endif
|
|||
|
|
|
|||
|
|
#ifdef DEC
|
|||
|
|
static unsigned short P[] = {
|
|||
|
|
0034500,0046473,0051374,0135174,
|
|||
|
|
0037777,0037566,0145712,0150321,
|
|||
|
|
0040722,0002426,0031543,0123107,
|
|||
|
|
0041356,0046513,0170752,0004346,
|
|||
|
|
0041562,0071553,0023536,0163343,
|
|||
|
|
0041542,0170221,0024316,0114216,
|
|||
|
|
0041237,0016454,0046611,0104602
|
|||
|
|
};
|
|||
|
|
static unsigned short Q[] = {
|
|||
|
|
/*0040200,0000000,0000000,0000000,*/
|
|||
|
|
0041160,0100260,0067736,0102424,
|
|||
|
|
0041645,0075552,0036563,0147072,
|
|||
|
|
0042134,0125025,0021132,0025320,
|
|||
|
|
0042231,0120211,0046030,0103271,
|
|||
|
|
0042126,0172241,0052151,0120426,
|
|||
|
|
0041556,0125702,0072116,0047103
|
|||
|
|
};
|
|||
|
|
#endif
|
|||
|
|
|
|||
|
|
#ifdef IBMPC
|
|||
|
|
static unsigned short P[] = {
|
|||
|
|
0x974f,0x6a5f,0x09a7,0x3f08,
|
|||
|
|
0x5a1a,0xd979,0xe7ee,0x3fdf,
|
|||
|
|
0x74c9,0xc66c,0x40a2,0x401a,
|
|||
|
|
0x411d,0x7e3d,0xc9a9,0x403d,
|
|||
|
|
0xdcdc,0x64eb,0x4e6d,0x404e,
|
|||
|
|
0xd312,0x2519,0x5e12,0x404c,
|
|||
|
|
0x3130,0x89b1,0xe3a5,0x4033
|
|||
|
|
};
|
|||
|
|
static unsigned short Q[] = {
|
|||
|
|
/*0x0000,0x0000,0x0000,0x3ff0,*/
|
|||
|
|
0xd0a2,0x0dfb,0x1016,0x402e,
|
|||
|
|
0x79c7,0x47ae,0xaf6d,0x4054,
|
|||
|
|
0x455a,0xa44b,0x9542,0x406b,
|
|||
|
|
0x10d7,0x2983,0x3411,0x4073,
|
|||
|
|
0x3423,0x2a8d,0xde94,0x406a,
|
|||
|
|
0xc9c8,0x4e89,0xd578,0x404d
|
|||
|
|
};
|
|||
|
|
#endif
|
|||
|
|
|
|||
|
|
#ifdef MIEEE
|
|||
|
|
static unsigned short P[] = {
|
|||
|
|
0x3f08,0x09a7,0x6a5f,0x974f,
|
|||
|
|
0x3fdf,0xe7ee,0xd979,0x5a1a,
|
|||
|
|
0x401a,0x40a2,0xc66c,0x74c9,
|
|||
|
|
0x403d,0xc9a9,0x7e3d,0x411d,
|
|||
|
|
0x404e,0x4e6d,0x64eb,0xdcdc,
|
|||
|
|
0x404c,0x5e12,0x2519,0xd312,
|
|||
|
|
0x4033,0xe3a5,0x89b1,0x3130
|
|||
|
|
};
|
|||
|
|
static unsigned short Q[] = {
|
|||
|
|
0x402e,0x1016,0x0dfb,0xd0a2,
|
|||
|
|
0x4054,0xaf6d,0x47ae,0x79c7,
|
|||
|
|
0x406b,0x9542,0xa44b,0x455a,
|
|||
|
|
0x4073,0x3411,0x2983,0x10d7,
|
|||
|
|
0x406a,0xde94,0x2a8d,0x3423,
|
|||
|
|
0x404d,0xd578,0x4e89,0xc9c8
|
|||
|
|
};
|
|||
|
|
#endif
|
|||
|
|
|
|||
|
|
#define SQRTH 0.70710678118654752440
|
|||
|
|
#define L102A 3.0078125E-1
|
|||
|
|
#define L102B 2.48745663981195213739E-4
|
|||
|
|
#define L10EA 4.3359375E-1
|
|||
|
|
#define L10EB 7.00731903251827651129E-4
|
|||
|
|
|
|||
|
|
#ifdef ANSIPROT
|
|||
|
|
extern double frexp ( double, int * );
|
|||
|
|
extern double ldexp ( double, int );
|
|||
|
|
extern double polevl ( double, void *, int );
|
|||
|
|
extern double p1evl ( double, void *, int );
|
|||
|
|
extern int isnan ( double );
|
|||
|
|
extern int isfinite ( double );
|
|||
|
|
#else
|
|||
|
|
double frexp(), ldexp(), polevl(), p1evl();
|
|||
|
|
int isnan(), isfinite();
|
|||
|
|
#endif
|
|||
|
|
extern double LOGE2, SQRT2, INFINITY, NAN;
|
|||
|
|
|
|||
|
|
double c_log10(x)
|
|||
|
|
double x;
|
|||
|
|
{
|
|||
|
|
VOLATILE double z;
|
|||
|
|
double y;
|
|||
|
|
#ifdef DEC
|
|||
|
|
short *q;
|
|||
|
|
#endif
|
|||
|
|
int e;
|
|||
|
|
|
|||
|
|
#ifdef NANS
|
|||
|
|
if( isnan(x) )
|
|||
|
|
return(x);
|
|||
|
|
#endif
|
|||
|
|
#ifdef INFINITIES
|
|||
|
|
if( x == INFINITY )
|
|||
|
|
return(x);
|
|||
|
|
#endif
|
|||
|
|
/* Test for domain */
|
|||
|
|
if( x <= 0.0 )
|
|||
|
|
{
|
|||
|
|
if( x == 0.0 )
|
|||
|
|
{
|
|||
|
|
mtherr( fname, SING );
|
|||
|
|
return( -INFINITY );
|
|||
|
|
}
|
|||
|
|
else
|
|||
|
|
{
|
|||
|
|
mtherr( fname, DOMAIN );
|
|||
|
|
return( NAN );
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
/* separate mantissa from exponent */
|
|||
|
|
|
|||
|
|
#ifdef DEC
|
|||
|
|
q = (short *)&x;
|
|||
|
|
e = *q; /* short containing exponent */
|
|||
|
|
e = ((e >> 7) & 0377) - 0200; /* the exponent */
|
|||
|
|
*q &= 0177; /* strip exponent from x */
|
|||
|
|
*q |= 040000; /* x now between 0.5 and 1 */
|
|||
|
|
#endif
|
|||
|
|
|
|||
|
|
#ifdef IBMPC
|
|||
|
|
x = frexp( x, &e );
|
|||
|
|
/*
|
|||
|
|
q = (short *)&x;
|
|||
|
|
q += 3;
|
|||
|
|
e = *q;
|
|||
|
|
e = ((e >> 4) & 0x0fff) - 0x3fe;
|
|||
|
|
*q &= 0x0f;
|
|||
|
|
*q |= 0x3fe0;
|
|||
|
|
*/
|
|||
|
|
#endif
|
|||
|
|
|
|||
|
|
/* Equivalent C language standard library function: */
|
|||
|
|
#ifdef UNK
|
|||
|
|
x = frexp( x, &e );
|
|||
|
|
#endif
|
|||
|
|
|
|||
|
|
#ifdef MIEEE
|
|||
|
|
x = frexp( x, &e );
|
|||
|
|
#endif
|
|||
|
|
|
|||
|
|
/* logarithm using log(1+x) = x - .5x**2 + x**3 P(x)/Q(x) */
|
|||
|
|
|
|||
|
|
if( x < SQRTH )
|
|||
|
|
{
|
|||
|
|
e -= 1;
|
|||
|
|
x = ldexp( x, 1 ) - 1.0; /* 2x - 1 */
|
|||
|
|
}
|
|||
|
|
else
|
|||
|
|
{
|
|||
|
|
x = x - 1.0;
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
|
|||
|
|
/* rational form */
|
|||
|
|
z = x*x;
|
|||
|
|
y = x * ( z * polevl( x, P, 6 ) / p1evl( x, Q, 6 ) );
|
|||
|
|
y = y - ldexp( z, -1 ); /* y - 0.5 * x**2 */
|
|||
|
|
|
|||
|
|
/* multiply log of fraction by log10(e)
|
|||
|
|
* and base 2 exponent by log10(2)
|
|||
|
|
*/
|
|||
|
|
z = (x + y) * L10EB; /* accumulate terms in order of size */
|
|||
|
|
z += y * L10EA;
|
|||
|
|
z += x * L10EA;
|
|||
|
|
z += e * L102B;
|
|||
|
|
z += e * L102A;
|
|||
|
|
|
|||
|
|
|
|||
|
|
return( z );
|
|||
|
|
}
|