Google Groups no longer supports new Usenet posts or subscriptions. Historical content remains viewable.
Dismiss

algorithm for computing sqrt()

36 views
Skip to first unread message

Aaron Wallace

unread,
Jan 30, 2002, 7:56:06 PM1/30/02
to
Hi everyone,

I know there is a way for calculating sqrt() by recursion:

y(n) = .5 * [ x(n)/y(n-1) + y(n-1) ]

where y(n) -> sqrt(A) as n -> inf, if x(n)=A and y(-1)=1.

I tried it on my HP48 and it takes 6 iterations to get to the answer 3
with A=9. I'm thinking there has to be a faster way than this. Does
anyone know if the HP internally employs this method of recursion or
does it use some other method? Does anyone know of a better one?

Just to see how fast the sqrt() is on the HP48, I timed %SQRT on the
argument 9. The TIM program reports 1.4ms. I'd be awfully impressed if
the 48 did 6 iterations of the above expression in 1.4ms.

--
Aaron

Ville Koskinen

unread,
Jan 31, 2002, 5:31:16 AM1/31/02
to
Aaron Wallace wrote:
> I know there is a way for calculating sqrt() by recursion:
...

> I tried it on my HP48 and it takes 6 iterations to get to the answer 3
> with A=9. I'm thinking there has to be a faster way than this. Does
> anyone know if the HP internally employs this method of recursion or
> does it use some other method? Does anyone know of a better one?

I used to wonder this too.

A good initial guess y(-1) would speed the algorithm up. Intuitively it
seems that you could choose or interpolate your guess from a table of
precalculated square roots. This might not be too fast, though.

On the other hand, if you're doing arbitrary precision floating point,
you could use the fixed precision value as a guess. I tried this once
and I was surprised how fast the algorithm converged.

But I don't think calculators use this algorithm. There is a paper on
this on the net... but I've forgotten where.

Regards,
Ville Koskinen

Peter Geelhoed

unread,
Jan 31, 2002, 6:08:54 AM1/31/02
to
Aaron Wallace wrote:

> Hi everyone,
>
> I know there is a way for calculating sqrt() by recursion:
>
> y(n) = .5 * [ x(n)/y(n-1) + y(n-1) ]

Ah yes, the square halving algorithm! Perhaps my favorite algoritm.
It is simple and fast.

There also is an algoritm that computes digits (like the miraculous
Bailey-Borwein-Plouffe
Pi Algorithm), which I *actually* learned in primary school in the mid 80's

(before the calculator age :-)
take a look at this site
http://www.interstice.com/~sdattalo/technical/theory/sqrt.html

--
This message was written with 100% recycled electrons

Pivo

Tanguy Briancon

unread,
Jan 31, 2002, 7:48:40 AM1/31/02
to
Ville Koskinen wrote:

If x_n is the sequence x_(n+1)=1/2(x_n/A+x_n), you have:
|x_n(+1)-sqrt(A)| < C|x_n-sqrt(A)|^2
so at each step of the algorithm the number of good digits is multiplying
by 2.

Peter Geelhoed

unread,
Jan 31, 2002, 9:59:49 AM1/31/02
to
Tanguy Briancon wrote:

> If x_n is the sequence x_(n+1)=1/2(x_n/A+x_n), you have:
> |x_n(+1)-sqrt(A)| < C|x_n-sqrt(A)|^2

let's make some definitions:
D=X(0)-sqrt(A)
E=X(1)-sqrt(A)=(X(0)^2+A)/(2*X(0))-sqrt(A))

and call X(0) X for short

so the error will decrease if D^2 is larger than E^2

so
D^2-E^2 > 0

(X-sqrt(A))^2-(X(0)^2+A)/(2*X(0))-sqrt(A))^2 > 0

and after some rewriting

(X+sqrt(A))*(3*X-sqrt(A))*(sqrt(A)-X)^2 > 0

which means that X should be larger than sqrt(A)/3
for the error to decrease.

Scott Hemphill

unread,
Jan 31, 2002, 12:12:45 PM1/31/02
to
Peter Geelhoed <p.f.ge...@student.tnw.tudelft.nl> writes:

