src/libm/e_log.c
author Ryan C. Gordon <icculus@icculus.org>
Thu, 21 Apr 2016 03:16:44 -0400
changeset 11729 d1ce8396c356
parent 11711 8a982ed61896
permissions -rw-r--r--
Initial shot at a renderer target for Apple's Metal API.

This isn't complete, but is enough to run testsprite2. It's currently
Mac-only; with a little work to figure out how to properly glue in a Metal
layer to a UIView, this will likely work on iOS, too.

This is only wired up to the configure script right now, and disabled by
default. CMake and Xcode still need their bits filled in as appropriate.
slouken@2756
     1
/*
slouken@2756
     2
 * ====================================================
slouken@2756
     3
 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
slouken@2756
     4
 *
slouken@2756
     5
 * Developed at SunPro, a Sun Microsystems, Inc. business.
slouken@2756
     6
 * Permission to use, copy, modify, and distribute this
slouken@2756
     7
 * software is freely granted, provided that this notice
slouken@2756
     8
 * is preserved.
slouken@2756
     9
 * ====================================================
slouken@2756
    10
 */
slouken@2756
    11
slouken@11711
    12
#if defined(_MSC_VER)           /* Handle Microsoft VC++ compiler specifics. */
slouken@11711
    13
/* C4723: potential divide by zero. */
slouken@11711
    14
#pragma warning ( disable : 4723 )
slouken@11711
    15
#endif
slouken@11711
    16
slouken@2756
    17
/* __ieee754_log(x)
slouken@2756
    18
 * Return the logrithm of x
slouken@2756
    19
 *
slouken@2756
    20
 * Method :
slouken@2756
    21
 *   1. Argument Reduction: find k and f such that
slouken@2756
    22
 *			x = 2^k * (1+f),
slouken@2756
    23
 *	   where  sqrt(2)/2 < 1+f < sqrt(2) .
slouken@2756
    24
 *
slouken@2756
    25
 *   2. Approximation of log(1+f).
slouken@2756
    26
 *	Let s = f/(2+f) ; based on log(1+f) = log(1+s) - log(1-s)
slouken@2756
    27
 *		 = 2s + 2/3 s**3 + 2/5 s**5 + .....,
slouken@2756
    28
 *	     	 = 2s + s*R
slouken@2756
    29
 *      We use a special Reme algorithm on [0,0.1716] to generate
slouken@2756
    30
 * 	a polynomial of degree 14 to approximate R The maximum error
slouken@2756
    31
 *	of this polynomial approximation is bounded by 2**-58.45. In
slouken@2756
    32
 *	other words,
slouken@2756
    33
 *		        2      4      6      8      10      12      14
slouken@2756
    34
 *	    R(z) ~ Lg1*s +Lg2*s +Lg3*s +Lg4*s +Lg5*s  +Lg6*s  +Lg7*s
slouken@2756
    35
 *  	(the values of Lg1 to Lg7 are listed in the program)
slouken@2756
    36
 *	and
slouken@2756
    37
 *	    |      2          14          |     -58.45
slouken@2756
    38
 *	    | Lg1*s +...+Lg7*s    -  R(z) | <= 2
slouken@2756
    39
 *	    |                             |
slouken@2756
    40
 *	Note that 2s = f - s*f = f - hfsq + s*hfsq, where hfsq = f*f/2.
slouken@2756
    41
 *	In order to guarantee error in log below 1ulp, we compute log
slouken@2756
    42
 *	by
slouken@2756
    43
 *		log(1+f) = f - s*(f - R)	(if f is not too large)
slouken@2756
    44
 *		log(1+f) = f - (hfsq - s*(hfsq+R)).	(better accuracy)
slouken@2756
    45
 *
slouken@2756
    46
 *	3. Finally,  log(x) = k*ln2 + log(1+f).
slouken@2756
    47
 *			    = k*ln2_hi+(f-(hfsq-(s*(hfsq+R)+k*ln2_lo)))
slouken@2756
    48
 *	   Here ln2 is split into two floating point number:
slouken@2756
    49
 *			ln2_hi + ln2_lo,
slouken@2756
    50
 *	   where n*ln2_hi is always exact for |n| < 2000.
slouken@2756
    51
 *
slouken@2756
    52
 * Special cases:
slouken@2756
    53
 *	log(x) is NaN with signal if x < 0 (including -INF) ;
slouken@2756
    54
 *	log(+INF) is +INF; log(0) is -INF with signal;
slouken@2756
    55
 *	log(NaN) is that NaN with no signal.
slouken@2756
    56
 *
slouken@2756
    57
 * Accuracy:
slouken@2756
    58
 *	according to an error analysis, the error is always less than
slouken@2756
    59
 *	1 ulp (unit in the last place).
slouken@2756
    60
 *
slouken@2756
    61
 * Constants:
slouken@2756
    62
 * The hexadecimal values are the intended ones for the following
slouken@2756
    63
 * constants. The decimal values may be used, provided that the
slouken@2756
    64
 * compiler will convert from decimal to binary accurately enough
slouken@2756
    65
 * to produce the hexadecimal values shown.
slouken@2756
    66
 */
slouken@2756
    67
slouken@6044
    68
#include "math_libm.h"
slouken@2756
    69
#include "math_private.h"
slouken@2756
    70
slouken@2756
    71
static const double
slouken@11683
    72
ln2_hi  =  6.93147180369123816490e-01,	/* 3fe62e42 fee00000 */
slouken@11683
    73
