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
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
> 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
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.
> 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.
> 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
> 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
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
> 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
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...
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
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?
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
> 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
> 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
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! =-----
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 ;-)
> >
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
>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.
>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
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'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
> 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.
> 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]
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