My Taylor Square Root C Program doesn't like me

[Orcinus]

Homework Statement

4. Implement a simple method to find the square root of a double precision floating point number x. A simple method is to consider the error produced by a “guess” y of the solution. Square the value y and compare with the value x. If y is correct, the error e=|y2-x| where || is the magnitude will be zero and as the value of y becomes more wrong, the error e will increase. Hence the solution is to home in on the minimum for the error function or some value close enough for our purposes. This problem actually has two solutions +/-y but we can just consider +y. Write a C program to solve this problem using a repetition method. Hint: use a while loop and devise a way to find when you get to the minimum for the error e.

This is in C by the way (Linux gcc) and obviously without using maths libraries.

The question itself is simple, I quickly devised the solution that the question's looking for: http://pastebin.com/cQGzGXN4
But I was curious of doing it another way, and after just learning about Taylor series in maths I wanted to try doing it that way. So far it's not working out so well.

Course: Engineering Programming 100, we've yet to even learn functions yet so yeah I'm going beyond the course and hence beyond my knowledge level.

Homework Equations

I don't know whether this is right, but a plot on www.desmos.com/calculator seems to work
f(x) = \sqrt{x}
<br /> f(x) = 1+ \sum_{n=1}^\infty\frac{\frac{(2n-2)!(-1)^{n-1}(x-1)^n}{2^{2n-1}(n-1)!}}{n!}<br />

The Attempt at a Solution

Here's what I've come up with:

#include <stdio.h>

double power(double base, int index) //Note that this function does not allow fractional indicies
{
	int x;
	int inverse = 0;
	double result = 1;
	
	if(index<0)
	{
		index = -index;
		inverse = 1;
	}
	if(index==0)
	{
		return 1;
	}
	for(x=0;x<index;x++)
	{
		result = base*result;
	}
	if(inverse == 1)
	{
		return (1/result);
	}
	else
	{
		return result;
	}
}

int factorial(int num)
{
	int x;
	int result = 1;
	for(x=1;x<=num;x++)
	{
		result = result * x;
	}
	return result;
}

double absolute( double num)
{
	double result;
	result = num;
	if(result < 0)
	{
		result = -result;
	}
	return result;
}

double squareroot(double square)
{
	int x = 1;
	double result = 1;
	while(absolute(result*result-square) > 0.000001)
	{ //hello taylor series! I have no idea whether this works.
		result = result + factorial(2*x-2)*power(-1, x-1)/power(2,2*x-1)/factorial(x-1)*power(square - 1, x)/factorial(x);
		x++;
		printf("%lf\n", result);
	}
	return result;
}

int main(void)
{
	double square;
	double result;
	printf("Type the number you'd like to find the square root of: ");
	//scanf("%lf", square);
square = 2; // for debugging
	result = squareroot(square);
	printf("The square root of %lf is %lf.\n", square, result);
	return 0;
}

Note that "printf("%lf\n", result);" in squareroot function and "square = 2; // for debugging" in main function are for debugging purposes. I'm writing this on my tablet so can't compile properly but am using http://codepad.org/W5yNE67K to debug. Hence I manually enter numbers into "square = 2" to test different numbers.

The problem being that my result is coming out nothing like \sqrt{x}, the while loop spits out a NaN for some reason after a while and this brings together so much new stuff for me that I don't know where the problem(s) might lie.

 

[DrClaude]

I think your problem is with the Taylor series. If you check WolframAlpha, you will find that you only get convergence for $|1-x|<1$.

 

[Orcinus]

Aaw really? Is that just a thing with square root or is it just my method isn't that great?

 

[DrClaude]

I'm not a mathematician, but from what I gather it won't work for the square root. Other functions could still be programmed that way.

If you're more interested in general about these kind of numerical methods, I can point two possibile difficulties with series expansion. First, if you implement it naively from $n=0$ to whatever, you are progressively adding smaller and smaller terms, and numerical errors can creep up because of that. This is especially important here, and it touches on my second point, where subsequent terms have different signs, and therefore can cancel each other out. This can also be a problem even when the terms in the series do not progressively decrease, but there is a maximum somewhere in the sum. You can again lose precision if some big term kill off some significant digits before this big term gets canceled out by another big term.

I don't know if you can make sense of what I just wrote. If not, just say so and I will try to come up with specific examples.

 

[DrClaude]

I should also say that I just had a quick look at your program, and their might be some problems in there too.

 

[cepheid]

Yeah for one thing, in the power function, saying return(1/result) will probably give you 0, even though the return type is double. I suspect that it will do the integer arithmetic first, and then cast the result of that to a double, which is not what you want.

Things like power and factorial really lend themselves to *recursion*. You should look it up. It's a technique where a function calls itself.

I could be wrong, but I think that statements like index = -index are inefficient. In any case, you can just write index *= -1.

Regarding the Taylor series: regardless of the value of x, you should be able to get your answer arbitrarily close to sqrt(x) as long as you take enough terms in the sum.

EDIT: oh I didn't think about convergence though. Would have to refresh my memory about that.

 

[Orcinus]

So the problem was a conceptual one then? That's a shame.
Yeah, I figured that imprecision with floats would come into effect with a series, but I was hoping that (-1)^n might bounce the number around a bit til it was close enough to the right value. But statistics doesn't really let that happen too easily.
So would: "result = 1/result; return result;" make any difference?
I thought x *= y was just the same as x=x*y, shouldn't the compiler recognise that and choose the quickest method?

Well the program conceptually doesn't work, so there's not much point "fixing it". But now I've had a go, do you know where I could find out how <math.h> does it?

 

[Mark44]

Yeah for one thing, in the power function, saying return(1/result) will probably give you 0, even though the return type is double. I suspect that it will do the integer arithmetic first, and then cast the result of that to a double, which is not what you want.

Not in this case. result is declared as a double, so in the division, 1 is promoted to a double value and then the division is performed. If result had been declared as an integral type (e.g., int, long, char, etc.), then 1/result would be calculated using integer division.

Things like power and factorial really lend themselves to *recursion*. You should look it up. It's a technique where a function calls itself.

I could be wrong, but I think that statements like index = -index are inefficient. In any case, you can just write index *= -1.

Speaking of inefficiency, recursive routines can be much more inefficient than loops, as there is the extra overhead of multiple calls to the recursive function. As far as index = -index vs. index *= -1, I doubt that there is much difference. Also, since this statement is executed only once (as opposed to being executed thousands or more times in a loop), the difference in execution speed would be insignificant.

Regarding the Taylor series: regardless of the value of x, you should be able to get your answer arbitrarily close to sqrt(x) as long as you take enough terms in the sum.

EDIT: oh I didn't think about convergence though. Would have to refresh my memory about that.