You might think that you have worked around the tanh() library bug with
this code, but you have actually done more than that. The expression
you are using
(1 - aa*aa/15)*aa/3
actually consists of the leading terms for the taylor expansion of the
exact function that you are evaluating, f(x)=coth(x)-1/x. There is no
cancellation of error at all for small arguments in that taylor
expansion, where of course there is in the tanh() expression. Therefore,
when you switch at x=2.e-4 from the expansion to the original tanh()
expression, you are actually introducing errors, not eliminating them.
It is correct that the truncated taylor expansion has errors, and for
sufficiently large arguments, those errors will begin to dominate. So
one might ask, when should the switchover occur? You can figure it out
by hand, but perhaps an even better demonstration of this is the
following program:
program cotherr
use iso_fortran_env
implicit none
real(real32) :: x, f, g
real(real128) :: xp, fp
x = epsilon(x)
do while (x<=128.0)
f = 1/tanh(x) - 1/x
g = (1-x*x/15)*x/3
xp = x
fp = 1/tanh(xp) - 1/xp
write(*,'(5es12.4)') x, f, g, (f-fp), (g-fp)
x = x + x
end do
end program
I'm "cheating" here by using 128-bit floating point as a reference. Of
course, there are cancellation errors there too, but they are small
enough compared to the 32-bit arithmetic so that they can be ignored.
The output from this program on my computer is:
1.1921E-07 0.0000E+00 3.9736E-08 -3.9736E-08 1.1842E-15
2.3842E-07 0.0000E+00 7.9473E-08 -7.9473E-08 2.3685E-15
4.7684E-07 0.0000E+00 1.5895E-07 -1.5895E-07 4.7370E-15
9.5367E-07 0.0000E+00 3.1789E-07 -3.1789E-07 9.4739E-15
1.9073E-06 0.0000E+00 6.3578E-07 -6.3578E-07 1.8948E-14
3.8147E-06 0.0000E+00 1.2716E-06 -1.2716E-06 3.7897E-14
7.6294E-06 0.0000E+00 2.5431E-06 -2.5431E-06 7.5801E-14
1.5259E-05 0.0000E+00 5.0863E-06 -5.0863E-06 1.5166E-13
3.0518E-05 0.0000E+00 1.0173E-05 -1.0173E-05 3.0380E-13
6.1035E-05 0.0000E+00 2.0345E-05 -2.0345E-05 6.1138E-13
1.2207E-04 0.0000E+00 4.0690E-05 -4.0690E-05 1.2531E-12
2.4414E-04 0.0000E+00 8.1380E-05 -8.1380E-05 2.7487E-12
4.8828E-04 2.4414E-04 1.6276E-04 8.1380E-05 7.4376E-12
9.7656E-04 3.6621E-04 3.2552E-04 4.0690E-05 1.2935E-12
1.9531E-03 6.7139E-04 6.5104E-04 2.0345E-05 1.0348E-11
3.9062E-03 1.3123E-03 1.3021E-03 1.0174E-05 -3.3633E-11
7.8125E-03 2.6093E-03 2.6042E-03 5.0969E-06 -3.6280E-11
1.5625E-02 5.2109E-03 5.2082E-03 2.6279E-06 1.7395E-10
3.1250E-02 1.0414E-02 1.0416E-02 -1.8650E-06 4.1294E-10
6.2500E-02 2.0828E-02 2.0828E-02 3.3707E-07 -1.9348E-09
1.2500E-01 4.1623E-02 4.1623E-02 -2.1284E-07 -6.3824E-08
2.5000E-01 8.2988E-02 8.2986E-02 9.7103E-08 -2.0561E-06
5.0000E-01 1.6395E-01 1.6389E-01 -1.0945E-07 -6.4527E-05
1.0000E+00 3.1304E-01 3.1111E-01 -3.5789E-08 -1.9242E-03
2.0000E+00 5.3731E-01 4.8889E-01 5.1878E-08 -4.8426E-02
4.0000E+00 7.5067E-01 -8.8889E-02 -2.1015E-09 -8.3956E-01
8.0000E+00 8.7500E-01 -8.7111E+00 1.3348E-08 -9.5861E+00
1.6000E+01 9.3750E-01 -8.5689E+01 -2.5328E-14 -8.6626E+01
3.2000E+01 9.6875E-01 -7.1751E+02 -3.2076E-28 -7.1848E+02
6.4000E+01 9.8438E-01 -5.8041E+03 0.0000E+00 -5.8051E+03
1.2800E+02 9.9219E-01 -4.6561E+04 0.0000E+00 -4.6562E+04
You switched expressions at x=2.4e-4, which is when the complete
cancellation stops in the tanh() expression, but as you can see in the
last column, the truncated taylor expansion is still more accurate even
there. You have to get up to x=0.25 before the tanh() expression becomes
more accurate than the truncated taylor expansion. Also note that at
x=4, the truncated expression has the wrong sign, so there is ample
evidence of truncation errors at that point even without knowing the
exact function values for comparison.
The next term in the taylor expansion is x**5*2/945 if you want to
experiment with including that. Assuming the compiler is eliminating the
divisions for you at compile time, that just costs one more multiply-add
sequence to include it.
For a long time there was a numerical analysis newsgroup where these
things were discussed. There was always a large overlap between posters
in that group and here in c.l.f., and even many crossposts to both
groups. Fortran programmers have always been concerned with numerical
analysis. However, over time that newsgroup degenerated into SPAM for
prescription drugs and textbook solution manuals, while thankfully
c.l.f. has mostly survived. So 20 or 30 years ago, I would have posted
this in that group in order to bring in other perspectives on this topic
and to avoid the slightly off-topic subject matter.
$.02 -Ron Shepard