ln2_lo  =  1.90821492927058770002e-10,	/* 3dea39ef 35793c76 */
slouken@11683
    74
two54   =  1.80143985094819840000e+16,  /* 43500000 00000000 */
slouken@11683
    75
Lg1 = 6.666666666666735130e-01,  /* 3FE55555 55555593 */
slouken@11683
    76
Lg2 = 3.999999999940941908e-01,  /* 3FD99999 9997FA04 */
slouken@11683
    77
Lg3 = 2.857142874366239149e-01,  /* 3FD24924 94229359 */
slouken@11683
    78
Lg4 = 2.222219843214978396e-01,  /* 3FCC71C5 1D8E78AF */
slouken@11683
    79
Lg5 = 1.818357216161805012e-01,  /* 3FC74664 96CB03DE */
slouken@11683
    80
Lg6 = 1.531383769920937332e-01,  /* 3FC39A09 D078C69F */
slouken@11683
    81
Lg7 = 1.479819860511658591e-01;  /* 3FC2F112 DF3E5244 */
slouken@11683
    82
slouken@11683
    83
static const double zero   =  0.0;
slouken@11683
    84
slouken@11683
    85
double attribute_hidden __ieee754_log(double x)
slouken@11683
    86
{
slouken@11683
    87
	double hfsq,f,s,z,R,w,t1,t2,dk;
slouken@11683
    88
	int32_t k,hx,i,j;
slouken@11683
    89
	u_int32_t lx;
slouken@11683
    90
slouken@11683
    91
	EXTRACT_WORDS(hx,lx,x);
slouken@11683
    92
slouken@11683
    93
	k=0;
slouken@11683
    94
	if (hx < 0x00100000) {			/* x < 2**-1022  */
slouken@11683
    95
	    if (((hx&0x7fffffff)|lx)==0)
slouken@11683
    96
		return -two54/zero;		/* log(+-0)=-inf */
slouken@11683
    97
	    if (hx<0) return (x-x)/zero;	/* log(-#) = NaN */
slouken@11683
    98
	    k -= 54; x *= two54; /* subnormal number, scale up x */
slouken@11683
    99
	    GET_HIGH_WORD(hx,x);
slouken@11683
   100
	}
slouken@11683
   101
	if (hx >= 0x7ff00000) return x+x;
slouken@11683
   102
	k += (hx>>20)-1023;
slouken@11683
   103
	hx &= 0x000fffff;
slouken@11683
   104
	i = (hx+0x95f64)&0x100000;
slouken@11683
   105
	SET_HIGH_WORD(x,hx|(i^0x3ff00000));	/* normalize x or x/2 */
slouken@11683
   106
	k += (i>>20);
slouken@11683
   107
	f = x-1.0;
slouken@11683
   108
	if((0x000fffff&(2+hx))<3) {	/* |f| < 2**-20 */
slouken@11683
   109
	    if(f==zero) {if(k==0) return zero;  else {dk=(double)k;
slouken@11683
   110
				 return dk*ln2_hi+dk*ln2_lo;}
slouken@11683
   111
	    }
slouken@11683
   112
	    R = f*f*(0.5-0.33333333333333333*f);
slouken@11683
   113
	    if(k==0) return f-R; else {dk=(double)k;
slouken@11683
   114
	    	     return dk*ln2_hi-((R-dk*ln2_lo)-f);}
slouken@11683
   115
	}
slouken@11683
   116
 	s = f/(2.0+f);
slouken@11683
   117
	dk = (double)k;
slouken@11683
   118
	z = s*s;
slouken@11683
   119
	i = hx-0x6147a;
slouken@11683
   120
	w = z*z;
slouken@11683
   121
	j = 0x6b851-hx;
slouken@11683
   122
	t1= w*(Lg2+w*(Lg4+w*Lg6));
slouken@11683
   123
	t2= z*(Lg1+w*(Lg3+w*(Lg5+w*Lg7)));
slouken@11683
   124
	i |= j;
slouken@11683
   125
	R = t2+t1;
slouken@11683
   126
	if(i>0) {
slouken@11683
   127
	    hfsq=0.5*f*f;
slouken@11683
   128
	    if(k==0) return f-(hfsq-s*(hfsq+R)); else
slouken@11683
   129
		     return dk*ln2_hi-((hfsq-(s*(hfsq+R)+dk*ln2_lo))-f);
slouken@11683
   130
	} else {
slouken@11683
   131
	    if(k==0) return f-s*(f-R); else
slouken@11683
   132
		     return dk*ln2_hi-((s*(f-R)-dk*ln2_lo)-f);
slouken@11683
   133
	}
slouken@11683
   134
}
slouken@11683
   135
slouken@11683
   136
/*
slouken@11683
   137
 * wrapper log(x)
slouken@11683
   138
 */
slouken@11683
   139
#ifndef _IEEE_LIBM
slouken@11683
   140
double log(double x)
slouken@11683
   141
{
slouken@11683
   142
	double z = __ieee754_log(x);
slouken@11683
   143
	if (_LIB_VERSION == _IEEE_ || isnan(x) || x > 0.0)
slouken@11683
   144
		return z;
slouken@11683
   145
	if (x == 0.0)
slouken@11683
   146
		return __kernel_standard(x, x, 16); /* log(0) */
slouken@11683
   147
	return __kernel_standard(x, x, 17); /* log(x<0) */
slouken@11683
   148
}
slouken@2756
   149
#else
slouken@11683
   150
strong_alias(__ieee754_log, log)
slouken@2756
   151
#endif
slouken@11683
   152
libm_hidden_def(log)