> Aaron Wallace wrote:
>
> > Hi everyone,
> >
> > I know there is a way for calculating sqrt() by recursion:
> >
> > y(n) = .5 * [ x(n)/y(n-1) + y(n-1) ]
>
> Ah yes, the square halving algorithm! Perhaps my favorite algoritm.
> It is simple and fast.
>
> There also is an algoritm that computes digits (like the miraculous
> Bailey-Borwein-Plouffe
> Pi Algorithm), which I *actually* learned in primary school in the mid 80's
>
> (before the calculator age :-)
> take a look at this site
> http://www.interstice.com/~sdattalo/technical/theory/sqrt.html

For architectures that don't have native division, or have slow native
division, it is easier to calculate the inverse square root of a number,
and then multiply by that by the number. The following sequence converges
to n^(-1/2) without any divisions:

y(k+1) = .5 * y(k) * (3 - n*y(k)*y(k))

This is how square roots of million-digit-long numbers is done.

Scott
--
Scott Hemphill hemp...@alumni.caltech.edu
"This isn't flying. This is falling, with style." -- Buzz Lightyear

Aaron Wallace

unread,
Jan 31, 2002, 1:52:20 PM1/31/02
to
Just wanted to thank everyone for their input. From what I've read, it
doesn't sound like I'll be speeding up this operation much. ;-)

--
Aaron

Dan Kalish

unread,
Jan 31, 2002, 9:11:19 PM1/31/02
to
Scott Hemphill <hemp...@hemphills.net> wrote in message news:<m3elk6a...@pearl.local>...
> Peter Geelhoed <p.f.ge...@student.tnw.tudelft.nl> writes:
>

> For architectures that don't have native division, or have slow native
> division, it is easier to calculate the inverse square root of a number,
> and then multiply by that by the number. The following sequence converges
> to n^(-1/2) without any divisions:
>
> y(k+1) = .5 * y(k) * (3 - n*y(k)*y(k))
>
> This is how square roots of million-digit-long numbers is done.
>

Neat!

BTW, Aaron, my HP48 took more than 6 iterations to get SQRT(9) to 12
decimal places using the half algorithm. I don't have a clue how the
hardware does it.

Dan

Aaron Wallace

unread,
Feb 1, 2002, 3:25:38 PM2/1/02
to
Hi Dan,

What did you use as your initial guess? I used 1. Maybe also we had different modes set and
mine rounded to 3 sooner or something. Also, there is a distinct possibility that I simply can
not count anymore ;-)

--
Aaron

Robert Tiismus

unread,
Feb 1, 2002, 5:31:12 PM2/1/02
to
Aaron Wallace wrote:

> I tried it on my HP48 and it takes 6 iterations to get to the answer 3
> with A=9. I'm thinking there has to be a faster way than this. Does
> anyone know if the HP internally employs this method of recursion or
> does it use some other method? Does anyone know of a better one?

