REXX math functions
John Brock
been reading about them here (but I'm new), and I don't see them
in the FAQ file.
I had originally thought that math functions (sqrt, log, sin, etc)
had been left out of REXX because it was too difficult to implement
them with REXX style infinite precision. Then I was asked to
write such a package, and I discovered it was not that difficult
after all. I wrote a CMS exec, entirely in REXX, which included
all the usual pocket calculator math functions, was reasonably
efficient (and could be compiled to increase efficiency by a factor
of 10 to 40), and returned all answers to arbitrary precision.
It can, for example, calculate pi to 1000 places in about 20 seconds
(VIRTCPU = 16 sec, for CMS people), and of course is much faster with
more reasonable precisions. A typical call might look like this:
x = 'RXMATH'(25, 'POW', 5, 3.8)
i.e., 5 to the power of 3.8 to 25 places (= 452.9872897985597071691338).
The reason I am writing is I just want to put in a plug for this
format for built-in REXX math functions. (Of course a built-in
math function would not need a precision argument, as it would know
the current numeric digits setting). I think it would be a
mistake for there to be a separate log function and square root
function and so on -- it just clutters up the name space and makes it
harder to add new math functions later. I think it would be much more
elegant to put all math functions in a single REXX function and
use one of the arguments to specify the function required. You
would of course have to depart from the REXX convention of single
letter arguments, but I think it would be worth it.
I really didn't do anything terribly clever writing RXMATH. I
took the specified precision, doubled it to guarantee that there would
be no rounding errors, then just calculated the standard infinite
series until the sum no longer changed when new elements were added
(because they had gotten too small), and handed back the result
rounded to the requested precision. As I said, this is reasonably
efficient (and extremely accurate), and I can see no reason why this
couldn't be part of the REXX interpreter.
I actually did do one clever thing people might find interesting.
I have a "calc" function where you could do something like this:
x = 'RXMATH'(25, 'CALC', 'sqrt(2 * sin(ln(2)) + 1)')
(the result is 1.509278818716829936251786). The third argument is
a string which is interpreted, can contain any REXX function
(including external function calls), and in which the REXX math
functions are evaluated directly. The clever thing was that I still
doubled the precision at the beginning, and kept all intermediate
results at the doubled precision, and only rounded down when I
returned the result. This does wonderful things for the accuracy
of complicated expressions. (When the third argument here is omitted
my exec goes into an interpreted "REXX shell", but I don't think
we would want that as standard REXX).
If anyone is interested in my code let me know. I would have to
get permission from someone here to post it to the net, but that
would probably be no problem. There are actually two execs, one
of which can be compiled and the other handles the interpretive
"calc" function. They will run under either CMS or TSO. I don't
know about UNIX, as we have no UNIX REXX interpreter here.
I am sure people have already thought about this. I have seen
some math packages, but none that were completely general. For
completeness, my functions are: fact(x), perm(x,y), comb(x,y),
sqrt(x), pow(x,y), log(x,y) [log of y to any base x], exp(x), ln(x),
pi(), sin(x), cos(x), tan(x), cot(x), sec(x), csc(x), arcsin(x),
arccos(x), arctan(x), arccot(x), arcsec(x), arccsc(x), calc(str).
RXMATH will print a result if called from the command line. If
called as a function and given bad arguments it returns nothing,
hopefully causing your exec to blow up, the way a built-in
function with bad arguments does.
I can be reached at the addresses below. I can send mail to Usenet
addresses but not to Internet addresses.
--
John Brock (212) 909-4662
The First Boston Corporation
12 East 49th Street
New York, NY 10017
uunet!csfb1!jbrock
jbr...@csfb1.fir.fbc.com
John Churchman
Please send me copies of your math function EXECs. Thanks....
Dave Gomberg
>
> numeric digits 1000
> numeric fuzz 1
> call time r
> parse value 1 1/sqrt(2) 1/4 1 0 with a b c x p
> do j=1 until p=p1
> parse value (a+b)/2 sqrt(a*b) c-x*(b-a)*(b-a)/4 x+x with a b c x
> parse value p (a+b)*(a+b)/(4*c) with p1 p
> end
> say p
> say time(e)
Wow, is that ever clever code. But on my machine it converges to 3.056 and
change at a couple of precisions. For accuracy I would prefer:
NUMERIC DIGITS 1000; RETURN 22/7
which is even faster, and on my machine, more accurate. Maybe there is a
typo above, Ian? ;^)> Dave
Dave Gomberg GOMBERG@UCSFVM Internet node UCSFVM.UCSF.EDU (415)731-7793
Seven Gateview Court, San Francisco CA 94116-1941
Ian Collier
>Wow, is that ever clever code. But on my machine it converges to 3.056 and
>change at a couple of precisions. For accuracy I would prefer:
>NUMERIC DIGITS 1000; RETURN 22/7
The program, as copied from your quotation of my article, worked fine for
me. I'm sending you a trace. Perhaps other people should check their
implementations with this program as well...
Of course, it may just be an operator precedence problem with REXX/imc, or a
sqrt() problem with your interpreter.
Ian Collier
Ian.C...@prg.ox.ac.uk | i...@ecs.ox.ac.uk
Bill Salisbury -- C2390S
Manager, Admin. Computing Center
UMKC SSB 102F -- 235-1482
The program worked fine for me on an IBM PS/2 Model 50 using Personal REXX.
I finally had to interrupt its execution for 1000 digits -- it was taking
too long (a relative judgment).
My timing tests (for what they're worth) are:
10 digit precision -- 6.76 seconds-- 3.141592656
20 digit precision -- 40.15 seconds-- 3.14159265354941579892
25 digit precision -- 63.0 seconds-- 3.141592653594157988802306
I agree with the assessment of WOW to the coding technique.
HAVE A GOOD DAY!
John Brock
for REXX. I've gotten permission to post it and will do so
today. The whole thing is about 1000 lines (two execs and
one help file) and about 32K in size. This isn't big enough
to break my News software, but I really don't know much
about these things, so if you see this note but not the
main post let me know, and I'll try to mail it to you
directly.
In addition, Ian Collier has raised some questions, which
I will address here.
My original post had three purposes: to request information
about the current state of REXX math functions, to argue
that arbitrary precision math functions should be built in
to the interpreter and that they should take a certain form,
and to let people know about my math package, which frankly
I am very proud of and would like to see other people using.
I think it is kind of silly to argue whether REXX ought to
have a built in math package. If it is not too difficult,
of course it should! One of the reasons FORTRAN is so
popular among scientists is the built in math functions.
Not everyone will use them, but then not everyone will ever
need to calculate 22 / 7 to 1000 places, and yet REXX lets
you do it if you need to. A math package for REXX that is not
arbitrary precision is simply out of keeping with the spirit
of REXX! As for efficiency, there is no reason someone
who is writing an interpreter should not use the hardware
datatypes when appropriate, as long as he switches over
transparently to the arbitrary precision code when greater
precision is requested. I didn't do this because I
wanted a pure REXX solution, and because I wanted to keep
the complexity of the program manageable.
As for the format, this is just a personal preference.
I think it makes sense to bundle all the math functions
into a single call simply because there are many people
who will never use them, and this will keep the functions
out of their way. At the very least, if math functions are
implemented individually, they should not be scattered
alphabetically throughout the manual along with the other
REXX functions! For the record, my ideal built in REXX
math function calls would look something like this:
x = math('POW', 5, 3.8) /* = 5 ** 3.8 */
y = math('LN', 3.5) /* = natural log of 3.5 */
z = math('PI') /* = 6 * arcsin(0.5) :-) */
RXMATH is optimized for accuracy, not for speed, and I am
sure there are things that could be done to speed it up.
However I do have a fairly strong math background (Columbia
physics ABD), and I paid close attention to the convergence
of the various series expansions I used. For example, if
RXMATH is asked for sqrt(1200000), what it actually calculates
is sqrt(1.2) * 1000, and similar adjustments are made for
other functions. I hard code values for pi and e, and
only recalculate them when the hard coded precision is
inadequate. And I found a better series expansion for
pi than the one listed in all the handbooks (which may
be interesting mathematically but is totally pathetic
convergence-wise). As a result the efficiency of RXMATH,
if not quite comparable to FORTRAN in a tight loop, is
still quite good, especially if the program is compiled.
(And EXECLOADING won't hurt if you're on VM).
It's been years since I benchmarked RXMATH, and I am not
going to redo it now. I will mention again what I noted
in my first post, that when I used RXMATH to calculate pi
to 1000 places it took 20 seconds elapsed and 16 seconds
VIRTCPU. This is the compiled version running on what
I believe to be the latest and greatest IBM big iron, so
this is probably as fast as it's going to go!
One the whole I think I have a very clean, convenient,
and exceptionally accurate REXX math package, which adheres
very closely to the "REXX spirit". I sincerely hope people
will find it useful.
David Jones
> My timing tests (for what they're worth) are:
>
> 10 digit precision -- 6.76 seconds-- 3.141592656
> 20 digit precision -- 40.15 seconds-- 3.14159265354941579892
> 25 digit precision -- 63.0 seconds-- 3.141592653594157988802306
>
Unfortunately, if these numbers are supposed to be approximations of pi,
then they're off.
pi = 3.1415926535897932384626433832795...
I'll stick to my assembly routines (which agree with the value of pi
generated by the Maple symbolic algebra package).
--
David Jones, 6730 Tooney Drive, Orleans, Ontario K1C 6R4 CANADA
email: d...@qpoint.ocunix.on.ca Fido: 1:163/109.8
AMIGA: Advanced Multimedia with Interactive Graphics and Audio
W. J. Metzger
>
>In article <REXXLIST%9210071...@UGA.CC.UGA.EDU> Bill Salisbury -- C2390S
><C23...@UMVMA.BITNET> writes:
>> My timing tests (for what they're worth) are:
>>
>> 10 digit precision -- 6.76 seconds-- 3.141592656
>> 20 digit precision -- 40.15 seconds-- 3.14159265354941579892
>> 25 digit precision -- 63.0 seconds-- 3.141592653594157988802306
>>
>
>Unfortunately, if these numbers are supposed to be approximations of pi,
>then they're off.
>
>pi = 3.1415926535897932384626433832795...
>
>I'll stick to my assembly routines (which agree with the value of pi
>generated by the Maple symbolic algebra package).
I just tried rxmath for pi and got the following
3.1415926535897932384626433832795029
Looks to me like it agrees with MAPLE
John E. Brock
What are you running RXMATH on? On my VM/CMS system RXMATH gives exactly
the same result as your Maple. In case anyone wants to check to more
places, on my machine:
rxmath(1000, 'PI') =
3.1415926535897932384626433832795028841971693993751058209749
445923078164062862089986280348253421170679821480865132823066
470938446095505822317253594081284811174502841027019385211055
596446229489549303819644288109756659334461284756482337867831
652712019091456485669234603486104543266482133936072602491412
737245870066063155881748815209209628292540917153643678925903
600113305305488204665213841469519415116094330572703657595919
530921861173819326117931051185480744623799627495673518857527
248912279381830119491298336733624406566430860213949463952247
371907021798609437027705392171762931767523846748184676694051
320005681271452635608277857713427577896091736371787214684409
012249534301465495853710507922796892589235420199561121290219
608640344181598136297747713099605187072113499999983729780499
510597317328160963185950244594553469083026425223082533446850
352619311881710100031378387528865875332083814206171776691473
035982534904287554687311595628638823537875937519577818577805
32171226806613001927876611195909216420199
John E. Brock
I posted a response to this a little while ago thinking that the
original post (Bill Salisbury's) refered to my RXMATH program.
Now I'm not sure. The thing is, in RXMATH I hard code a value for pi,
and don't need to recalculate it unless I am asked for more than 25 places
of precision, so all three of the above timing tests should have returned
the same time, which should be very short, since there is actually no
calculation involved. Do Mr. Salisbury's numbers refer to another
math package? My News reader has been screwed up, so I'm afraid I
don't really know what is going on right now. Sorry!
Bill Salisbury -- C2390S
Manager, Admin. Computing Center
UMKC SSB 102F -- 235-1482
my calculations was caused by a SQRT function using 12 digit precision for
all calculations and thus introducing rounding errors no matter how precise
the PI calculation specified -- my fault, not the algorithm's!
HAVE A GOOD DAY!
Ian Collier
>> 10 digit precision -- 6.76 seconds-- 3.141592656
>> 20 digit precision -- 40.15 seconds-- 3.14159265354941579892
>> 25 digit precision -- 63.0 seconds-- 3.141592653594157988802306
>Unfortunately, if these numbers are supposed to be approximations of pi,
>then they're off.
>I'll stick to my assembly routines (which agree with the value of pi
>generated by the Maple symbolic algebra package).
For the record, the above numbers seem to have been generated by the
algorithm which I posted earlier (copied below), and are probably due to the
use of a sqrt() routine which gives only about 10 digits of precision. I
have said before and will say again that "my" algorithm is known to be
correct. I didn't invent it - I heard it at a recreational maths lecture.
It has quadratic convergence: that is, the number of correct figures doubles
at each iteration. It is therefore much better than any series summation.
I strongly recommend that Mr Jones tries it out. I didn't post a sqrt
function last time. I will this time, and thus all implementations of REXX
are more-or-less guaranteed to return the correct value of pi.
After 8 iterations, pi = 3.141592653589793238462643383279502884197169399
375105820974944592307816406286208998628034825342117069
16 seconds elapsed
After 10 iterations, pi = 3.14159265358979323846264338327950288419716939
937510582097494459230781640628620899862803482534211706798214808651328230
664709384460955058223172535940812848111745028410270193852110555964462294
895493038196442881097566593344612847564823378678316527120190914564856692
346034861045432664821339360726024914127372458700660631558817488152092096
2829254091715364367892590360011330530548820466521384146951941511610
214 seconds elapsed
parse arg d /* required number of digits */
if d='' then d=20
d=d+1
numeric digits d
numeric fuzz 1
parse value 1 1/sqrt(2) 1/4 1 0 with a b c x p
do j=1 until p=p1
parse value (a+b)/2 sqrt(a*b) c-x*(b-a)*(b-a)/4 x+x with a b c x
parse value p (a+b)*(a+b)/(4*c) with p1 p
end
numeric digits d-1
say "After" j "iterations, pi =" format(p)
exit
sqrt: /* You need this if and only if there is no builtin
arbitrary-precision sqrt() function */
procedure expose d
parse arg a
if a=0 then b=0
else if a<0 then do
say "sqrt("a"): Incorrect call to routine"
return
end
else do
numeric digits d+1
numeric form scientific
parse value format(a,1,,,0) with before '.' after 'E' e
if pos('E',before)>0 then parse var before before 'E' e
a=before||after /* a is now pointless */
if e='' then e=0 /* e is the exponent */
if e//2 then e=e-1
else a='0'a /* exponent is odd, and decimal comes after 2 digits */
b='';b1=0;c=0
numeric digits d%2+2
do i=1 by 2 until length(b)=d+1
b1=trunc(b1)||substr(a,i,2,0)
/* find least x s.t. (c||x)*x>b1. Find (c||x-1)*(x-1). */
d2=c||1
d1=0
x=1
do while d2<=b1
d1=d2
d2=d2+(c||x)
x=x+1
d2=d2+x
end
x=x-1
b=b||x
b1=b1-d1
c=trunc((c||x)+x) /* never in exponential form */
end
b=left(b,1)'.'substr(b,2)'E'e/2
end
numeric digits d
return format(b)
Ian Collier
Ian.C...@prg.ox.ac.uk | i...@ecs.ox.ac.uk
Bill Salisbury -- C2390S
Manager, Admin. Computing Center
UMKC SSB 102F -- 235-1482
it might be of interest to mainframers to see how the timing on a PS/2 model
50 might compare. As I said earlier, Ian pointed out that my results seemed to
have been corrupted by a "not so precise" SQRT function -- absolute right on!
I'm sorry my "contribution" was the cause of shifting the original focus of
the discussion -- which I perceived to be discussion of variable precision
algorithms for REXX mathematical functions. Ian responded to my posting with
his correct analysis of my inaccurate results via the REXXLIST last week. I
assumed that everyone interested in the topic would have caught the gist of
that message and that no response from me was necessary. Apparently I was
wrong. My apologies are formally extended to Ian Collier for causing the
accuracy of the algorithm to be questioned. Maybe all of this fuss does
point out a need for people to more precisely identify sources of error in
numeric calculations -- I referenced a SQRT function that I have stored
previously that calculated with a precision of 12 digits -and- used it in a
routine that was designed to permit variable precision calculation. I err'd
in not considering the "imbedded" restriction of my SQRT routine. It seems
that this type of oversight is a possible consequence of "modular" programming
languages and is a lesson that we must sometimes be reminded of. I learned
something from this episode. I hope other people did too. But please, if
you don't like the accuracy of the results produced by my run, don't blame the
algorithm or Ian Collier -- blame me in my programming oversight.
HAVE A GOOD DAY!
Ian Collier
| sqrt: /* You need this if and only if there is no builtin
| arbitrary-precision sqrt() function */
...
| if e//2 then e=e-1
I did, of course, mean:
if e//2 ^= 0 then e=e-1
|Ian Collier
|Ian.C...@prg.ox.ac.uk | i...@ecs.ox.ac.uk
Otto Stolz
on Tue, 13 Oct 92 15:16:45 GMT Ian Collier said:
> [...] I didn't post a sqrt
> function last time. I will this time, [...]
> sqrt: /* You need this if and only if there is no builtin
> arbitrary-precision sqrt() function */
> [42 more lines of code suppressed]
Why all this ado about formating the radicand, parsing mantissa &
exponent out of it, fiddling with the numeric settings, etc?
For me, the following simple Newton iteration did the job:
sqrt: procedure /* variant 1 */
parse value arg(1) 1 0 with x r r1
do while r ^= r1
parse value r (x/r+r)/2 with r1 r
end
return r
It requires that numeric fuzz be at least 1. A sqrt for an arbitrary
caller would read thus (and take more CPU time, of course):
sqrt: procedure /* variant 2 */
parse numeric nd .
numeric digits nd+1
numeric fuzz 1
parse value arg(1) 1 0 with x r r1
do while r ^= r1
parse value r (x/r+r)/2 with r1 r
end
return r
Do not worry about the altered numeric settings. Excerpt from CALL
HELPREXX:
: NUMERIC settings (the DIGITS, FUZZ, and FORM of arithmetic operations,
: described in the NUMERIC instruction) are saved and are then restored
: on RETURN. A subroutine may therefore set the precision, etc., that it
: needs to use without affecting the caller.
My times to calculate pi with Ian Collier's original program and various
sqrt routines (specified as internal functions, each) on a Comparex 7/88
under VM/SP HPO and CMS Rel 5 were:
significant ----------------- virt CPU times ------------------
places (nd) my variant 1 my variant 2 Ian Collier's variant
10 0.04 sec 0.05 sec 0.16 sec
32 0.17 sec 0.24 sec 0.68 sec
101 2.39 sec 2.58 sec 7.20 sec
319 30.16 sec 56.42 sec 73.67 sec
1009 496.12 sec
Note that the shorter times vary considerably due to system load
variations. The times for 100 significant places (and more) of pi
apparently vary less than 10% between repeated measurements.
Happy programming,
Otto Stolz <RZO...@DKNKURZ1.Bitnet>
<RZO...@nyx.uni-konstanz.de>
Scott Ophof
...
>I think it is kind of silly to argue whether REXX ought to
>have a built in math package. If it is not too difficult,
>of course it should! One of the reasons FORTRAN is so
>popular among scientists is the built in math functions.
REXX is a small and complete language for its purpose, easy to learn
& understand, and probably why REXX is a success with many different
types of users (not only scientists). So I argue against adding
BIFs that aren't really necessary (portability might demand some,
though I hope it can be done differently, more in keeping with the
basic philosophy behind REXX's design).
>I think it makes sense to bundle all the math functions
>into a single call simply because there are many people
>who will never use them, and this will keep the functions
>out of their way.
Math functions as BIFs for the scientists means text processing
function, database functions, graphics functions, etc. should also
become BIFs, their respective user groups would argue...
Bundling a set of related functions into an external function
package, as you and others have done, keeps them together and still
be easily invokable by users. The concept of a *small*, basic set
of BIFs with the rest relegated to external function packages seems
to me as a very sensible way of not wasting resources, and *still*
be easily useable by the various types of users.
I could argue that LEFT() and RIGHT() aren't really needed; they
are subsets of SUBSTR(str,'L',n) resp. SUBSTR(str,'R',n), if "L"
resp. "R" were recognized as valid values for the starting position.
But as my own devil's advocate, readability & easy understandability
suffer more than would be gained by the reduction of language size.
This is imho not the case if the ARG instruction were to be dropped
(at least become undocumented); it's nothing but pure shorthand for
PARSE UPPER ARG (which is more readily understandable than ARG).
I'd rather see LOWER added as keyword than the ARG instruction
retained. Same goes for IBM-only FIND() vs. generic REXX WORDPOS().
I'd like to see IBM undocument FIND() in its implementations (and
move DIAG() & such to an external system function package).
I'd also like the DOS, hardware, windowing, and miscellaneous
functions in Personal REXX moved to external function packages.
(Dave, any comments here?)
BTW, I would very much appreciate a list of BIFs for the Amiga and
OS/2 implementations of REXX. But please do not post it to REXXlist
or comp.lang.rexx; thanks in advance. :-)
>(And EXECLOADING won't hurt if you're on VM).
On multi-user opsyses where the users can share in-mem interpreted/
compiled executables, EXECLOADing or the equiv. is *the* way to go.
I wish all opsyses had this sharing capability...
Regards.
$$/
John Brock
Math functions do *not* mean text processing and et cetera for everyone.
(This is the old slippery slope argument). Math is different. It is
universal and extremely well defined, and REXX already has partial
math functionality which cries out to be completed. None of these
things are true of your other examples. The full set of "calculator"
functions is quite small: powers and logs, trig functions and their
inverses, and the combinatorials, about 20 in all. On the other hand,
there is no simple universal way to put database or graphics functionality
into REXX, and any such functionality would be entirely new.
The reason for putting math functions into the language rather than a
function package is portability. If they were part of the language then
anyone writing a REXX program could always count on having them available,
and would feel safe using them. This isn't such a strong argument that
it would justify the enormous effort of putting database or graphics
primitives into REXX, but it's strong enough to justify math functions,
which simply aren't that hard to do.
I have to admit, it just bugs me that REXX lets you say x = 3.3 ** 2 but
not x = 3.3 ** 2.2 !
--
John Brock
uunet!csfb1!jbrock
jbr...@csfb1.fir.fbc.com
Ian Collier
>Hi gang,
Hi, Otto! :-)
>on Tue, 13 Oct 92 15:16:45 GMT Ian Collier said:
>> sqrt:
>Why all this ado about formating the radicand, parsing mantissa &
>exponent out of it, fiddling with the numeric settings, etc?
Basically, it's because I have always found that to be the best method. Not
that I have done any scientific tests (until now, that is)...
I apologise now for the rest of my article, to all those who aren't
interested in the gory details of the square root function. Those who are
interested can read on... ;-) I hope that includes everyone who has written
an arbitrary-precision sqrt function recently...
>For me, the following simple Newton iteration did the job:
[2 variants omitted]
[by the way, could you tell me what data the "parse numeric" instruction
parses? Sorry to ask such a dumb question...]
>Do not worry about the altered numeric settings. Excerpt from CALL
>HELPREXX:
>: NUMERIC settings (the DIGITS, FUZZ, and FORM of arithmetic operations,
>: described in the NUMERIC instruction) are saved and are then restored
>: on RETURN. A subroutine may therefore set the precision, etc., that it
>: needs to use without affecting the caller.
It may say that in CALL HELPREXX, but that's got nothing to do with my
implementation - which is broken ;-) I'll just make a note to fix that...
>My times to calculate pi with Ian Collier's original program and various
>sqrt routines (specified as internal functions, each) on a Comparex 7/88
>under VM/SP HPO and CMS Rel 5 were:
>significant ----------------- virt CPU times ------------------
>places (nd) my variant 1 my variant 2 Ian Collier's variant
> 10 0.04 sec 0.05 sec 0.16 sec
> 32 0.17 sec 0.24 sec 0.68 sec
> 101 2.39 sec 2.58 sec 7.20 sec
> 319 30.16 sec 56.42 sec 73.67 sec
> 1009 496.12 sec
Here are my measurements, on a Sun SPARCstation SLC. Times indicate the sum
of user+system CPU, as measured by the "set time" feature of csh. For
completeness, I have also included "imc2", which is the C version of my sqrt
variant, compiled into an external function package.
digits iterations imc Otto1 Otto2 imc2
50 6 5.4 3.4 3.5 1.2
100 7 15.6 14.6 14.9 4.0
300 8 110.2 164.8 166.3 36.4
500 9 308.2 506.1 507.0 111.1
Now it wasn't until I saw your impressive figures that I considered that
REXX370 gives your sqrt a big advantage - since most of the work for one
division can be done in a single assembly-language instruction (namely DP).
On most other systems it takes several lines of C code. Hence, on most
systems Newton's method loses on high-precision calculations.
The Newton-Kantorowich theorem guarantees that we can solve F(x)=0, where
F(x) = x**2-N, by iteration of the form
x.0 = max(1,N/2)
x.(n+1) = (x.n+n/x.n)/2 (n>=0)
provided N>0 and N\=2, and that the result will be in the interval
(0,2*x.0). It also gives that the convergence will be quadratic; more
precisely,
|x.n-y| / (x.(n-1) -y) ** 2 -> 1/(2*y) as n -> infinity
where y is the solution, i.e. sqrt(N).
The theorem does not guarantee anything if the initial guess is less than
sqrt(N/2), so starting the iteration at 1 is dubious even though it works.
Closer examination shows that if N>2 and x.0=1 then x.1=(1/N+N)/2 - which is
slightly greater that but approximately equal to N/2. Hence starting at 1
does work but is slightly slower (especially when N is small but greater
than 2).
Rearranging the above gives that e.n/s approaches (e.n/y)**2/2, where e.n is
the error in the nth term, i.e. |x.n-y|. Since e.n/y is a measure of "how
many digits of x.n are correct", it means that if the error is small, we can
expect to double the number of correct digits at each approximation.
However, if e.n/y is greater than sqrt(2) then the above formula would make
the iteration diverge if "->" were replaced by "approimately equal".
Instead we find that for large errors there is another formula, namely:
e.n <= e.(n-1)/2 (and the two sides are approximately equal for large e)
So for large errors the convergence is linear, not quadratic. This leads me
to the following observation:
-> For large values of N, you must separate the mantissa from the exponent,
or else the initial error (which is N/2-y, much larger than y) will be
so large as to make the linear portion of the convergence outweigh the
quadratic part.
You can either find the square root of the mantissa and deal with the
exponent separately, or simply use the exponent to make a better initial
guess. But not doing either will make the routine excessively slow for
large inputs.
My second observation is to calculate the complexity of each algorithm.
Division takes O(d*d) time, where d is the number of digits [Note 1]. The
important part of Newton's method is the division, and Newton will iterate
approximately log2(d) times when the initial estimate is small enough, so
Newton is O(d*d*log d).
On the other hand, "imc" does one iteration for each digit; the iteration
contains a number of constant-time operations plus an inner loop of at most
10 additions (each taking O(d) time). "imc" has a complexity of O(d*d) -
which is asymptotically faster, although the constant involved is likely to
be higher than for Newton's method.
Here are some timings (again on a SPARCstation, this time using time('E') to
time the results), which demonstrate these observations. Firstly, a test
which is unfriendly to Newton's method, using large numbers:
num='.'copies('3',d-1)
do i=1 to 20
s=sqrt(i||num'E'i)
end
digits imc Otto1 Otto2 imc2
25 4.92 7.47 8.02 0.25
50 11.06 26.74 27.89 0.62
100 28.68 105.08 107.11 1.48
200 84.04 420.20 426.22 4.16
and secondly, a test which is better suited to Newton's method:
do i=1 to 10 by .5
s=sqrt(i+num)
end
digits imc Otto1 Otto2 imc2
25 4.51 2.06 2.18 0.23
50 10.65 7.99 8.29 0.43
100 27.97 33.57 34.19 1.20
200 81.49 144.74 146.16 3.93
The large numbers in the first example give "imc" an advantage even at low
precision, while in the second test "imc" only wins out at about 75 digits.
Given that on a 370-like system division will run perhaps an order of
magnitude faster (relative to other operations) than on another system, it
may be that several thousand digits are required to give "imc" the advantage
for low input numbers on a 370-like system.
Finally, since I have noted that "imc" has complexity O(d*d) and so does
division, I found it interesting to compare the relative speeds of my
division operator and my sqrt function when both written in C. I ran the
following program and came out with surprising results (note that I do not
recommend the practice of calling time(r) with r unquoted...):
a=10/3; b=20/7
call time r
do 10;n=a/b;end
say "Division: "time(r)
do 10;n=a*b;end
say "Multiplication: "time(r)
do 10; call sqrt a;end
say "Sqrt: "time(r)
digits div mul sqrt
25 0.10 0.06 0.11
50 0.33 0.21 0.24
100 1.22 0.80 0.63
200 4.71 3.15 2.12
500 29.05 19.64 11.31
1000 116.40 78.80 45.65
This shows that, as long as division is not a machine-level operation, "imc2"
can carry out a square root in less time than it takes to do a division. It
would therefore be impossible for Newton's method to be better.
Ian Collier
Ian.C...@prg.ox.ac.uk | i...@ecs.ox.ac.uk
[Note 1] Non-trivial algorithms exist for multiplication and division which
have complexity between O(d) and O(d*d), but I do not consider
those. I would not expect that they are often used. Feel free to
correct me!
Ian Collier
Ian.C...@prg.ox.ac.uk | i...@ecs.ox.ac.uk
Otto Stolz
> sqrt:
In article <REXXLIST%92101616403467@DEARN>, I wrote:
> Why all this ado about formating the radicand, parsing mantissa &
> exponent out of it, fiddling with the numeric settings, etc?
On Wed, 21 Oct 92 14:41:57 GMT Ian Collier said:
> Basically, it's because I have always found that to be the best method.
...
> by the way, could you tell me what data the "parse numeric" instruction
> parses?
PARSE NUMERIC is CMS specific (sorry for posting system-dependent code;
I should have consulted the manuals first):
: The current numeric controls (as set by the NUMERIC instruction) are
: made available. These controls are in the order DIGITS FUZZ FORM.
: Example:
: 9 0 SCIENTIFIC
In article <REXXLIST%92101616403467@DEARN>, I wrote:
> Do not worry about the altered numeric settings. Excerpt from CALL
> HELPREXX:
...
On Wed, 21 Oct 92 14:41:57 GMT Ian Collier said:
> It may say that in CALL HELPREXX, but that's got nothing to do with my
> implementation -
I was referring to my own example procedures that set the numeric
settings but do not bother to reset them (which is indeed not necessary,
as REXX does it).
Ian's procedure does reset the settings; hence, the final result (arith-
metic expression in the return statement) is based on the original
settings.
On Wed, 21 Oct 92 14:41:57 GMT Ian Collier said:
...
Ian, thank you for the exhaustive explanation. I hope, tho term "ado"
did not hurt anybody (it was not meant that way).
Best wishes,
Otto Stolz <RZO...@DKNKURZ1.Bitnet>
<RZO...@nyx.uni-konstanz.de>
Ian Collier
[About square roots]
>Ian, thank you for the exhaustive explanation. I hope, tho term "ado"
>did not hurt anybody (it was not meant that way).
Incidentally, since I wrote that article it has come to my attention that
with Newton's method you do not have to use the same precision throughout
the operation. So if you rewrote your sqrt function to manipulate NUMERIC
DIGITS cleverly (so that k steps before the final result, the precision is
D*2**-k where D is the required precision - or something similar) then you
could in fact make its complexity equal to that of the divide instruction
(albeit multiplied by a fairly large constant). So then, although
digit-by-digit would always be faster on a machine without decimal-divide
instructions, it may well be advantageous to use Netwon's method on a
machine with decimal-divide. (Of course I have still avoided the issue of
time-less-than-d-squared divide operations).
Ian Collier
Ian.C...@prg.ox.ac.uk | i...@ecs.ox.ac.uk