/* 
 * mycos1 and mysin1: ordinary taylor expansions of sine and cosine.
 * End the expansion when a term is smaller than 10e-18 of the accumulated series
 * Note: pi is given in 21 decimals in myfunc.h
*/

#include "myfunc.h"
#include <stdio.h>
#include <math.h>

double mycos1(double x) {
	int i = x/PI/2;
	double sum=0;
	double term=1;
        int fac=0;
	x= x -2*PI*i;
	if (x>PI) x -=2*PI;
	if (x<-PI) x+=2*PI;
/* Now, x is between plus and minus PI. Obviously, one could even make x smaller than PI/4
by using cosine and sine identities */
	while (fabs(term) > fabs(sum)*1e-18) {
		sum+= term;
		term*= x/++fac;
		term*= -x/++fac;
	}
	return sum+term;
}

double mysin1(double x) {
	int i = x/PI/2;
	double sum=0;
	double term;
        int fac=1;
	x= x -2*PI*i;
	if (x>PI) x -= 2*PI;
	if (x<-PI) x+=2*PI;
	if (x>PI/2) x=PI-x;
	if (x<-PI/2) x=-PI-x;
/* x is betweeen +PI/2 and -PI/2 now */
	term=x;
	while (fabs(term) > fabs(sum)*1e-18) { 
		sum+=term;
		term*=x/++fac;
		term*=-x/++fac;
	}
	return sum+term;
}

/*
 * mycos2, mysin2: float equivalent of the routines mycos1 and mysin1
 * 
*/

float mycos2(float x) {
	int i = x/2/PI;
	float sum=0;
	float term=1;
        int fac=0;
	x= x -2*PI*i;
	if (x>PI) x -=2*PI;
	if (x<-PI) x+=2*PI;
	while (fabs(term) > fabs(sum)*1e-18) {
		sum+= term;
		term*=x/++fac;
		term*=-x/++fac;
	}
	return sum+term;
}

float mysin2(float x) {
	int i = x/2/PI;
	float sum=0;
	float term;
        int fac=1;
	x= x -2*PI*i;
	if (x>PI) x -= 2*PI;
	if (x<-PI) x+=2*PI;
	if (x>PI/2) x=PI-x;
	if (x<-PI/2) x=-PI-x;
	term=x;
	while (fabs(term) > fabs(sum)*1e-18) {
		sum+=term;
		term*=x/++fac;
		term*=-x/++fac;
	}
	return sum+term;
}

/*
  mycos3, mysin3 : now each term is calculated with pow(x,y). and a call to fac()
*/

inline double fac(int i) {
	double faculty=1.0;
	while (i>1) faculty*=i--;
	return faculty;
}

double mycos3(double x) {
	int i = x/PI/2;
	double sum=0;
	double term=1;
	x= x -2*PI*i;
	if (x>PI) x -=2*PI;
	if (x<-PI) x+=2*PI;
/* Now, x is between plus and minus PI. Obviously, one could gain even more by using cosine and
 sine identities. 
*/
        i=0;
	while (fabs(term) > fabs(sum)*1e-18) {
		i++;
		sum+= term;
		term=pow(x,2.0*i)/fac(2*i);
		if (i%2) term=-term;
        }
	return sum+term;
}

double mysin3(double x) {
	int i = x/PI/2;
	double sum=0;
	double term;
	x= x -2*PI*i;
	if (x>PI) x -=2*PI;
	if (x<-PI) x+=2*PI;
/* Now, x is between plus and minus PI. Obviously, one could gain even more by using cosine and
 sine identities. 
*/
        i=1;
	term=x;
	while (fabs(term) > fabs(sum)*1e-18) {
		sum+= term;
		i++;
		term=pow(x,2.0*i-1.0)/fac(2*i-1);
		if (!(i%2)) term=-term;
        }
	return sum+term;
}


void powercount() {
/* used to see the time needed by a call to pow() */
	int i;
	double x=0;
	for (i=0; i<10000000; i++) {
		x+=pow(2.3, i/40000.0);
	}
	fprintf(stdout,"powercount : %e\n",x);
}