I have read somewhere that HP employs the CORDIC algorithm to calculate
rational functions,
sin, cos, tan, exp, ln, sqrt. It's close relative to algorithm whitch is
used in floating point units of x86 processors. One old Dr. Dobbs journal
explains the algorithm in details. Sad thing is that to use the algorithm,
one has to have tables of arctan(2^-n), and one other function, at hand :(

Best wishes,

Robert Tiismus

Aaron Wallace

unread,
Feb 2, 2002, 1:11:09 PM2/2/02
to
Actually, I knew that the HP used the cordic algorithm to compute sines and
cosines and other trig functions, plus log and ln, but I didnt know it could
be used for square root too. Apparently, the "HP Journal" talks about this,
but I've never read it.

Tom Ferrio

unread,
Feb 3, 2002, 10:21:42 AM2/3/02
to
An article about CORDIC techniques with a nice list of original references
can be found at
http://education.ti.com/product/tech/30xa/faqs/faq83086.html.

The description likely applies to HP products as well.

Best regards, Tom Ferrio
"Aaron Wallace" <awal...@utdallas.edu> wrote in message
news:3C5C2CC4...@utdallas.edu...

Aaron Wallace

unread,
Feb 3, 2002, 11:16:26 AM2/3/02
to
Thanks for the reference. Very interesting.

Werner Huysegoms

unread,
Feb 4, 2002, 2:15:34 AM2/4/02
to
Robert Tiismus <Robert....@mail.ee> wrote in message > I have read somewhere that HP employs the CORDIC algorithm to calculate

> rational functions,
> sin, cos, tan, exp, ln, sqrt. It's close relative to algorithm whitch is
> used in floating point units of x86 processors. One old Dr. Dobbs journal
> explains the algorithm in details. Sad thing is that to use the algorithm,
> one has to have tables of arctan(2^-n), and one other function, at hand :(
>
> Best wishes,
>
> Robert Tiismus

Hi Robert.
The CORDIC algorithms are used for the trig and exp functions alright,
but not for sqrt.
The sqrt algorithm implemented is the decimal equivalent of an
algorithm that
finds each bit of the answer successively.
Barring a few details, it goes like this.
Suppose you want to compute sqrt(a), a=2, and you have an estimate x=1
and the quantity d = a-x*x = 1.
Then the next digit of x is found adding 1 to the next decimal
position, and
inspecting d:
if x'=x+10^-i
then d'=a-x'*x'=d-(x+x')*10^-i

To make the computation of d even easier, it is multiplied by 10
every time a new decimal digit is to be found:

x d
1.00000000 10.0000000000
1.10000000 7.9 (10 - 1 - 1.1)
1.20000000 5.6 (7.9 - 1.1 - 1.2)
1.30000000 3.1
1.40000000 0.4
1.5 <0

We added one too many, so we subtract it again, multiply d by 10 and
start with the next digit:

x d
1.4 4
1.41 1.19
1.42 <0

x d
1.41 11.9
1.411 9.079
1.412 6.256
1.413 3.431
1.414 0.604
1.415 <0

...

You get the idea.
The process takes only slightly longer than a multiplication.

Best Regards,
Werner Huysegoms

Wolfgang Rautenberg

unread,
Feb 4, 2002, 6:53:39 AM2/4/02
to
Werner Huysegoms wrote:
> The CORDIC algorithms are used for the trig and exp functions alright,
> but not for sqrt.
> The sqrt algorithm implemented is the decimal equivalent of an
> algorithm that
> finds each bit of the answer successively...

According to your description, essentially the same
algorithm we still learned in the school :-) This is
also the way the Babylonians computed sqr(2) long before
Greek mathematics was at its height - and quickly fell
into a deep crisis because it wasn't understood why the
real world is as it it and that the diagonal of the
unit square isn't "rational". Babylonians often used

1 + 25/60 + 51/60^2 + 51/60^3

as an approximation for sqrt(2) which, recalculated
as a decimal fraction, is exact up to the 6th decimal.
They needed this precision not only in land measuring
but also in their highly developped astronomy.

The standard algorithm described by Werner was known
also in indian and arabian mathematics. It came to Europe
only due to the "Liber abbaci" by Leonardo from Pisa, also
called Fobinacci. The book was first published in 1202. In
this basic work also the modern algorithms for adding and
multiplying natural numbers were presented, all in Latin.
Leonardo was a protegee of the Roman-German Emperor
Friedrich II and hence had enough time to do math ...

- Wolfgang

PS - Question. Is sqrt(2) first computed as a binary
fraction, i.e., as 1.0101000100101...b and then
transformed into a decimal fraction?

Wolfgang Rautenberg

unread,
Feb 4, 2002, 7:08:55 AM2/4/02
to
Sorry, Fibonacci, of course. The same who also
defined his famous sequence. Did you know that
Fibonacci also was the first who described an
iterative procedure to get the cubic root of a
real number? He explained his procedure with a
single example, XROOT(3,47), but his decription
is instructive enough to get the general idea.
His procedure converges linearly and is as good
as the latest school procedure :-)

- Wolfgang

Aaron Wallace

unread,
Feb 4, 2002, 3:08:13 PM2/4/02
to
Hi Werner,

