Number: 167
Date: 18-Mar-84 15':45':45
Submitter: Sannella.PA
Source: JonL.pa
Subject: Improved accuracy, speed for transcendental fns
Assigned To: Jonl.pa
Attn:
Status: Closed
In/By: Carol
Problem Type: Performance
Impact: Moderate
Difficulty:
Frequency:
Priority: Perhaps
System: Language Support
Subsystem: Arithmetic
Machine:
Disk:
Lisp Version:
Source Files:
Microcode Version:
Memory Size:
File Server:
Server Software Version:
Disposition: '
["Sannella.PA" "20-Aug-84 17':18':56" Attn': Status':(Fixed->Closed) In/By':]
Description: '
Date': 18 Mar 84 00':05 PST'
From': JonL.pa'
Subject': Lisp': New AARITH, and speed-up for some generic arith fns'
To': LispSupport.pa'
cc': JonL.pa, vanBuer%USC-ECL@SRI-NIC'
Lisp-System-Date': 7-Mar-84 13':56':45'
Machine-Type': Dorado'
'
I''ve incorporated some polynomial approximations suggested by Darrel vanBuer into the functions of AARITH, and put a substantially new file on <LispCore>Sources> -- the old linear-interpolation kludge is now gone, along with its infamous inaccuracy.'
'
Note': there was a bug in the former ARCTAN2 in that (ARCTAN2 -3 1) returned a value out of the range. Also the arguments X and Y were reversed, so I fixed them to accord with the manual.'
'
There is an undocumented function ATAN on AARITH which looks like someone''s botched attempt to copy the PDP-10 MacLisp version. I''ve supplied a correct translation (from the machine language, glaaag!), but I suspect that this one just ought to be deleted (Masterscope finds no system callers of ATAN). Aside from everything else, it appears to duplicate the functionality of ARCTAN2.'
'
I''ve also fixed a glitch in several of the generic arithmetic functions (ABS, PLUS, MINUS, and TIMES) whereby they forced the "hard case" of a microcoded floating point operation. By changing, e.g. (FDIFFERENCE 0 X) into (FDIFFERENCE 0.0 X), forms like (ABS 3.4) and (TIMES <fp-number1> <fp-number2>) were speed up by a factor of from 2 to 5. [This doesn''t apply to the DLion yet, since it doesn''t have floating-point in ucode].'
'
-----'
'
Date': 28 Mar 84 00':20 PST'
From': Masinter.pa'
Subject': new transcendental functions'
To': JonL'
cc': LispSupport'
'
What accuracy are they to? Are we within the range of what is possible given our floating point format?'
'
I tried in Interlisp-10, Interlisp-D (on a Dorado), and Franz on GSLVAX, (ARCSIN (SIN x)) and (ARCCOS (COS x)) for x = .001 and .027.'
'
fn maxc vax D'
sin .001 .001 .0009999999999999999 .009972'
cos .001 0.0 .001000000000000064311 .02799 (!!!!)'
sin .027 .027 .027 .026993'
cos .027 .02797 .02700000000000001 .044 (!!!!)'
'
I think there may be a problem with ARCCOS or COS.'
'
Is 32-bit IEEE floating point really that bad? (I guess Maxc is 36-bit, and the VAX defaults to double precision.)'
'
-----'
'
Date': 30 MAR 84 23':44 PST'
From': JONL.PA'
Subject': Re': new transcendental functions'
To': Masinter'
cc': LispSupport, JONL'
'
In response to your message sent 28 Mar 84 00':20 PST'
'
First off, you misplaced a decimal point on the first reported Interlisp-D'
result, so that should relieve a lot of anxiety.'
'
Now, the problems you point out have nothing to do, per se, with the "new"'
transcendental functions -- the "old" ones are much worse. In essence, yes, '
32-bit IEEE isn''t very good for this test -- (ARCCOS (COS X)) -- which '
is numerically unstable to a very high degree near zero.'
'
I might point out that this standard test is usually done with the angles '
given in radians -- thus you should be comparing the various values of '
X1 = (ARCCOS (COS .001 T) T)'
X2 = (ARCCOS (COS X1 T) T)'
. . . '
X(n+1) = (ARCCOS (COS Xn T) T)'
If you make this comparison, things don''t look quite so drastically bad; by'
allowing the .001 to be taken in degrees, you get the effect of starting'
out at about 1.7e-5 in radians. And for a problem *this* unstable, that'
makes all the difference.'
'
'
'
Finally, I''m a little surprised that the relative error for (ARCSIN (SIN .001))'
is so large (nearly 0.3%) -- ARCSIN isn''t that unstable in this range, '
it''s ARCCOS that''s the loser -- so I did a little mathematical research '
and came up with another approximation that''s significantly better near '
the endpoint of the interval [it''s also less "consy"!] with a relative '
error back down near 1e-7. This version, along with a slight improvement '
to SIN for arguments less than PI/4, is in the new AARITH on <LispCOre>Sources>.'
'
-----'
'
From': MASINTER.PA'
Date': 31-Mar-84 12':40':19 PST'
Subject': Re': new transcendental functions'
In-reply-to': JONL''s message of 30 MAR 84 23':44 PST'
To': JONL'
cc': Masinter, LispSupport'
'
Thank you for your research. I realized after I sent in the message that ARCCOS of COS couldn''t be particularly stable for values near 0. '
'
Kelly Roach has a whole bunch of functions in his MATHFNS package. Can you document HORNERIFY enough to allow him to write his own versions?'
'
-----'
'
Date': 2 APR 84 12':37 PST'
From': JONL.PA'
Subject': Re': new transcendental functions'
To': MASINTER'
cc': JONL, Roach'
'
In response to your message sent 31-Mar-84 12':40':19 PST'
'
Out of curiosity, how did you "realize that ARCCOS of COS couldn''t be'
particularly stable for values near 0" after you sent out the msg? After'
all, it could well have been a bug -- the situation with ARCSIN of SIN,'
for example, was similar! and it takes a little analysis to determine'
which is the Poor-Approximation and which is the Inherent-Instability.'
'
'
HORNERIFY is merely a shorthand way of writing out a polynomial evaluation --'
it exists as a non-exported macro in AARITH and takes a form like'
(HORNERIFY a[n] X a[n-1] . . . a2 a1 a0)'
into the appropriate composition of FPLUS and FTIMES to evaluate the polynomial'
a[n]X↑n + a[n-1]X↑(n-1) + . . . + a1X + a0'
using, of course, horners rule for order of calculations.'
'
'
The brunt of the work for numerical approximation is to find a transformation'
of the argument domain into some limited domain under which a polynomial'
approximation has desirable "closeness" to the function. There exist methods'
for generating polynomials with minimum "max error" over the range, and there'
exist methods for generating polynomials with minimum integral of error over'
the domain. As you can imagine, a few points with large error can still be'
quite acceptable to the latter method.'
'
'
Workaround:
Test Case:
Edit-By: Sannella.PA
Edit-Date: 20-Aug-84 17':19':00