This is a possible new implementation for flatanh (which has a FIXME)
that uses log1p.

I computed the relative error over the interval (0,1) between the
computed result and the exact result from my computational reals package.

The exact result is not representible as a floating-point number, so
there are two ways to measure error:

1.  Measuring the relative error ("relative exact error") as

     computed-result - exact-result
     ------------------------------
         exact-result

2.  Measuring the relative error ("relative double error") as

     computed-result - correctly-rounded-exact-result
     ------------------------------------------------
          correctly-rounded-exact-result

There are accuracy results for the following routines:

system:                Ubuntu 18.04 atanh()
SRFI 144:              Current SRFI 144 routine
proposed with fllog1+: Proposed routine using SRFI 144 fllog1+
proposed with fllog1p: Proposed routine using system log1p

with the following measures of accuracy:

max relative exact error (with the associated argument)
max relative double error (with the associated argument)
RMS (root-mean-square) relative exact error
RMS (root-mean-square) relative double error

10,000 random real arguments were sampled in (0,1). All tested functions
were applied to each random sample.

Results:

system                samples= 10000 max-error-exact=
2.0566384272197089e-16 max-error-exact-arg= .12930892545177056
max-error-double= 2.2180158801390325e-16 max-error-double-arg=
.4625478839059267 RMS-error-exact= 6.048172223622568e-17
RMS-error-double= 6.872328229179074e-17
SRFI 144              samples= 10000 max-error-exact=
2.0953065632974514e-16 max-error-exact-arg= .1248274466404945
max-error-double= 2.218788527766463e-16 max-error-double-arg=
.24509421130167228 RMS-error-exact= 6.931105800579686e-17
RMS-error-double= 8.086239830220818e-17
proposed with fllog1+ samples= 10000 max-error-exact=
3.518261541328893e-16 max-error-exact-arg= .19513743524173893
max-error-double= 4.212358395187254e-16 max-error-double-arg=
.19513743524173893 RMS-error-exact= 8.719986211409874e-17
RMS-error-double= 9.477245741092292e-17
proposed with fllog1p samples= 10000 max-error-exact=
2.495093999329481e-16 max-error-exact-arg= .13142086550023874
max-error-double= 2.399275082584831e-16 max-error-double-arg=
.2273244935282261 RMS-error-exact= 6.4621088191398e-17 RMS-error-double=
7.278969295852278e-17

The proposed routine is more accurate when using the system log1p than
when using fllog1+ from SRFI 144.

Using log1p, the proposed routine has a bit smaller RMS errors than the
current SRFI 144 routine and a bit larger maximum errors.

Using fllog1+ from SRFI 144, the proposed routine has noticeably larger
maximum errors than the current SRFI 144 routine and slightly larger RMS
errors.

I don't know if the larger errors using SRFI 144's fllog1+ is a deal
breaker.

On the other hand, it does eliminate a 25-term Taylor series.  I don't
know what the trade-off of speed versus accuracy should be.

It passes all the SRFI 144 tests:

../larceny/larceny --path . --r7rs --program tests/scheme/run/flonum.sps
Trying /home/lucier/programs/srfi-144/srfi/144.constants.scm
Trying /home/lucier/programs/srfi-144/srfi/144.body0.scm
Trying /home/lucier/programs/srfi-144/srfi/144.body.scm
Trying /home/lucier/programs/srfi-144/srfi/144.special.scm
Trying /home/lucier/programs/srfi-144/srfi/144.ffi.scm
Running tests for (scheme flonum)
1276 tests passed

And using either the system log1p or the SRFI 144 fllog1+, the proposed
routine gives the right answers for all the special arguments:

 > (map atanh-fllog1+ '(-inf.0 -2.0 -1.0 -0.5 -0.0 +0.0 +0.5 +1.0 +2.0
+inf.0 +nan.0))
(+nan.0 +nan.0 -inf.0 -.5493061443340549 -0. 0. .5493061443340549 +inf.0
+nan.0 +nan.0 +nan.0)
 > (map atanh-fllog1p '(-inf.0 -2.0 -1.0 -0.5 -0.0 +0.0 +0.5 +1.0 +2.0
+inf.0 +nan.0))
(+nan.0 +nan.0 -inf.0 -.5493061443340548 -0. 0. .5493061443340548 +inf.0
+nan.0 +nan.0 +nan.0)

Just BTW, the system routine gives
 > (flatanh 0.5)

.5493061443340548

Brad