It appears that cordic can be used to compute a square root. I looked in the HP's rom and saw what you are describing and one of the
links posted previously describes the binary "digit-by-digit" process that the HP uses. The same process can be used with base 10,
which is how most students learn to compute square roots and it sounds like the same process that Wolfgang was talking about. It's
interesting to see all the various methods.

--
Aaron

Werner Huysegoms

unread,
Feb 5, 2002, 1:13:23 AM2/5/02
to
Hi Aaron

> It appears that cordic can be used to compute a square root.

Yes, but not in the 48/49 ROM

> I looked in the HP's rom and saw what you are describing and one of the
> links posted previously describes the binary "digit-by-digit" process that the HP uses. The same process can be used with base 10,

Well, that's what I said.. The CORDIC implementation would probably be
a lot slower as it typically needs double precision to arrive at
correctly rounded
12-digit values for all arguments. The code as it is is very well
written;
I cannot improve it.

Cheers, Werner

Werner Huysegoms

unread,
Feb 5, 2002, 1:16:55 AM2/5/02
to
Hello, Wolfgang!

> PS - Question. Is sqrt(2) first computed as a binary
> fraction, i.e., as 1.0101000100101...b and then
> transformed into a decimal fraction?

No, no! all floating point computations are carried out in BCD.

Werner

John H Meyers

unread,
Feb 5, 2002, 8:15:41 PM2/5/02
to Werner Huysegoms
What ever happened to the "odd integer" method,
which was in use starting with the HP35?

That method is like ordinary division, except that
you increase the "current divisor" by 2 before
each subtraction; it works because 1+3+5+...(2*n-1)=n^2

The few extra speed tricks used in the HP35 were described
in the very first HP Journal article in June 1972 about
"The Powerful Pocketful"

Credit for the HP35 algorithms was given to HPs then current
corporate "cash manager" -- and ever since that time,
bean counters have continued to rule the roost :)

http://www.hpmuseum.org/journals/journals.htm
http://www.vcalc.net/hp-jrnl.htm

[r->] [OFF]


-----= Posted via Newsfeeds.Com, Uncensored Usenet News =-----
http://www.newsfeeds.com - The #1 Newsgroup Service in the World!
-----== Over 80,000 Newsgroups - 16 Different Servers! =-----

Dan Kalish

unread,
Feb 5, 2002, 9:55:12 PM2/5/02
to
Hi, Aaron:

6 iterations it is.

There are three kinds of mathematicians: those who can count and those
who can't. <g>

Aaron Wallace <awal...@utdallas.edu> wrote in message news:<3C5AFAC8...@utdallas.edu>...


> Hi Dan,
>
> What did you use as your initial guess? I used 1. Maybe also we had different modes set and
> mine rounded to 3 sooner or something. Also, there is a distinct possibility that I simply can
> not count anymore ;-)
> >

Werner Huysegoms

unread,
Feb 6, 2002, 12:36:13 AM2/6/02
to
John H Meyers <jhme...@miu.edu> wrote in message news:<3C6083B...@miu.edu>...

> What ever happened to the "odd integer" method,
> which was in use starting with the HP35?
>
> That method is like ordinary division, except that
> you increase the "current divisor" by 2 before
> each subtraction; it works because 1+3+5+...(2*n-1)=n^2
>

Hi John!
It actually is still that. Notice in my example that the first time,
you subtract 2.1 (1.0+1.1), then 2.3 (1.1+1.2) etc.
Werner

Jonathan Busby

unread,
Feb 6, 2002, 3:58:29 AM2/6/02
to
On 3 Feb 2002 23:15:34 -0800, werner-h...@freegates.be (Werner
Huysegoms) wrote:

>Hi Robert.
>The CORDIC algorithms are used for the trig and exp functions alright,
>but not for sqrt.
>The sqrt algorithm implemented is the decimal equivalent of an
>algorithm that
>finds each bit of the answer successively.
>Barring a few details, it goes like this.
>Suppose you want to compute sqrt(a), a=2, and you have an estimate x=1
>and the quantity d = a-x*x = 1.
>Then the next digit of x is found adding 1 to the next decimal
>position, and
>inspecting d:
>if x'=x+10^-i
>then d'=a-x'*x'=d-(x+x')*10^-i
>
>To make the computation of d even easier, it is multiplied by 10
>every time a new decimal digit is to be found:

