Javascript squareroot algorithm and machine epsilon

[Jamin2112]

So I made a little application to show the steps of approximating a squareroot using the Babylonian Method:

http://jaminweb.com/squarerootCalc.html

It's not working all the time.

When I do the squareroot of 25 starting with an initial guess of 3.1, it works:


Applying x_n+1=(x_n²-S)/(2x_n) with x₀=3.1 and S=25 ...

x₁=5.582258064516129

x₂=5.030366246935904

x₃=5.000091654256142

x₄=5.000000000840035

x₅=5

Awesome.

Then the page crashes when I try 128371923 and 19. So there's something wrong.

Here's the loop I use to generate the output:

		for (var i = 1, xi = x0; Math.abs(Math.pow(xi,2) - S) >= eps; ++i)
		{
			xi = xi - (Math.pow(xi,2) - S)/(2*xi);
			thisStr += "<p>x<sub>" + i + "</sub>=<b>" + xi + "</b></p>";
		}

In the above code, eps is gotten by the function

function getMacheps()
{
	var temp1 = 1.0, temp2, mchEps;
	do {
		mchEps = temp1
		temp1 /= 2
		temp2 = 1.0 + temp1
	} while (temp2 > 1.0)
	return mchEps;
}

Are there any glaring fallacies?

 

[phinds]

128371923 is probably too big for the number system

 

[micromass]

128371923 is probably too big for the number system

19 isn't, so I doubt that is the reason.

Anyway, to solve your question. Could you post your Machine-Epsilon? And could you also post every xi that you get when you try this algorithm on a case that doesn't work (try 19).

 

[SteamKing]

I tried the routine to find the square root of 100,000 (10^5) using various starting points, but the routine failed.
If I use 99,999 (10^5 - 1), the routine worked with various starting values. This is where your problem starts.

 

[Mark44]

Jamin,
In your for loop, you use Math.pow() twice to calculate the square of a number. It's much more efficient to calculate x * x than it is to calculate Math.pow(x, 2). Raising a number to a power takes a lot more machine cycles than multiplication.

 

[Mark44]

I've never seen this called "Babylonian method" but have heard it called Newton's method. In a numerical methods class I took back in the last century, we learned that Newton's method can fail to converge to a root if you start off too far from the root or at some number that causes the approximations to go into an infinite loop. Either situation is not good. Certainly starting at 128371923 is way far away from the root (5).

 

[SteamKing]

Jamin,
In your for loop, you use Math.pow() twice to calculate the square of a number. It's much more efficient to calculate x * x than it is to calculate Math.pow(x, 2). Raising a number to a power takes a lot more machine cycles than multiplication.

Why are you using x(n+1) = x(n) - (x(n)^2 - S) / (2*x(n))

when the equivalent expression

x(n+1) = 0.5* (x(n) + S/x(n))

is simpler to evaluate? With the latter, you don't have to square any possible large numbers.

 

[Borg]

The short answer to your problem is that the getMacheps() is casting the variables as doubles which results in an extremely low number that won't allow the for loop to exit in random cases. If S has a square root that is an exact int, it always works. However, it fails for other cases due to rounding issues.

Change your getMacheps() to this:

function getMacheps()
{
	var temp1 = 1, temp2, mchEps;
	do {
		mchEps = temp1
		temp1 /= 2
		temp2 = 1 + temp1
	} while (temp2 > 1)
	return mchEps;
}

Now the long answer:

I don't understand why you have a method like getMacheps() that returns the same value every time. Do you know what it returns? I rewrote the method in a Java class and because of the 1.0 values in getMacheps(), the return value has to be cast as a double. The return value is eventually mchEps: 2.220446049250313E-16 (after 53 times through the loop). However, if I change the 1.0 values to 1, the javascript will cast them as ints which results in a return value of mchEps: 1. Which value are you expecting?

Also, just because you can initialize multiple variable declarations inside a for loop's declaration section doesn't mean that you should. In the real world, others may have to debug your code and trying to debug that for loop is a pain. Your for loop is just a overly complex while loop that is initializing all of your variables. And while it looks like the variable 'i' has something to do with the equation, it doesn't. It's just a counter for your output. Your for loop is really nothing more than this (written in Java for variable clarity):

