Fast Square Root Calculation?

[Alphanos]

Can anyone make this faster?

Public Function SquareRoot(ByVal Number As Double) As Double
 
If Number < 0 Then Number = Number * -1
 
Dim i As Long
Dim Temp As Double
Dim Power As Long
Dim Estimate As Double
 
Temp = 1
'Calculate an initial Estimate
If Number >= 1 Then 
    Do While Number > Temp
        Temp = Temp * 10
        Power = Power + 1
    Loop
Else
    Do While Number < Temp
        Temp = Temp / 10
        Power = Power - 1
    Loop
End If
 
Power = Power / -2
Estimate = 10 ^ Power * Number
 
'Calculate actual square root
Dim OldEstimate As Double
 
Do
    OldEstimate = Estimate
    Estimate = 0.5 * (Estimate + Number/Estimate)
Loop Until OldEstimate = Estimate
 
SquareRoot = Estimate
End Function

(Wrote this just now, so I haven't debugged it; I'll go check it now but I want to post it first just in case I get a BSOD.)

[Alphanos]

Seems to work okay (produces same results as windows calculator for numbers I've checked, whole, fractional, and tiny), but I realized that I hadn't actually used i, LOL. I'm so used to using i as a counter that I declared it for the loop, even though it wasn't needed.

[Yonatan]

There's already a VB function that approximates the square root, I think it's called Sqr. Unless you really have to use your own function, I suggest you use it, since it takes advantage of the FPU.

[Alphanos]

If VB didn't have a square root function, you could still use ^0.5.

This is actually a precursor to developing (relatively) fast functions that will do basic math operations on arbitrarily accurate numbers. I intend (at some point in the next few months, lol, lack of time) to create a class capable of storing a number to any desired accuracy, for use in testing my method of calculating pi. I need to be able to do fast, accurate square roots for some trig identities.

PS: Thanks to you guys who responded to that sine calculation thread, you helped me to look at the problem in a different way such that the same problem could be solved a bit differently.

[Yonatan]

Originally posted by Alphanos
I need to be able to do fast, accurate square roots for some trig identities.

Oh... Well, you should definitely not be using VB, then.

[A$$Bandit]

Public Function SquareRoot(ByVal Number As Double) As Double
    Dim i As Integer
    Dim X As Double
    X = Number / 2
    For i = 1 To 12
        X = (X + (Number / X)) / 2
    Next i
    SquareRoot = X
End Function

Now that HAS to be quicker. Innaccuracy > 1 after about Number = 1 million, but of course you can just change the number of iterations of the For loop to suit the number you're squarerooting. More iterations, more accurate etc you get the idea. Now I'm fairly sure that's the quickest way of doing it.

[Alphanos]

I think the difference in speed would depend on what type of numbers you're using. The main part of my function is essentially the same, except that it compares the new value of x to the previous and only stops when they are equal (in other words, when it can't get any more accurate).

The long part of the code at the beginning is to calculate a close initial estimate, something that is wasteful with normal numbers, but time-saving for numbers such as 5.67 * 10^153. I had another idea to maybe improve speed in most cases that's different from this one, I'll check it out and post later.

[jim mcnamara]

This isn't VB, but it is what VB calls --

Here is a page which has a lot of the best known algorithms to find square roots:

http://www.azillionmonkeys.com/qed/sqroot.html

[Guv]

It might be difficult to beat the Sqr function. The code is probably written in assembler, allowing for trickery. Division by two can be done by subtracting one from the floating point exponent instead of actually dividing. A good first approximation can be obtained by using the number itself with the half the floating point exponent.

Methods like x^.5 or logs should not be considered. I never heard of exponential or log functions being available when Square Root was not. A Square Root function will be faster.

The Newton iteration is so fast that you are likely to waste time trying to get a very good first approximation.

A first approximation like (1 + Number)/4 is a good compromise. It is not terrible for numbers greater than one or less than one.

Requiring two successive approximations to be equal might cause an infinite loop for some values due to unlucky round off error. I would not use that test without some numerical analysis to make certain that an infinite loop would not occur.

If you are looking for speed as well as accuracy, remember that the number of correct digits doubles each time you iterate. When the first 5 digits of two successive approximations agree, the next approximation will be accurate to ten digits. It is probably worthwhile to check for an good approximation and then do one more iteration without a test for convergence.

[Yonatan]

Originally posted by Guv
It might be difficult to beat the Sqr function. The code is probably written in assembler, allowing for trickery.

The Sqr function, as I already said, uses the FPU. It's just a wrapper around the FPU instruction fsqrt.
Here is the disassembly, commented by me:

; Taken from MSVBVM60.DLL

; Data at virtual address 66030248

Zero dq 0.0

; Export rtcSqr, virtual address 660E1B2E

                fld     qword ptr [esp+4] ; Load float parameter to FPU stack
                fcomp   [Zero] ; Compare parameter to zero
                fnstsw  ax ; Store FPU status word in AX without checking for FPU exceptions
                sahf       ; Store some bits of AH in some of EFLAGS
                jnb     Continue ; Jump if CF=0
                           ; Note - the above jump occurs if the parameter is bigger than or equal to 0.0
                push    5 ; Parameter for ErrRaise
                call    ErrRaise ; Err.Raise 5 (Invalid procedure call or argument, I think)
Continue:
                fld     qword ptr [esp+4] ; Load float parameter to FPU stack (it was unloaded by fcomp)
                fsqrt      ; Calculate the square root
                retn    8  ; Return

As you can see, all VB really does is the input validity checking. The FPU does the work of the square root. It is guaranteed that the FPU does it fastest.

[Alphanos]

A good first approximation can be obtained by using the number itself with the half the floating point exponent.

That's basically what I was doing, although my way is somewhat inefficient.