<snip>

Actually, I think it is done a bit differently in the 48/49 ROM.
( or maybe these are the details you left out ;)

Here's how I derived it :

Assuming y is the integer quantity whose square root is to be taken, x
is the guess of this root in successive iterations and R is the
remainder, then :

y - ( x + 10^n )^2 = y - x^2 - 2*x*10^n - 10^(2*n) =
R - 2*10^n*( x + 10^n ) + 10^(2*n) ( adding and subtracting 10^(2*n) )

Since the ratio of the remainder and x^2 is all that matters we can
divide by 2 :

R/2 - 10^n*( x + 10^n ) + 5*10^(2*n-1)

In the process of testing successive values of x's nth digit the above
can be rewritten :

R/2 - 10^n*( x + (a-1)*10^n + 5*10^(n-1) )

where a is the digit value.


The above final version of the expression is what the 48/49 ROM sqrt
routine is based on.


Here is an integer implementation of the above algorithm. It really
only differs from the floating point version in that the guess gets
successively shifted right and the remainder is not multiplied by ten
when starting on a new digit position in the guess. In the following
A.W is assumed to contain the argument ( < 1E14 ) . The result
( IP(SQRT(A.W)) ) is in C.W :

SETDEC
ASL W

C=A W
A=A+A W
A=A+A W
A=A+C W
ASR W

C=0 W
P= 13
LC(1) 5


- CSR WP
C=C-1 P

-- C=C+1 P
A=A-C W
GONC --

A=A+C W
CSR W

P=P-1
P=P-1
GONC -

SETHEX


-------------------------------------------------------------------------------

Jonathan Busby - <j...@SNMAPOhouston.rr.com>

Remove the random permutation of "NOSPAM" from my e-mail address
before replying.


Jonathan Busby

unread,
Feb 6, 2002, 4:04:24 AM2/6/02
to
On Sun, 3 Feb 2002 09:21:42 -0600, "Tom Ferrio" <tfe...@yahoo.com>
wrote:

>An article about CORDIC techniques with a nice list of original references
>can be found at
>http://education.ti.com/product/tech/30xa/faqs/faq83086.html.
>
>The description likely applies to HP products as well.
>
>Best regards, Tom Ferrio


Another good reference on CORDIC algorithms is :
http://www.fpga-guru.com/cordic.htm

John H Meyers

unread,
Feb 6, 2002, 5:18:26 AM2/6/02
to
More about the "odd integer" method for square roots:

Simplified explanation of the general practical algorithm:
(either of the following URLs should point to the same article)
http://groups.google.com/groups?oi=djq&selm=an_162734155 (Deja View :)
http://groups.google.com/groups?selm=4r2jn4%24oio%40news.iastate.edu

A "cheat" UserRPL method for HP48/49 user binary integers (#nnnnn):
http://groups.google.com/groups?selm=692an3%2441a%241%40news.iastate.edu

An old ML version for HP48 user binary integers, by Stefan Odelfalk:
http://groups.google.com/groups?selm=68vvbl%24e1u%241%40zingo.tninet.se
http://groups.google.com/groups?selm=69vb61%249d2%241%40cubacola.tninet.se
(should be adjusted to accept any "short" hex strings, which are legal)

Jonathan Busby's new ML for integer/hex values already in registers
looks simpler, however:
http://groups.google.com/groups?selm=4mn16uc2l02etion4...@4ax.com


But still no unlimited-precision square root of new 49G integers, right?

Jonathan Busby

unread,
Feb 6, 2002, 5:51:25 AM2/6/02
to
On Wed, 06 Feb 2002 04:18:26 -0600, John H Meyers <jhme...@miu.edu>
wrote:


>Jonathan Busby's new ML for integer/hex values already in registers
>looks simpler, however:
>http://groups.google.com/groups?selm=4mn16uc2l02etion4...@4ax.com

Shameless plug ;) : This is for BCD values of course. If you want a
fast sqrt routine based on the similar binary algorithm try
http://www.hpcalc.org/hp48/math/numeric/fisqrt.zip .
Hmmm... Looking at that code for the first time in 2 years I think I
can improve it ;). Here's the new version which works for hex strings
up to 16 nibbles in length :

