improve accuracy of tanh function (#2565)

improve accuracy of tanh function

using the implementation of https://www.math.utah.edu/~beebe/software/ieee/tanh.pdf

by changing from (/ (- 1 exp2z) (+ 1 exp2z)) to (- 1 (/ 2 (+ 1 exp2z)) the accuracy after rounding is increased (I was comparing with bftanh) and removes the fluctuations around z=18.35
using the polynomial for z ϵ(1.290e-8 to 0.549) seems to increase the accuracy after rounding even further
see comparison: http://pasterack.org/pastes/48436

especially the fact that  (< (tanh 18.36)(tanh 18.37)) ;=> #t was tripping me up
the two extra conditions (z . < . 1.29e-8) and (z . < . 0.549) are optional to solve this

racket/collects/racket/math.rkt

@@ -65,18 +65,29 @@
           [else                    (/ (+ (exp z) (exp (- z))) 2)]))
 
   (define (tanh z)
+    ;implementation based on https://www.math.utah.edu/~beebe/software/ieee/tanh.pdf
     (unless (number? z) (raise-argument-error 'tanh "number?" z))
     (cond [(= z 0) z]  ; preserve 0, 0.0, -0.0, 0.0f0, 0.0+0.0i, etc.
           [(real? z)
            (let loop ([z z])
              (cond [(z . < . 0)   (- (loop (- z)))]
-                   [(z . < . 20)  (define exp2z (exp (* 2 z)))
-                                  (/ (- exp2z 1) (+ exp2z 1))]
-                   [(z . >= . 20) (if (single-flonum? z) 1.0f0 1.0)]
+                   [(z . < . 1.29047841397589243466D-08) z]
+                   [(z . < . 0.54930614433405484570D+00)
+                    (define p0 -0.16134119023996228053D+04)
+                    (define p1 -0.99225929672236083313D+02)
+                    (define p2 -0.96437492777225469787D+00)
+                    (define q0 0.48402357071988688686D+04)
+                    (define q1 0.22337720718962312926D+04)
+                    (define q2 0.11274474380534949335D+03)
+                    (define g (* z z))
+                    (define R
+                      (/ (* g (+ (* (+ (* p2 g) p1) g) p0))
+                         (+ (* (+ (* (+ g q2) g) q1) g) q0)))
+                    (+ z (* z R))]
+                   [(z . < . 19.06154746539849600897D+00) (- 1 (/ 2 (+ 1 (exp (* 2 z)))))]
+                   [(z . >= . 19.06154746539849600897D+00) (if (single-flonum? z) 1.0f0 1.0)]
                    [else          z]))]  ; +nan.0 or +nan.f
-          [else
-           (define exp2z (exp (* 2 z)))
-           (/ (- exp2z 1) (+ exp2z 1))]))
+          [else (- 1 (/ 2 (+ 1 (exp (* 2 z)))))]))
 
   ;; angle conversion
   (define (degrees->radians x)