BTW, to reiterate, I'm not trying to make a function to beat VB's Sqr using standard floats or doubles, etc. I'm trying to determine fairly fast methods which I will later alter for use in special-purpose classes designed for storing very large/small numbers without floating points, to very high accuracy.

That's the purpose of requiring two sucessive iterations to produce the same result: ensuring that the operation only ends when the result is as accurate as possible using however many decimal places it is using. Since it won't be using floating-point, I didn't think rounding error would be a problem?

The other option I was considering instead of my slow first approximation thingy is to start with an estimate of one, run it through the standard equation once, and then run it a few times through a modified version. ie.

x = 0.125 * ( X + 7*Y/X)

That would essentially halve it thrice for numbers when 1 is very far away from the actual value.

ie.
1 initially as a guess for the root of 65536 would become 32768.5 from the first normal equation, then skip down to about 4097.81, then 526.22. Maybe after that run it twice in quarters instead of eigths, changing the multiplier to 3.

What do you think of this idea? It certainly adds steps for some numbers, but I thought it might help for most.

[PT Exorcist]

i am afraid vb6 doesnt support ^ with floating types..i tryed
dim i as double
i = 9 ^ 0,5

and it didnt work!

Originally posted by Alphanos
If VB didn't have a square root function, you could still use ^0.5.

This is actually a precursor to developing (relatively) fast functions that will do basic math operations on arbitrarily accurate numbers. I intend (at some point in the next few months, lol, lack of time) to create a class capable of storing a number to any desired accuracy, for use in testing my method of calculating pi. I need to be able to do fast, accurate square roots for some trig identities.

PS: Thanks to you guys who responded to that sine calculation thread, you helped me to look at the problem in a different way such that the same problem could be solved a bit differently.

[Yonatan]

Use . not ,

[Guv]

PT Exorcist: My Test Application indicated that VB ^ operator works fine for most any values. The following shows the code without supporting subroutines and the form.

Dim A As Double
Dim B As Double
Dim Result As Double
 
A = MakeDouble(txtA.Text)
txtA.Text = Dformat(A)
B = MakeDouble(txtB.Text)
txtB.Text = Dformat(B)
Result = A ^ B
txtC.Text = Dformat(Result)

Alphanos:

That's the purpose of requiring two sucessive iterations to produce the same result: ensuring that the operation only ends when the result is as accurate as possible using however many decimal places it is using. Since it won't be using floating-point, I didn't think rounding error would be a problem?

Note the following Newton method successive approximations for the square root of 4 using 4 as first estimate.

• 4.000 000 00000
• 2.500 000 00000
• 2.050 000 00000
• 2.000 609 75610
• 2.000 000 09292
• 2.000 000 00000

Note that you need not require two consecutive estimates to be equal . Once you get close, the precision doubles with each iteration.

Suppose you wanted 100 digits of precision and speed was critical. You could loop until you had two estimates which matched in the first 25-26 digits (a bit more than one quarter of 100 digits). Then calculate two more estimates without testing for convergence. The last two calculations would be faster due to no testing. If you are programming extra precision arithmetic, comparing the initial parts of consecutive estimates is not difficult.

Note that (n + 1)/2 is the second estimate, if you start with n as the first estimate. In the above example, the third estimate is 2.050 which is a fairly good estimate. A very fast way to get a good initial value is to do 2-3 Newton calculations without testing for convergence, starting with the number itself as the initial estimate.

I think that the further the number is from one, the worse it is for use as a first approximation. The closer the number is to one, the better it is. This suggests that three algorithms might be useful for getting an initial estimate.

• One method for very small numbers.
• The second for numbers in the neighborhood of one.
• The third for very large numbers.

[Alphanos]

A very fast way to get a good initial value is to do 2-3 Newton calculations without testing for convergence, starting with the number itself as the initial estimate.

I think that the further the number is from one, the worse it is for use as a first approximation.

Using 1 as a first estimate is equivalent to using the number itself.
ie. number = x

x as estimate:
y = 0.5 * (x + x/x)
= 0.5 * (x + 1)

1 as estimate:
y = 0.5 * (1 + x/1)
= 0.5 * (1 + x)

So with an initial estimate of either 1 or x, the second estimate will be 0.5 * (x+1) .

Thanks for the idea about checking only the first bunch of numbers for convergence. Do you know how far out this effect works with very large numbers of digits? IE. If I have 1000 digits, and the first 64 are matching currently, will they require 4 more estimates to match? Will it require another 24 to match 10^9 after 64 match already?

I'm reasonably sure that the idea of starting with either the number or one as an intial estimate, then doing one normal estimate, 2 at 1/8, and 2 at 1/4 (all before testing for convergence/beginning main loop) will improve performance in most cases. I'll have to try some tests.

[Guv]

Alphanos: I think the precision approximately doubles for each iteration, assuming you have at least 2-4 digits correct. It might converge even faster as you get closer.

For 1000 digits of precision, getting 64 correct and doing four more iterations should get about 1024 digits correct, which would assure that 1000 digits were okay.

Visualize a graph of the function y = x² - n

Square roots are intersections of this curve and the X-Axis.

The newton method is equivalent to the following geometric method.

• Pick a point (x₀, y₀) on the curve as the first estimate of the intersection.
• Draw a tangent at that point and determine where it crosses the X-Axis.
• Use the x-value at the interection as x₁ and determine the point (x₁, y₁) on the curve.
• Draw a tangent at that point and determine where it crosses the X-Axis, to get the next estimate for x.
• Keep this up until the value of y for some point is as close to zero as you require it to be.

The distance between the tangent intercept, and the curve intercept might approach zero faster than being halved with each iteration.

If you make an accurate drawing, the speed of convergence is vividly indicated.

[Alphanos]

Thanks