::
0LASTOWDOB!
CK1NOLASTWD
CK&DISPATCH1
#B
CODE
POPHXS EQU #53F8D ( 48 only )
GOSBVL POPHXS
C=0 W
P= 15
LC(1) 4
B=0 W
SB=0
- B=B+C W
A=A-B W
GONC +
A=A+B W
B=B-C W
B=B-C W
+ B=B+C W
BSRB
CSRB
CSRB
?SB=0
GOYES -
A=B W
GOVLNG =PUSHhxsLoop
ENDCODE

;

Werner Huysegoms

unread,
Feb 6, 2002, 10:40:05 AM2/6/02
to
Hi Jonathan!
That's precisely what I meant with 'barring a few details'.
The algorithm is essentially the same, but I didn't want to go into
the gory details - that I studied in detail. And could not improve..

Werner

Robert Tiismus

unread,
Feb 8, 2002, 10:22:52 AM2/8/02
to
Werner Huysegoms wrote:

> Hi Robert.
> The CORDIC algorithms are used for the trig and exp functions alright,
> but not for sqrt.

Thank you for this information. As I understand, the cordic can be used to calculate sqrt, but there are better ways to do it as HP
is showing it by not using cordic?

Best wishes,

Robert Tiismus.

John H Meyers

unread,
Feb 27, 2002, 1:50:37 PM2/27/02
to
I wrote:

> But still no unlimited-precision square root of new 49G integers, right?

Is there a built-in 49G function
for the exact integer square root of an exact integer of any size?

There doesn't seem to be, for as of 1.19-6, ^ZINTSQRT seems interested
only in perfect squares (and otherwise factors the input, sort of),
while ^ZSQRT seems to rely upon an "extended real" approximation,
which is useful only for a limited argument size.

If not hiding there in ROM someplace,
is there then a program for the desired general function?

All right, then go away and try to make one, as a Mini-Challenge,
just to see how easy it is, or try this "Newton's method" kluge
(must be entered or downloaded in *exact* mode)

\<< R\->I \<< :^IDIV2: #3E9006h FLASHEVAL DROP \>>
\-> q \<< DUP -3 CF @ next 2 lines are commented
{ @ DUP SIZE R\->I 2 q EVAL 12 - 0 MAX OVER 100 PICK3 ^ q EVAL
@ I\->R \v/ CEIL R\->I 10 ROT ^ *
DUP :^ZSQRT: #E0006h FLASHEVAL DROP 1 + @ replaces commented lines
DO DUP2 q EVAL DUP ROT + 2 q EVAL SWAP
UNTIL OVER - DUP SIZE == END NIP } IFT \>> \>> 'ISQRT' STO

For example, try 5E100 as an argument, then do 5 * (in exact mode),
and see that the "golden ratio" 5*(1+sqrt(5))/10 must be
1.61803398874989484820458683436563811772030917980576...

For all I know, this program may still have
an obscure bug, caused by truncation in ^ZSQRT
(note that CEIL was necessary in my original version, where
IP wouldn't suffice, and that 1 + is necessary after ^ZSQRT),
but the idea in today's software factories is to
"get it out the door," and see whether customers
will buy it anyway, for that's where the "bottom line" is ;-)

[r->] [OFF]

Steen Schmidt

unread,
Feb 27, 2002, 6:43:22 PM2/27/02
to
> Is there a built-in 49G function
> for the exact integer square root of an exact integer of any size?
[...]

> ^ZSQRT seems to rely upon an "extended real" approximation,
> which is useful only for a limited argument size.

You're right, ^ZSQRT isn't good at very large arguments. I guess you must
program your way out of it.

This is from my memory, but there may be an easier solution after all. If I
find it, I'll post it here.

Regards
Steen

0 new messages