private void approximationLoop(double xi, double S) {
	// eps is always 1 or 2.220446049250313E-16 based on getMacheps() using ints or doubles
	double eps = getMacheps();
	int i = 1;
	String thisStr = "";

	while(Math.abs(Math.pow(xi, 2) - S) >= eps){
		xi = xi - (Math.pow(xi, 2) - S) / (2 * xi);
		thisStr += "<p>x<sub>" + i++ + "</sub>=<b>" + xi + "</b></p>";
	}
}

Starting with x0 = 19 and S = 128371923.

If I force getMacheps() to use and return ints, I get the following result after 13 iterations:

S: 128371923
eps: 1.0
xi: 11330.133406125462
Math.pow(xi, 2) result: 1.2837192300060016E8
Math.pow(xi, 2) - S result: 6.001591682434082E-4
Status of while term (Math.abs result >= eps): false
<p>x₁₃=<b>11330.133406125462</b></p>

If I let getMacheps() use doubles, the for loop never ends because the Math.abs result stabilizes around 1.4901161193847656E-8 which is greater than the eps value of 2.220446049250313E-16.

After 22278 passes through the loop...

S: 128371923
eps: 2.220446049250313E-16
xi: 11330.133406098978
Math.pow(xi, 2) result: 1.2837192300000001E8
Math.pow(xi, 2) - S result: 1.4901161193847656E-8
Status of while term (Math.abs result >= eps): true
<p>x₂₂₂₇₈=<b>11330.133406098978</b></p>

Using SteamKing's example of 99,999, even the double version of getMacheps() works because the Math.pow(xi, 2) result equals S exactly which allows the loop to exit:

S: 99999
eps: 2.220446049250313E-16
xi: 316.226184874055
Math.pow(xi, 2) result: 99999.0
Math.pow(xi, 2) - S result: 0.0
Status of while term (Math.abs result >= eps): false
<p>x₉=<b>316.226184874055</b></p>

But, a much smaller pair of starting values fails (x0 = 2 and S = 5):

S: 5
eps: 2.220446049250313E-16
xi: 2.23606797749979
Math.pow(xi, 2) result: 5.000000000000001
Math.pow(xi, 2) - S result: 8.881784197001252E-16
Status of while term (Math.abs result >= eps): true
<p>x₃₄₄₀₄=<b>2.23606797749979</b></p>

 

[Jamin2112]

Jamin,
In your for loop, you use Math.pow() twice to calculate the square of a number. It's much more efficient to calculate x * x than it is to calculate Math.pow(x, 2). Raising a number to a power takes a lot more machine cycles than multiplication.

I was going by the "Optimize readability" and "There's no need to optimize JavaScript code" advice from Stack Overflow.

 

[Mark44]

I was going by the "Optimize readability" and "There's no need to optimize JavaScript code" advice from Stack Overflow.

I would maintain that x * x is eminently more readable than Math.pow(x, 2), as well as being quicker to type, not to mention way faster to execute. As far as there being no need to optimize JavaScript code, that's an opinion, not something that is chiseled into stone. I'm not talking about optimizing the entire code, but using pow() to do a simple multiplication is definitely not optimal, especially inside a loop that can run many times.

 

[Jamin2112]

I've never seen this called "Babylonian method".

I got that from Wikipedia and thought it sounded cooler than "Newton's Method."

http://en.wikipedia.org/wiki/Methods_of_computing_square_roots

In a numerical methods class I took back in the last century, we learned that Newton's method can fail to converge to a root if you start off too far from the root or at some number that causes the approximations to go into an infinite loop. Either situation is not good. Certainly starting at 128371923 is way far away from the root (5).

Dang it. I might have to add some code so that the program "times out" if it hasn't found a solution.

 

[Jamin2112]

The short answer to your problem is that the getMacheps() is casting the variables as doubles which results in an extremely low number that won't allow the for loop to exit in random cases. If S has a square root that is an exact int, it always works. However, it fails for other cases due to rounding issues.

