tanf.h   [plain text]

/*	This is an implementation of tanf.  It is written in standard C except that
	float must have an IEEE 754 single-precision implementation, and its
	exponent field is accessed via a union and a struct containing bit fields.
	It should be good enough to serve as the libm tanf with tolerable
	performance.  (The exponent field access could be replaced with ilogbf to
	make this a standard C routine, but additional branches would be required
	and performance would be impaired.)

	When this file is compiled, a definition must be supplied for the
	preprocessor symbol Name.  This both becomes the name of the routine and
	specifies the function to perform.  The only supported name is:

		tanf	Tangent function.
*/


// Define constants for preprocessor conditionals.
#define	Tangent	1

// Temporarily define symbols so we can use preprocessor tests.
#define	tanf			Tangent

#if Name == Tangent
	#define	Function	Tangent
#else
	#error "Unknown function."
#endif

// Remove temporary definitions.
#undef	tanf

#define	ExponentFold	1


// Include math.h to ensure we match the declarations of tanf.
#include <math.h>


// Describe the destination floating-point type.
#define	SignificandBits	24
#define	ExponentBits	8


/*	float Name(float x).

	Notes:

		Citations in parentheses below indicate the source of a requirement.

		"C" stands for ISO/IEC 9899:TC2.

		The Open Group specification (IEEE Std 1003.1, 2004 edition) adds no
		requirements since it defers to C and requires errno behavior only if
		we choose to support it by arranging for "math_errhandling &
		MATH_ERRNO" to be non-zero, which we do not.

	Return value:

		For tangent of +/- zero, return zero with same sign (C F.9 12 and
		F.9.1.7).

		For +/- infinity, return NaN (C F.9.1.7).

		For a NaN, return the same NaN (C F.9 11 and 13).

		Otherwise:

			If the rounding mode is round-to-nearest, return tangent(x)
			faithfully rounded.
		  
			Not implemented:  In other rounding modes, return tangent(x)
			possibly with slightly worse error, not necessarily honoring the
			rounding mode (Ali Sazegari narrowing C F.9 10).

	Exceptions:

		Raise underflow for a denormal result (C F.9 7 and Draft Standard for
		Floating-Point Arithmetic P754 Draft 1.2.5 9.5).  If the input is the
		smallest normal, underflow may or may not be raised.  This is stricter
		than the older 754 standard.

		May or may not raise inexact, even if the result is exact (C F.9 8).

		Raise invalid if the input is a signalling NaN (C 5.2.4.2.2 3, in spite
		of C 4.2.1)  or an infinity (C F.9.1.7) but not if the input is a quiet
		NaN (C F.9 11).

		May not raise exceptions otherwise (C F.9 9).

	Properties:

		Monotonic.  Not yet proven.
*/
float Name(float x)
{
	if (isnan(x))
		return x;

	#if defined __STDC__ && 199901L <= __STDC_VERSION__ && !defined __GNUC__
		// GCC does not currently support FP_CONTRACT.  Maybe someday.
		#pragma STDC FP_CONTRACT OFF
	#endif

	typedef struct { double m0, m1; } Double2;

	/* This table was generated by tanfMakeTable.c. */

	static const Double2 ReductionTable[] =
	{
		{ 0x1.45f306ep-1, -0x1.b1bbead603d8bp-32 },	// Exponent -126.
		{ 0x1.45f306ep-1, -0x1.b1bbead603d8bp-32 },	// Exponent -125.
		{ 0x1.45f306ep-1, -0x1.b1bbead603d8bp-32 },	// Exponent -123.
		{ 0x1.45f306ep-1, -0x1.b1bbead603d8bp-32 },	// Exponent -121.
		{ 0x1.45f306ep-1, -0x1.b1bbead603d8bp-32 },	// Exponent -119.
		{ 0x1.45f306ep-1, -0x1.b1bbead603d8bp-32 },	// Exponent -117.
		{ 0x1.45f306ep-1, -0x1.b1bbead603d8bp-32 },	// Exponent -115.
		{ 0x1.45f306ep-1, -0x1.b1bbead603d8bp-32 },	// Exponent -113.
		{ 0x1.45f306ep-1, -0x1.b1bbead603d8bp-32 },	// Exponent -111.
		{ 0x1.45f306ep-1, -0x1.b1bbead603d8bp-32 },	// Exponent -109.
		{ 0x1.45f306ep-1, -0x1.b1bbead603d8bp-32 },	// Exponent -107.
		{ 0x1.45f306ep-1, -0x1.b1bbead603d8bp-32 },	// Exponent -105.
		{ 0x1.45f306ep-1, -0x1.b1bbead603d8bp-32 },	// Exponent -103.
		{ 0x1.45f306ep-1, -0x1.b1bbead603d8bp-32 },	// Exponent -101.
		{ 0x1.45f306ep-1, -0x1.b1bbead603d8bp-32 },	// Exponent -99.
		{ 0x1.45f306ep-1, -0x1.b1bbead603d8bp-32 },	// Exponent -97.
		{ 0x1.45f306ep-1, -0x1.b1bbead603d8bp-32 },	// Exponent -95.
		{ 0x1.45f306ep-1, -0x1.b1bbead603d8bp-32 },	// Exponent -93.
		{ 0x1.45f306ep-1, -0x1.b1bbead603d8bp-32 },	// Exponent -91.
		{ 0x1.45f306ep-1, -0x1.b1bbead603d8bp-32 },	// Exponent -89.
		{ 0x1.45f306ep-1, -0x1.b1bbead603d8bp-32 },	// Exponent -87.
		{ 0x1.45f306ep-1, -0x1.b1bbead603d8bp-32 },	// Exponent -85.
		{ 0x1.45f306ep-1, -0x1.b1bbead603d8bp-32 },	// Exponent -83.
		{ 0x1.45f306ep-1, -0x1.b1bbead603d8bp-32 },	// Exponent -81.
		{ 0x1.45f306ep-1, -0x1.b1bbead603d8bp-32 },	// Exponent -79.
		{ 0x1.45f306ep-1, -0x1.b1bbead603d8bp-32 },	// Exponent -77.
		{ 0x1.45f306ep-1, -0x1.b1bbead603d8bp-32 },	// Exponent -75.
		{ 0x1.45f306ep-1, -0x1.b1bbead603d8bp-32 },	// Exponent -73.
		{ 0x1.45f306ep-1, -0x1.b1bbead603d8bp-32 },	// Exponent -71.
		{ 0x1.45f306ep-1, -0x1.b1bbead603d8bp-32 },	// Exponent -69.
		{ 0x1.45f306ep-1, -0x1.b1bbead603d8bp-32 },	// Exponent -67.
		{ 0x1.45f306ep-1, -0x1.b1bbead603d8bp-32 },	// Exponent -65.
		{ 0x1.45f306ep-1, -0x1.b1bbead603d8bp-32 },	// Exponent -63.
		{ 0x1.45f306ep-1, -0x1.b1bbead603d8bp-32 },	// Exponent -61.
		{ 0x1.45f306ep-1, -0x1.b1bbead603d8bp-32 },	// Exponent -59.
		{ 0x1.45f306ep-1, -0x1.b1bbead603d8bp-32 },	// Exponent -57.
		{ 0x1.45f306ep-1, -0x1.b1bbead603d8bp-32 },	// Exponent -55.
		{ 0x1.45f306ep-1, -0x1.b1bbead603d8bp-32 },	// Exponent -53.
		{ 0x1.45f306ep-1, -0x1.b1bbead603d8bp-32 },	// Exponent -51.
		{ 0x1.45f306ep-1, -0x1.b1bbead603d8bp-32 },	// Exponent -49.
		{ 0x1.45f306ep-1, -0x1.b1bbead603d8bp-32 },	// Exponent -47.
		{ 0x1.45f306ep-1, -0x1.b1bbead603d8bp-32 },	// Exponent -45.
		{ 0x1.45f306ep-1, -0x1.b1bbead603d8bp-32 },	// Exponent -43.
		{ 0x1.45f306ep-1, -0x1.b1bbead603d8bp-32 },	// Exponent -41.
		{ 0x1.45f306ep-1, -0x1.b1bbead603d8bp-32 },	// Exponent -39.
		{ 0x1.45f306ep-1, -0x1.b1bbead603d8bp-32 },	// Exponent -37.
		{ 0x1.45f306ep-1, -0x1.b1bbead603d8bp-32 },	// Exponent -35.
		{ 0x1.45f306ep-1, -0x1.b1bbead603d8bp-32 },	// Exponent -33.
		{ 0x1.45f306ep-1, -0x1.b1bbead603d8bp-32 },	// Exponent -31.
		{ 0x1.45f306ep-1, -0x1.b1bbead603d8bp-32 },	// Exponent -29.
		{ 0x1.45f306ep-1, -0x1.b1bbead603d8bp-32 },	// Exponent -27.
		{ 0x1.45f306ep-1, -0x1.b1bbead603d8bp-32 },	// Exponent -25.
		{ 0x1.45f306ep-1, -0x1.b1bbead603d8bp-32 },	// Exponent -23.
		{ 0x1.45f306ep-1, -0x1.b1bbead603d8bp-32 },	// Exponent -21.
		{ 0x1.45f306ep-1, -0x1.b1bbead603d8bp-32 },	// Exponent -19.
		{ 0x1.45f306ep-1, -0x1.b1bbead603d8bp-32 },	// Exponent -17.
		{ 0x1.45f306ep-1, -0x1.b1bbead603d8bp-32 },	// Exponent -15.
		{ 0x1.45f306ep-1, -0x1.b1bbead603d8bp-32 },	// Exponent -13.
		{ 0x1.45f306ep-1, -0x1.b1bbead603d8bp-32 },	// Exponent -11.
		{ 0x1.45f306ep-1, -0x1.b1bbead603d8bp-32 },	// Exponent -9.
		{ 0x1.45f306ep-1, -0x1.b1bbead603d8bp-32 },	// Exponent -7.
		{ 0x1.45f306ep-1, -0x1.b1bbead603d8bp-32 },	// Exponent -5.
		{ 0x1.45f306ep-1, -0x1.b1bbead603d8bp-32 },	// Exponent -3.
		{ 0x1.45f306ep-1, -0x1.b1bbead603d8bp-32 },	// Exponent -1.
		{ 0x1.45f306ep-1, -0x1.b1bbead603d8bp-32 },	// Exponent 1.
		{ 0x1.45f306ep-1, -0x1.b1bbead603d8bp-32 },	// Exponent 3.
		{ 0x1.45f306ep-1, -0x1.b1bbead603d8bp-32 },	// Exponent 5.
		{ 0x1.45f306ep-1, -0x1.b1bbead603d8bp-32 },	// Exponent 7.
		{ 0x1.45f306ep-1, -0x1.b1bbead603d8bp-32 },	// Exponent 9.
		{ 0x1.45f306ep-1, -0x1.b1bbead603d8bp-32 },	// Exponent 11.
		{ 0x1.45f306ep-1, -0x1.b1bbead603d8bp-32 },	// Exponent 13.
		{ 0x1.45f306ep-1, -0x1.b1bbead603d8bp-32 },	// Exponent 15.
		{ 0x1.45f306ep-1, -0x1.b1bbead603d8bp-32 },	// Exponent 17.
		{ 0x1.45f306ep-1, -0x1.b1bbead603d8bp-32 },	// Exponent 19.
		{ 0x1.45f306ep-1, -0x1.b1bbead603d8bp-32 },	// Exponent 21.
		{ 0x1.45f306ep-1, -0x1.b1bbead603d8bp-32 },	// Exponent 23.
		{ 0x1.17cc1b7p-3, 0x1.391054a7f09d6p-34 },	// Exponent 25.
		{ 0x1.7cc1b72p-7, 0x1.c882a53f84ebp-37 },	// Exponent 27.
		{ 0x1.7cc1b72p-7, 0x1.c882a53f84ebp-37 },	// Exponent 29.
		{ 0x1.f306dcap-9, -0x1.bbead603d8a83p-40 },	// Exponent 31.
		{ -0x1.9f246c7p-14, 0x1.054a7f09d5f48p-46 },	// Exponent 33.
		{ -0x1.9f246c7p-14, 0x1.054a7f09d5f48p-46 },	// Exponent 35.
		{ 0x1.836e4e4p-16, 0x1.054a7f09d5f48p-46 },	// Exponent 37.
		{ -0x1.f246c6fp-18, 0x1.529fc2757d1f5p-52 },	// Exponent 39.
		{ 0x1.b727221p-23, -0x1.5ac07b1505c16p-53 },	// Exponent 41.
		{ 0x1.b727221p-23, -0x1.5ac07b1505c16p-53 },	// Exponent 43.
		{ 0x1.b727221p-23, -0x1.5ac07b1505c16p-53 },	// Exponent 45.
		{ -0x1.236377dp-25, -0x1.6b01ec5417056p-55 },	// Exponent 47.
		{ -0x1.1b1bbebp-28, 0x1.4fe13abe8fa9ap-59 },	// Exponent 49.
		{ 0x1.c9c882ap-29, 0x1.4fe13abe8fa9ap-59 },	// Exponent 51.
		{ -0x1.b1bbeadp-32, -0x1.80f62a0b82b2dp-62 },	// Exponent 53.
		{ 0x1.391054ap-34, 0x1.fc2757d1f534ep-64 },	// Exponent 55.
		{ -0x1.8ddf56bp-35, -0x1.ec54170565912p-71 },	// Exponent 57.
		{ 0x1.c882a54p-37, -0x1.ec54170565912p-71 },	// Exponent 59.
		{ -0x1.bbead6p-40, -0x1.ec54170565912p-71 },	// Exponent 61.
		{ 0x1.1054a7fp-42, 0x1.3abe8fa9a6eep-75 },	// Exponent 63.
		{ -0x1.df56b02p-43, 0x1.3abe8fa9a6eep-75 },	// Exponent 65.
		{ 0x1.054a7f1p-46, -0x1.8a82e0acb223fp-76 },	// Exponent 67.
		{ -0x1.f56b01fp-47, 0x1.d5f47d4d37703p-78 },	// Exponent 69.
		{ 0x1.529fc27p-52, 0x1.5f47d4d377037p-82 },	// Exponent 71.
		{ 0x1.529fc27p-52, 0x1.5f47d4d377037p-82 },	// Exponent 73.
		{ -0x1.5ac07b1p-53, -0x1.4170565911f92p-83 },	// Exponent 75.
		{ -0x1.6b01ec5p-55, -0x1.05c1596447e49p-85 },	// Exponent 77.
		{ -0x1.ac07b15p-57, -0x1.70565911f924fp-91 },	// Exponent 79.
		{ 0x1.4fe13acp-59, -0x1.70565911f924fp-91 },	// Exponent 81.
		{ 0x1.3f84ebp-61, -0x1.70565911f924fp-91 },	// Exponent 83.
		{ 0x1.fc2757dp-64, 0x1.f534ddc0db629p-96 },	// Exponent 85.
		{ -0x1.ec5417p-71, -0x1.596447e493ad5p-101 },	// Exponent 87.
		{ -0x1.ec5417p-71, -0x1.596447e493ad5p-101 },	// Exponent 89.
		{ -0x1.ec5417p-71, -0x1.596447e493ad5p-101 },	// Exponent 91.
		{ -0x1.ec5417p-71, -0x1.596447e493ad5p-101 },	// Exponent 93.
		{ 0x1.3abe8fbp-75, -0x1.96447e493ad4dp-105 },	// Exponent 95.
		{ 0x1.3abe8fbp-75, -0x1.96447e493ad4dp-105 },	// Exponent 97.
		{ 0x1.d5f47d5p-78, -0x1.6447e493ad4cep-109 },	// Exponent 99.
		{ -0x1.505c159p-81, -0x1.911f924eb5336p-111 },	// Exponent 101.
		{ -0x1.505c159p-81, -0x1.911f924eb5336p-111 },	// Exponent 103.
		{ -0x1.4170566p-83, 0x1.bb81b6c52b328p-113 },	// Exponent 105.
		{ -0x1.05c1596p-85, -0x1.11f924eb53362p-115 },	// Exponent 107.
		{ -0x1.7056591p-91, -0x1.f924eb53361dep-123 },	// Exponent 109.
		{ -0x1.7056591p-91, -0x1.f924eb53361dep-123 },	// Exponent 111.
		{ -0x1.7056591p-91, -0x1.f924eb53361dep-123 },	// Exponent 113.
		{ -0x1.c159644p-93, -0x1.f924eb53361dep-123 },	// Exponent 115.
		{ 0x1.f534ddcp-96, 0x1.b6c52b3278872p-129 },	// Exponent 117.
		{ -0x1.596447ep-101, -0x1.24eb53361de38p-131 },	// Exponent 119.
		{ -0x1.596447ep-101, -0x1.24eb53361de38p-131 },	// Exponent 121.
		{ -0x1.596447ep-101, -0x1.24eb53361de38p-131 },	// Exponent 123.
		{ -0x1.65911f9p-103, -0x1.275a99b0ef1bfp-134 },	// Exponent 125.
		{ -0x1.96447e5p-105, 0x1.b14acc9e21c82p-135 },	// Exponent 127.
	};

	// Define a structure for accessing the internal components of a float.
	typedef union
	{
		float f;
		struct
		{
			#if defined __BIG_ENDIAN__
				unsigned int sign     : 1;
				unsigned int exponent : ExponentBits;
				unsigned int fraction : SignificandBits-1;
			#else	// defined __BIG_ENDIAN__
				unsigned int fraction : SignificandBits-1;
				unsigned int exponent : ExponentBits;
				unsigned int sign     : 1;
			#endif	// defined __BIG_ENDIAN__
		} s;
	} Float;

	Float X = { x };
	
	/*	Look up 2**ArcBits * (u/p % 1)/u in the reduction table.  See
		tanfMakeTable.c for more information.  The properties we need are:

			The multiplication of residue.m0 by x is exact.

			Let ki be residue.m0 * x rounded to an integer.

			Let xr be residue.m0 * x - ki + residue.m1 * x.

			If ki is even, the tangent of xr is nearly equal to the tangent of
			x.

			If ki is odd, the negative of the tangent of xr is nearly equal to
			the tangent of x.
	*/
	Double2 residue = ReductionTable[X.s.exponent >> ExponentFold];

	double p = x * residue.m0;

	// Round p to an integer in k.
	int ki = (int) lrint(p);
	double k = ki;
		/*	Alternatives to lrint:

				Use nearbyint.  This will not raise the inexact exception.
				However, we already raise inexact in the version-two polynomial
				even when x is zero.

				Use lround.  This will take longer but will round to nearest
				regardless of rounding mode.
		*/

	double f0 = p - k;
	double f1 = x * residue.m1;

	// xr is x reduced to a final residue.
	double xr = f0 + f1;

	// Approximate tangent with a polynomial.

	double x2 = xr*xr;	// Square xr.

	double t = 
			(   ((x2 - 0.34431297545117363561) * x2 + 0.75337415977492632831)
			* ((x2 + 0.91301179133375347365) * x2 + 0.58807104225499649417)
			)
		*
			(   ((x2 - 1.53503935362160698123) * x2 + 0.92783629599735333058)
			* ((   + 5.27860261872543679865) * x2 + 3.82127246397188494330)
			)
		*
			xr;

	// If low bit of ki is 1, return t or -1/t according to low bit of ki.
	return (float) (ki & 1 ? -1/t : t);
}