#include <cmath>
#include <iostream>

// Define a number of constants

// atan(1) == pi/4
// By IEEE 754, atan(1) must return an anser correct to 0.6 ulps
// Multiplication by 2 does not decrease precision

const double PI_BY_2 = 2 * std::atan( 1.0 );
const double PI      = 2 * PI_BY_2;
const double TWO_PI  = 2 * PI;

/*
 * Function Name:  my_sin
 * File:           sin.cpp
 *
 * Douglas Wilhelm Harder
 * 2009-04-20
 * Version 1.0         Initial implementation.
 *
 * Calculuate sin(x) where x is a double-precision floating-point
 * number using Taylor series and trigonometric relations.
 *
 * No pre- or postconditions.
 *
 * Parameters:
 *        double x             A real value assumed to be radians.
 *
 * Returns:
 *        double               The sine of x.
 *
 * Bugs:           Does not work for NaN or +/- infinity.
 * Todo:           Add the numbers in reverse order to reduce error
 *                 due to numeric round off.
 */

double my_sin( double x ) {
	// Divide the interval into five intervals, solving it on
	//
	//        I       II     III         IV            V
	//            .-------.       o---------------.
	// ...--------o       o-------.               o--------...
	// ---+-------+-------+-------+-------+-------+-------+---
	//    pi              pi             3pi             5pi
	//  - --      0       --      pi     ---     2pi     ---
	//    2               2               2               2
	// Region          Formula
	//   I            sin(x) = -sin(-x)
	//   II           Taylor series
	//   III          sin(x) = sin(pi - x)
	//   IV           use sin(x) = -sin(x - pi)
	//   V            sin(x) = sin(x - 2pi floor(x/(2pi)))

	if ( x < 0 ) {
		// Region I
		// For x < 0, use sin(-x) = -sin(x)

		return -sin( -x );
	} else if ( x <= PI_BY_2 ) {
		// Region II
		//               pi
		// For 0 <= x <= --, use the Taylor series
		//               2
		// 
		//              N
		//            -----     2i + 1        2N + 3
		//             \       x            xi
		//  sin(x) =    )     ---------  +  ---------
		//             /      (2i + 1)!     (2N + 3)!
		//            -----
		//            i = 0
		//                              -53        pi
		// The error term is less than 2   for x = -- when N =  11
		//                                         2

		double term = x;
		double result = term;
		double nxx = -x*x;     // each iteration, the term is
		                       // is multiplied by -x^2/(i - 1)/i

		for ( int i = 3; i <= 23; i += 2 ) {
			term *= nxx/static_cast<double>(i*(i - 1));
			result += term;
		}

		return result;
	} else if ( x <= PI ) {
		// Region III
		//     pi
		// For -- < x < pi, use sin(x) = sin(pi - x)
		//     2

		return sin( PI - x );
	} else if ( x <= TWO_PI ) {
		// Region IV
		// For pi < x < 2pi, use sin(x) = -sin(x - pi)

		return -sin( x - PI );
	} else {
		// Region V
		// Otherwise, x > 2pi, and therefore
		// use sin(x) = sin(x - 2pi floor(x/(2pi))) where
		// where 0 <= x - 2pi floor(x/(2pi)) < 2pi

		return sin( x - TWO_PI*std::floor( x/TWO_PI ) );
	}
}

/*
 * int main()
 *
 * Douglas Wilhelm Harder
 * 2009-04-20
 * Version:        1.0         Initial implementation.
 *
 * Unit test function for double my_sin( double x ) by
 * testing sin(n) for n = -100, -99, ..., -1, 0, 1, ..., 99, 100.
 * If the relative error is too large, some diagnostic information
 * is printed.
 *
 * No pre- or postconditions.
 * No parameters.
 *
 * Todo: Does not test x = NaN or x = +/- infinity.           
 */