Change your getMacheps() to this:

function getMacheps()
{
	var temp1 = 1, temp2, mchEps;
	do {
		mchEps = temp1
		temp1 /= 2
		temp2 = 1 + temp1
	} while (temp2 > 1)
	return mchEps;
}

Now the long answer:

I don't understand why you have a method like getMacheps() that returns the same value every time. Do you know what it returns? I rewrote the method in a Java class and because of the 1.0 values in getMacheps(), the return value has to be cast as a double. The return value is eventually mchEps: 2.220446049250313E-16 (after 53 times through the loop). However, if I change the 1.0 values to 1, the javascript will cast them as ints which results in a return value of mchEps: 1. Which value are you expecting?

Also, just because you can initialize multiple variable declarations inside a for loop's declaration section doesn't mean that you should. In the real world, others may have to debug your code and trying to debug that for loop is a pain. Your for loop is just a overly complex while loop that is initializing all of your variables. And while it looks like the variable 'i' has something to do with the equation, it doesn't. It's just a counter for your output. Your for loop is really nothing more than this (written in Java for variable clarity):

private void approximationLoop(double xi, double S) {
	// eps is always 1 or 2.220446049250313E-16 based on getMacheps() using ints or doubles
	double eps = getMacheps();
	int i = 1;
	String thisStr = "";

	while(Math.abs(Math.pow(xi, 2) - S) >= eps){
		xi = xi - (Math.pow(xi, 2) - S) / (2 * xi);
		thisStr += "<p>x<sub>" + i++ + "</sub>=<b>" + xi + "</b></p>";
	}
}

Starting with x0 = 19 and S = 128371923.

If I force getMacheps() to use and return ints, I get the following result after 13 iterations:

S: 128371923
eps: 1.0
xi: 11330.133406125462
Math.pow(xi, 2) result: 1.2837192300060016E8
Math.pow(xi, 2) - S result: 6.001591682434082E-4
Status of while term (Math.abs result >= eps): false
<p>x₁₃=<b>11330.133406125462</b></p>

If I let getMacheps() use doubles, the for loop never ends because the Math.abs result stabilizes around 1.4901161193847656E-8 which is greater than the eps value of 2.220446049250313E-16.

After 22278 passes through the loop...

S: 128371923
eps: 2.220446049250313E-16
xi: 11330.133406098978
Math.pow(xi, 2) result: 1.2837192300000001E8
Math.pow(xi, 2) - S result: 1.4901161193847656E-8
Status of while term (Math.abs result >= eps): true
<p>x₂₂₂₇₈=<b>11330.133406098978</b></p>

Using SteamKing's example of 99,999, even the double version of getMacheps() works because the Math.pow(xi, 2) result equals S exactly which allows the loop to exit:

S: 99999
eps: 2.220446049250313E-16
xi: 316.226184874055
Math.pow(xi, 2) result: 99999.0
Math.pow(xi, 2) - S result: 0.0
Status of while term (Math.abs result >= eps): false
<p>x₉=<b>316.226184874055</b></p>

But, a much smaller pair of starting values fails (x0 = 2 and S = 5):

S: 5
eps: 2.220446049250313E-16
xi: 2.23606797749979
Math.pow(xi, 2) result: 5.000000000000001
Math.pow(xi, 2) - S result: 8.881784197001252E-16
Status of while term (Math.abs result >= eps): true
<p>x₃₄₄₀₄=<b>2.23606797749979</b></p>

Thanks for the informative answer. I'll make the appropriate changes later today.

My buddy, a very experienced programmer, uses multiple declarations in for-loops. I always thought it was "fancy" to fully exploit for loops, have lines like

node* thisPtr = root;
for (; thisPtr; thisPtr = thisPtr->next);

in the middle of one's code.

 

[Borg]

Thanks for the informative answer. I'll make the appropriate changes later today.

My buddy, a very experienced programmer, uses multiple declarations in for-loops. I always thought it was "fancy" to fully exploit for loops, have lines like

