Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Javascript squareroot algorithm and machine epsilon

  1. Apr 4, 2014 #1
    So I made a little application to show the steps of approximating a squareroot using the Babylonian Method:

    http://jaminweb.com/squarerootCalc.html [Broken]

    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 xn+1=(xn2-S)/(2xn) with x0=3.1 and S=25 ...

    x1=5.582258064516129

    x2=5.030366246935904

    x3=5.000091654256142

    x4=5.000000000840035

    x5=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:

    Code (Text):
            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

    Code (Text):
    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?
     
    Last edited by a moderator: May 6, 2017
  2. jcsd
  3. Apr 4, 2014 #2

    phinds

    User Avatar
    Gold Member
    2016 Award

    128371923 is probably too big for the number system
     
  4. Apr 4, 2014 #3

    micromass

    User Avatar
    Staff Emeritus
    Science Advisor
    Education Advisor
    2016 Award

    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).
     
  5. Apr 4, 2014 #4

    SteamKing

    User Avatar
    Staff Emeritus
    Science Advisor
    Homework Helper

    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.
     
  6. Apr 5, 2014 #5

    Mark44

    Staff: Mentor

    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.
     
  7. Apr 5, 2014 #6

    Mark44

    Staff: Mentor

    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).
     
  8. Apr 5, 2014 #7

    SteamKing

    User Avatar
    Staff Emeritus
    Science Advisor
    Homework Helper

    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.
     
    Last edited: Apr 5, 2014
  9. Apr 5, 2014 #8

    Borg

    User Avatar
    Gold Member

    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:
    Code (Text):

    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):

    Code (Text):

    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<sub>13</sub>=<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<sub>22278</sub>=<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<sub>9</sub>=<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<sub>34404</sub>=<b>2.23606797749979</b></p>
     
    Last edited: Apr 5, 2014
  10. Apr 5, 2014 #9
    I was going by the "Optimize readability" and "There's no need to optimize JavaScript code" advice from Stack Overflow.
     
  11. Apr 5, 2014 #10

    Mark44

    Staff: Mentor

    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.
     
  12. Apr 5, 2014 #11
    I got that from Wikipedia and thought it sounded cooler than "Newton's Method."

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



    Dang it. I might have to add some code so that the program "times out" if it hasn't found a solution.
     
  13. Apr 5, 2014 #12

    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.
     
  14. Apr 5, 2014 #13

    Borg

    User Avatar
    Gold Member

    So what value did you expect to get from the getMacheps() function?
     
  15. Apr 5, 2014 #14

    Mark44

    Staff: Mentor

    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.
    Code (Text):
    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:
    Code (Text):
    for (; thisPtr; thisPtr = thisPtr->next)
    {
       // empty
    }
     
  16. Apr 5, 2014 #15

    D H

    User Avatar
    Staff Emeritus
    Science Advisor

    Two, both here:

    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 216 and 216 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.
     
    Last edited: Apr 5, 2014
  17. Apr 5, 2014 #16

    Borg

    User Avatar
    Gold Member

    I wish that were true.
     
  18. Apr 5, 2014 #17

    AlephZero

    User Avatar
    Science Advisor
    Homework Helper

    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!)
     
    Last edited: Apr 5, 2014
  19. Apr 5, 2014 #18

    SteamKing

    User Avatar
    Staff Emeritus
    Science Advisor
    Homework Helper


    I tried to tell the OP about this back in Post #7. Might as well been whistlin in a graveyard.
     
  20. Apr 5, 2014 #19
    I'm gonna change it when I get on my desktop computer tomorrow.
     
  21. Apr 7, 2014 #20
    To whoever was asking about the value of macheps, its

    2.220446049250313e-16
     
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook




Similar Discussions: Javascript squareroot algorithm and machine epsilon
  1. Javascript problem (Replies: 6)

Loading...