Writing my own math library from scratch

[fruitloop2222]

I was wondering, is the time required worth the experience-gain to write your own math library from scratch? One thing I challenged myself to is to at least make the trigonometry functions, but I'm not even half way there! I discovered that I need a robust- power/exponent function for this, and even this may or may not require (I'm still trying to figure out) a floor function, which may require a less-robust power/exponent function. By robustness, I mean that it either can or cannot handle decimal numbers.

Could some one teach me how a robust power/exponent function works. I thought it would use 1 % exponent, and apparently you can't do modulo with decimal numbers (I tried). So then I decided that I could use the hack for modulo, which requires a floor function and blah blah blah... Math functions are quite complex!

[D80DSL3A]

Doing so could be a great programming exercise for you. Go for it!!
Here's something to start with. This program finds sine() and cosine() using a power series expansion. The SineN2 and CosineN2 functions are more efficient. See notes in the code.
These functions even find the maximum error (easy to find because these are alternating series). I recommend studying power series in depth before trying to write more functions like these but power series expansions exist for all the major math functions.
Study the code and then if you have any questions about it (after studying it) feel free to ask.

// Evaluates power series for sin(x) and cos(x) to n terms in 2 ways
// values entered are not checked so enter carefully x = -PI to +PI
// and n not too high. Keep eye on value of n! vs. DBL_MAX in output
// NOTE: n=12 gives about 10 digit accuracy (series converges rapidly)
// so no need for high n values

#include <iostream>// for use of cin, cout
#include <iomanip>// for use of setprecision
#include <cstdio>// for use of gets(), scanf()
#include <cmath>// for use of sin, cos, pow functions
#include <float.h>// for use of the constant DBL_MAX
using namespace std;

double factorial(double n)// not needed by SineN2, CosineN2 functions
{
	if (n<=1.0)
		return(1.0);
	else
		n=n*factorial(n-1.0); // notice we just keep calling this over until we get to 1!
	return(n);
}

double SineN(const double& x, const int& n )// replaces first
{
	double sum = 0.0;

	if(n > 0)	
		for(int j=0; j<n; j++)
		{	// successive terms are to alternate sign
			if(j%2)// j is odd
				sum -= pow(x,(double)(2*j+1))/factorial((double)(2*j+1));
			else// j is even
				sum += pow(x,(double)(2*j+1))/factorial((double)(2*j+1));
		}

	return sum;
}// end of SineN()

// this method uses the recursion relation for the series
// to generate the next coefficient from the last.
// Advantages: 1) No need to find n! for each term.
// 2) No need to call pow(x,n) each term.
double SineN2(const double& x, const int& n, double* p_err )
{
	double a = 1.0*x;// 1st term
	double sum = 0.0;

	if(n > 0)
	{
		for(int j=1; j<=n; j++)// for sine
		{
			sum += a;
			a *= -1.0*x*x/( (2*j+1)*(2*j) );// next term from last
		//	sum += a;
		}
		*p_err = a;// max. error = 1st term not taken for alternating series
	}
	else
		cout << " n must be > 0";

	return sum;
}// end of SineN2()

double CosineN(const double& x, const int& n )// replaces second
{
	double sum = 1.0;

	if(n > 0)	
		for(int j=1; j<n; j++)
		{	
			if(j%2)// j is odd
				sum -= pow(x,(double)(2*j))/factorial((double)(2*j));
			else// j is even
				sum += pow(x,(double)(2*j))/factorial((double)(2*j));
		}

	return sum;
}// end of CosineN()

double CosineN2(const double& x, const int& n, double* p_err )
{
	double a = 1.0;// 1st term
	double sum = 0.0;

	if(n > 0)
	{
		for(int j=0; j<n; j++)// for cosine
		{
			sum += a;
			a *= -1.0*x*x/( (2*j+1)*(2*j+2) );// next term from last
		//	sum += a;
		}
		*p_err = a;// max. error = 1st term not taken for alternating series
	}
	else
		cout << " n must be > 0";

	return sum;
}// end of CosineN2()

int main()
{	
	int n = 0;// number of terms in series	
	double x = 0.0;// series parameter. Give value in radians
	double errorAmount = 0.0;// new functions find this value
	char repeat = 'y';
		
	do// until user is done
	{		
		cout << "Enter x (in radians): " ;
		cin >> x;
		cout << "Enter n: " ;// the # of terms in the series. Try entering 5 - 15 or so
		cin >> n;
		
		cout << "n! = " << factorial((double)n) << endl;
		cout << "DBL_MAX = " << DBL_MAX << endl << endl;
		cout << setprecision (10);// display result to 10 decimal places
		cout << "SineN  = "  << SineN(x,n) << endl;// using a power series
		cout << "SineN2 = "  << SineN2(x,n,&errorAmount) << endl;// using recursion relation		
		cout << "sin(x) = " << sin(x) << endl;// check against math library function
		cout << "errorAmount = " << errorAmount << endl << endl;// display max. error
		cout << "CosineN = " << CosineN(x,n) << endl;
		cout << "CosineN2 = " << CosineN2(x,n,&errorAmount) << endl;
		cout << "cos(x) = " << cos(x) << endl;
		cout << "errorAmount = " << errorAmount << endl << endl;

		cout << "repeat (y/n)? " << flush;
		cin >> repeat;
	}while( repeat=='y');
	return 0;
}