node* thisPtr = root;
for (; thisPtr; thisPtr = thisPtr->next);

in the middle of one's code.

So what value did you expect to get from the getMacheps() function?

 

[Mark44]

Thanks for the informative answer. I'll make the appropriate changes later today.

My buddy, a very experienced programmer, uses multiple declarations in for-loops. I always thought it was "fancy" to fully exploit for loops, have lines like

node* thisPtr = root;
for (; thisPtr; thisPtr = thisPtr->next);

in the middle of one's code.

What multiple declarations are you talking about? The code you show isn't a declaration. What it does is advance thisPtr to the end of the linked list.

Inside the parentheses, the initialization expression is empty. The test expression -- thisPtr -- is considered true as long as the value of thisPtr is not null. The increment section -- thisPtr = thisPtr -> next -- advances the pointer to the next node in the list.

The body of the for loop is empty, as indicated by the semicolon immediately following the closing parenthesis in the for loop header. Having the semicolon on the same line as the for loop header is considered by some to be a bit dangerous, as it's easy to overlook it and think that a following line is part of the loop body.

One way to make things more obvious is to put the semicolon on a different line.

for (; thisPtr; thisPtr = thisPtr->next)
   ;

I'm reasonably sure that many software firms have coding standards that mandate using braces for all for loops, even when the loop bodies are empty, like so:

for (; thisPtr; thisPtr = thisPtr->next)
{
   // empty
}

 

[D H]

Are there any glaring fallacies?

Two, both here:

for (var i = 1, xi = x0; Math.abs(Math.pow(xi,2) - S) >= eps; ++i)
{
   xi = xi - (Math.pow(xi,2) - S)/(2*xi);
   thisStr += "<p>x<sub>" + i + "</sub>=<b>" + xi + "</b></p>";
}

What's wrong? The biggest problem is your termination condition. You might never exit your loop.

Do you understand how floating point numbers work? Machine epsilon is the difference between the smallest representable number that is greater than one and one itself. The smallest representable number that is greater than 2¹⁶ and 2¹⁶ itself? That's not machine epsilon. It's one. You have to either somehow scale your tolerance or scale the comparison. The former would involve comparing |x^2-S| to S*epsilon, the latter would involve comparing |x^2/S-1.0| to epsilon (but now you have to worry about S=0).

There's still a problem here: you don't want to compare |x^2-S| to S*epsilon. You still might never exit the loop. Suppose your $S$ is perfectly ill-suited, that the two representable numbers closest to $\sqrt{S}$ are $\sqrt{S}(1-\epsilon/2)$ and $\sqrt{S}(1+\epsilon/2)$. Your algorithm will eventually reach one of these two values and then will alternate between them. In both cases, $|x^2/S-1.0|$ will be machine epsilon. Compare $|x^2/S-1.0|$ to machine epsilon and you never exit the loop.

What needs to be done is to exit the loop with a tolerance larger than machine epsilon, say 100*S*epsilon, and then perform one last iteration outside of the loop. Suppose you exit the loop with $x_n=\sqrt{S}(1+\Delta)$, where $\Delta$ is some small number. The next refinement will yield $x_{n+1}=\sqrt{S}(1-\Delta^2/4-O(\Delta^3))$. So long as $\Delta$ is smaller than $\sqrt{\epsilon}$, the next step will bring you within machine epsilon.


The other problem is your iteration. You used $x_i = x_i - \frac{x_i^2-S} {2x_i}$. It's easy to see that this is equivalent to $x_i = \frac 1 2 \left(x_i + \frac S {x_i}\right)$. This is a much more stable computation; it is what you should be using.

 

[Borg]

I'm reasonably sure that many software firms have coding standards that mandate using braces for all for loops, even when the loop bodies are empty, like so:

for (; thisPtr; thisPtr = thisPtr->next)
{
   // empty
}

I wish that were true.

 

[AlephZero]

The other problem is your iteration. You used $x_i = x_i - \frac{x_i^2-S} {2x_i}$. It's easy to see that this is equivalent to $x_i = \frac 1 2 \left(x_i + \frac S {x_i}\right)$. This is a much more stable computation; it is what you should be using.

