Author 
Message 
lcw1964 science forum beginner
Joined: 16 Jul 2006
Posts: 2

Posted: Sun Jul 16, 2006 7:22 am Post subject:
Correction factor in computing exp()?



I have recently developed an interest in Cody's approximations of
erf() and erfc(). Here is a snippet of the Fortran code from
http://www.netlib.org/specfun/erf:
YSQ = AINT(Y*SIXTEN)/SIXTEN
DEL = (YYSQ)*(Y+YSQ)
RESULT = EXP(YSQ*YSQ) * EXP(DEL) * RESULT
The point of this is to simply compute exp(y*y) and multiply to the
result of a rational approximation calculation. Instead of doing it
straight ahead, the code multiplies y by 16, takes the integer part,
multiplies the 16 back, calls it YSQ, then computes some curious factor
DEL. Then the routine takes exp() of these separated pieces and
multiplies the whole thing together.
Does anyone know the computational point of doing this? Is or was there
some bug in the Fortran implementation of exp() that required some sort
of correction? Someone familiar with the GSL in C advises me that this
turns up in that code too. I have Cody's original 1969 paper on which
this code is based (indeed, Cody is the author of the code too), and
there is no mention of this little adjustment, or correction, or
whatever it is.
Many thanks in advance. This is driving me nuts with head
scratching....
Les 

Back to top 


lcw1964 science forum beginner
Joined: 16 Jul 2006
Posts: 2

Posted: Sun Jul 16, 2006 7:37 am Post subject:
Re: Correction factor in computing exp()?



lcw1964 wrote:
Quote: 
Many thanks in advance. This is driving me nuts with head
scratching....

I have even just written a little code in Maple the carries out the
computation to whatever precision I like, then compares it to a direct
calculation of exp(y*y). Gives the exact same answer for a wide range
of argumentssomething I could''ve readily concluded if I had just
worked thru the basic algebra, I gather, but I am a bit slothful.
But still, I don't know why the programmer would go to this trouble.
Evidently, it is intended to improve the accuracy or precision of the
exp() calculation.
Still eager for wisdom. Fill me, please!
Les 

Back to top 


dave_and_darla@Juno.com science forum beginner
Joined: 17 Oct 2005
Posts: 49

Posted: Sun Jul 16, 2006 3:42 pm Post subject:
Re: Correction factor in computing exp()?



lcw1964 wrote:
Here is a snippet of the Fortran code from
Quote:  http://www.netlib.org/specfun/erf:
YSQ = AINT(Y*SIXTEN)/SIXTEN
DEL = (YYSQ)*(Y+YSQ)
RESULT = EXP(YSQ*YSQ) * EXP(DEL) * RESULT
Does anyone know the computational point of doing this?

It enhances accuracy on computers that don't have guard digits. On such
a machine, typically, Y*Y will be accurate to about 1/2 ULP (unit in
the last place). But YSQ*YSQ+DEL should have more accuracy because
YSQ*YSQ and YYSQ can be computed exacly, and YYSQ will be smaller
than Y if abs(Y) > 1/16. Thus, the 1/2 ULP introduced in Y+YSQ will be
multiplied by a small number.
Example: Suppose the machine uses IEEE single precision and Y =
4.00xxx. Then YSQ will be 4.0 and YYSQ = 0.00xxx, which is less than
Y/400. Then DEL < (Y/400)*(2*Y) = Y*Y/200, and YSQ*YSQ+DEL will have 7
to 8 bits more accuracy than Y*Y. You can check this by comparing Y*Y 
YSQ*YSQ to DEL.
Dave 

Back to top 


Google


Back to top 



The time now is Tue Mar 26, 2019 3:55 am  All times are GMT