Have fun!!

[ne555]

If you use Taylor Series, fit the value to be calculated in the range [-pi, pi) x=fmod(x, PI);

[fruitloop2222]

Okay, thank you. How does the pow function work?

[ne555]

If the exponent is a float
z = x ^ y
ln z = ln (x^y)
ln z = y * ln x
z = exp (y*ln(x))

If the exponent is an integer use divide and conquer a^n = a^(n/2) * a^(n/2) * a^(n%2)

[jsmith]

If you are looking to learn the language, I would say no, because to write a math library requires a lot
of research into efficient algorithms to compute even simple things.

[kfmfe04]

+1 jsmith

If you are writing a math library as a throw-away, it's ok to try a small piece, but for production, what jsmith says is 100% correct - don't bother using your own.

For example, creating a Decimal class would be a nice programming exercise for learning C++, but for real production code, depend on researchers who spend a lot of time researching and thinking about these things:

http://speleotrove.com/decimal/

Another example of trying for fun, but using someone else's for real is linear algebra matrix operations.

I don't even know what to say for an entire math library, whatever that really means in terms of scope.

Could some one teach me how a robust power/exponent function works. I thought it would use 1 % exponent, and apparently you can't do modulo with decimal numbers (I tried). So then I decided that I could use the hack for modulo, which requires a floor function and blah blah blah... Math functions are quite complex!

BTW, learn the math before you try to use code. Use someone else's code (so you can test against) before you try to write your own. Make sure you understand the math first, because no matter how hard the math is, it's always easier to go from the math to the code - trying to deduce the math from the code still requires a deep understanding of the math; you also have to assume that the code is bug-free which is especially dangerous if you don't understand the math.

What the heck do you mean with modulo with decimal numbers and how is that related to power/exponent? Write down what you want to do in math terms and what result you expect to get back (an example or a test case). Then, we can try to take it from there.

[D80DSL3A]

The pow() is the only function used in the example which I didn't try to write my own function for.
It is a "black box" method which is supplied for us to use. I don't know how the math library version of it works.
I think jsmith has a valid point in that tackling a project like this would distract from learning the language - which is more important to do.
Also, any replacement functions I might come up with are not going to be nearly as good as the library versions.
I enjoy projects like this though - hence my goofy username!
Are you asking how we could write our own version of pow() ?

[fruitloop2222]

Actually, I've learned most of C++, me and my C++ reference book will do sufficiently. So, I was just doing this to get better at math... ;)

[D80DSL3A]

kfmfe04 said:

If you are writing a math library as a throw-away, it's ok to try a small piece...

Of course, it's just for playing around + satisfaction of seeing a series solution work.

+1 ne555

z = exp (y*ln(x)) is exactly the relation to use for the pow() here!
So I had to write functions for exp(x) and ln(x) first:

exp(x) is easy since the interval of convergence for the series is for x^2 < infinity:

// Finds e^x. x>0 or x<0 OK
double myExp(const double& x, const int& n )// e^x - no error finding
{	
	double sum = 0.0;// retVal
	double a = 1.0;// 1st term

	if(n > 0)	
		for(int j=1; j<=n; j++)
		{
			sum += a;
			a *= x/(double)j;	
		}		

	return sum;
}// end of myExp() 

The natural log function was harder due to the need to keep the series parameter (u in the function) within the interval of convergence = -1 to +1. Also, the convergence is slow so the # of terms (n) must be 50+ to get a result accurate to 10 figures:

// Give x > 0.0 only
double myNatLog(const double& x, const int& n )// ln(x) - no error finding
{	
	double   u = 1.0;// to substitute for x
	double   a = 0.0;// 
	double sum = 0.0;//

	if( (n > 0) && (x > 0.0) )
	{
		if( x < 2.0 )
		{
			u = 1 - x;// u will remain within interval of convergence ( u^2 < 1 )
			sum = a = u;// ist term
			for(int j=2; j<=n; j++)//
			{
				a *= u;
				sum += a/(double)j;			
			}
			return ( -1.0*sum );
		}
		else// x > 2.0
		{
			u = (x-1)/x;// u will remain within interval of convergence ( u^2 < 1 )
			sum = a = u;// ist term
			for(int j=2; j<=n; j++)
			{
				a *= u;
				sum += a/(double)j;			
			}
			return ( sum );
		}
	}

	return 0.0;// How to handle bad input?
}// end of myNatLog() 

Then the pow() can be used:

// Finds base^x. Must input base >0 . This works best for "smallish" base and exponent (<10 or so).
double myPow(const double& base, const double& x, const int& n )// base^x - no error finding
{
	double u = 1.0;
	if( (base>0.0) && (n>0) )
	{
		u = x*myNatLog(base,n);
		return ( myExp(u,n) );
	}
	else
		return 0.0;// error
}// end of myPow() 

These functions are just for playing around with power series. Do not use them as substitutes for the much more efficient (and safer) library functions. They are also crude 1st versions. The next challenge would be to improve them.
OK - I'll quit now and let the thread die.

[fruitloop2222]

Sorry I haven't checked for quite a while, but I have to say... WOWWWW!