int main() {
	// The maximum allowable error for our approximation
	const double MAX_ERROR = 5e-14;

	// This array stores the values answers[i] = sin(i) for i = 0, ..., 100.
	// The results were calculated in Maple to 20 decimal digits of precision.

	double answers[101] = {
		 0.0,                       0.84147098480789650665,
		 0.90929742682568169540,    0.14112000805986722210,
		-0.75680249530792825137,   -0.95892427466313846889,
		-0.27941549819892587281,    0.65698659871878909040,
		 0.98935824662338177781,    0.41211848524175656976,
		-0.54402111088936981340,   -0.99999020655070345705,
		-0.53657291800043497167,    0.42016703682664092187,
		 0.99060735569487030788,    0.65028784015711686583,
		-0.28790331666506529478,   -0.96139749187955685726,
		-0.75098724677167610375,    0.14987720966295232975,
		 0.91294525072762765438,    0.83665563853605603187,
		-0.0088513092904038759217, -0.84622040417517063524,
		-0.90557836200662384514,   -0.13235175009777302890,
		 0.76255845047960273752,    0.95637592840450301343,
		 0.27090578830786901999,   -0.66363388421296750215,
		-0.98803162409286178999,   -0.40403764532306500605,
		 0.55142668124169055066,    0.99991186010726714573,
		 0.52908268612002382083,   -0.42818266949615100441,
		-0.99177885344311573684,   -0.64353813335699946069,
		 0.29636857870938531739,    0.96379538628408775326,
		 0.74511316047934878699,   -0.15862266880470898710,
		-0.91652154791563378590,   -0.83177474262859828821,
		 0.017701925105413577808,   0.85090352453411842486,
		 0.90178834764880918503,    0.12357312274522400406,
		-0.76825466132366679904,   -0.95375265275947181836,
		-0.26237485370392878591,    0.67022917584337473449,
		 0.98662759204048529659,    0.39592515018183418150,
		-0.55878904885161624582,   -0.99975517335861983660,
		-0.52155100208691188019,    0.43616475524782495908,
		 0.99287264808453711817,    0.63673800713913788077,
		-0.30481062110221670563,   -0.96611777000839294702,
		-0.73918069664922286728,    0.16735570030280692153,
		 0.92002603819679068335,    0.82682867949010346771,
		-0.026551154023966794464,  -0.85551997897532225900,
		-0.89792768068929126040,   -0.11478481378318722055,
		 0.77389068155788909779,    0.95105465325437463666,
		 0.25382336276203627307,   -0.67677195688730762215,
		-0.98514626046824737085,   -0.38778163540943043773,
		 0.56610763689818032361,    0.99952015858073124387,
		 0.51397845598753521170,   -0.44411266870750836851,
		-0.99388865392337518973,   -0.62988799427445387857,
		 0.31322878243308515263,    0.96836446110018540435,
		 0.73319032007329216636,   -0.17607561994858707696,
		-0.92345844700405980260,   -0.82181783663082254487,
		 0.035398302733660683625,   0.86006940581245322684,
		 0.89399666360055789052,    0.10598751175115685002,
		-0.77946606961580468855,   -0.94828214126994723213,
		-0.24525198546765432522,    0.68326171473612098370,
		 0.98358774543434485761,    0.37960773902752169648,
		-0.57338187199042288495,   -0.99920683418635369443,
		-0.50636564110975879366
	};

	// set the precision of the output to 18 decimal digits
	std::cout.precision( 18 );

	// For each value i = 0, ..., 100, calcluate sin(i).
	// Find the relative error of our approximation by comparing
	// the answer to answers[i] and print some diagnostic information
	// if the relative error is greater than MAX_ERROR.

	for ( int i = -100; i <= 100; ++i ) {
		double x = static_cast<double>( i );
		double approx = my_sin( x );
		double answer = (i >= 0) ? answers[i] : -answers[-i];

		double error;

		error = std::fabs( (approx - answer)/answer );

		// If the relative error is greater than MAX_ERROR, then print
		// diagnoistic information.  For example:
		//
		//     sin(44) = 0.01770192510541358
		//     Approximation:     0.01770192510541529
		//     Standard Library:  0.01770192510541529
		//     Answer (Maple):    0.01770192510541358
		//     Relative Error:    9.682036185662886e-14

		if ( error > MAX_ERROR ) { 
			std::cout << "sin(" << i << ") "                    << std::endl;
			std::cout << "Approximation:     " << approx        << std::endl;
			std::cout << "Standard Library:  " << std::sin( x ) << std::endl;
			std::cout << "Answer (Maple):    " << answer        << std::endl;
			std::cout << "Relative Error:    " << error         << std::endl
			                                                    << std::endl;
		}
	}

	return 0;
}