Not only is it more stable, it works for a bigger range of values.

64-bit floating point numbers can represent values between about $\pm 10^{\pm 308}$, and zero.

Your algorithm will fail for $x > 10^{154}$ (because $x^2$ will cause an overflow) or $x < 10^{-154}$ (because $x^2$ will underflow and return zero).

D.H's version will work for the full range of $x$ - or at least a lot closer to the full range.

You might not think this matters - but if you ever use 32-bit floating point where the exponent limit is only $10^{\pm 38}$, not being able to square root $10^{20}$ could be a bug, not just a feature.

This isn't just a quibble. In the past computer companies which should to have known better implemented arithmetic for complex numbers by calculating nonsense like $$\frac{1}{x + iy} = \frac{x - iy}{x^2 + y^2}.$$ Spot the errors (plural!)

 

[SteamKing]

The other problem is your iteration. You used $x_i = x_i - \frac{x_i^2-S} {2x_i}$. It's easy to see that this is equivalent to $x_i = \frac 1 2 \left(x_i + \frac S {x_i}\right)$. This is a much more stable computation; it is what you should be using.

Not only is it more stable, it works for a bigger range of values.

64-bit floating point numbers can represent values between about $\pm 10^{\pm 308}$, and zero.

Your algorithm will fail for $x > 10^{154}$ (because $x^2$ will cause an overflow) or $x < 10^{-154}$ (because $x^2$ will underflow and return zero).

D.H's version will work for the full range of $x$ - or at least a lot closer to the full range.

I tried to tell the OP about this back in Post #7. Might as well been whistlin in a graveyard.

 

[Jamin2112]

I tried to tell the OP about this back in Post #7. Might as well been whistlin in a graveyard.

I'm going to change it when I get on my desktop computer tomorrow.

 

[Jamin2112]

To whoever was asking about the value of macheps, its

2.220446049250313e-16

 

[Jamin2112]

New procedure:

	var S = parseFloat(document.getElementById("babyS").value);
	var x0 = parseFloat(document.getElementById("babyX0").value);
	var thisStr = new String();
	if (S == NaN || x0 == NaN)
	{
		alert("Invalid input.");
		return;
	}
	else if (x0 == 0)
	{
		alert("For obvious reasons, the algorithm does not work with an initial guess of 0!");
	}
	else
	{	
		thisStr = "<p>Applying x<sub>n+1</sub>=(x<sub>n</sub><sup>2</sup>-S)/(2x<sub>n</sub>) with x<sub>0</sub>=<b>" + x0.toString() + "</b> and S=<b>" + S.toString() + "</b> ...</p>";
		var exitval = 100 * S * macheps, i = 1, xi = x0;
		for (; Math.abs(xi * xi - S) >= exitval; ++i)
		{
			xi = (xi + S/xi)/2;
			thisStr += "<p>x<sub>" + i + "</sub>=<b>" + xi + "</b></p>";
		}
		xi = (xi + S/xi)/2;
		thisStr += "<p>x<sub>" + i + "</sub>=<b>" + xi + "</b></p>";
	}
	document.getElementById("babsteps").innerHTML = thisStr;

where

macheps

is a global variable created on page load.

(And I know someone's going to tell me I should be using JQuery for the element modifications. I'm trying to survive without familiarizing myself with the JQuery library)

 

[D H]

You might want to protect against negative inputs.

 

[Jamin2112]

You might want to protect against negative inputs.

I could just add an "i" at the end of each line of the output

 

[D H]

Try feeding your algorithm a negative value for S.

Try feeding S=-1 and x0=-1. You should get a divide by zero.
Try feeding S=-19 and x0=-1. Your string will overflow because of failure to converge.

Another nasty test case is S=1e-300 and x0=1e300. Your string may overflow because of too many iterations.
An even nastier one is S=1e300 and x0=1e-300. Your string will overflow because the iteration will get stuck at infinity.

You will find that it is a good idea to limit the number of iterations and to protect against NaNs and infinities.