PasteRack.org
Paste # 48436
2019-03-27 08:58:53

Fork as a new paste.

Paste viewed 78 times.


Embed:

tanh test

  1. #lang racket/base
  2.  
  3. (require (prefix-in M: math)
  4.          racket/math
  5.          (only-in racket/list first last)
  6.          math/bigfloat
  7.          )
  8.  
  9. ;(bf-precision 64);auto=128
  10.  
  11. #|
  12. implementation of https://www.math.utah.edu/~beebe/software/ieee/tanh.pdf
  13. changing from (/ (- 1 exp2z) (+ 1 exp2z)) to (- 1 (/ 2 (+ 1 exp2z))
  14.  increases the accuracy (94%->96%) and removes the fluctuations aroun z=18.35
  15. using the polinomial for z ϵ(1.290e-8 0.549) increases the accuracy to ~98%
  16. see comparrison:
  17. |#
  18. #;(define (tanh z)
  19.   (unless (number? z) (raise-argument-error 'tanh "number?" z))
  20.   (cond [(= z 0) z]  ; preserve 0, 0.0, -0.0, 0.0f0, 0.0+0.0i, etc.
  21.         [(real? z)
  22.          (let loop ([z z])
  23.            (cond [(z . < . 0)   (- (loop (- z)))]
  24.                  [(z . < . 1.29047841397589243466D-08) z]
  25.                  [(z . < . 0.54930614433405484570D+00)
  26.                   (define p0 -0.16134119023996228053D+04)
  27.                   (define p1 -0.99225929672236083313D+02)
  28.                   (define p2 -0.96437492777225469787D+00)
  29.                   (define q0 0.48402357071988688686D+04)
  30.                   (define q1 0.22337720718962312926D+04)
  31.                   (define q2 0.11274474380534949335D+03)
  32.                   (define g (* z z))
  33.                   (define R
  34.                     (/ (* g (+ (* (+ (* p2 g) p1) g) p0))
  35.                        (+ (* (+ (* (+ g q2) g) q1) g) q0)))
  36.                   (+ z (* z R))]
  37.                  [(z . < . 19.06154746539849600897D+00) (- 1 (/ 2 (+ 1 (exp (* 2 x)))))]
  38.                  [(z . >= . 19.06154746539849600897D+00) (if (single-flonum? z) 1.0f0 1.0)]
  39.                  [else          z]))]  ; +nan.0 or +nan.f
  40.         [else (- 1 (/ 2 (+ 1 (exp (* 2 z)))))]))
  41.  
  42. (define (mytanh0 x)
  43.   (cond
  44.     [(< x 0)
  45.      (- (mytanh0 (- x)))]
  46.     [(< x 19.06154746539849600897D+00)
  47.      ;(define temp (- 0.5 (/ 1 (+ 1 (exp (* 2 x))))))(+ temp temp)
  48.      (- 1 (/ 2 (+ 1 (exp (* 2 x)))))
  49.      ]
  50.     [else 1.0]))
  51. (define (mytanh1 x)
  52.   (cond
  53.     [(< x 0)
  54.      (- (mytanh1 (- x)))]
  55.     [(< x 19.06154746539849600897D+00)
  56.      (define temp (- 0.5 (/ 1 (+ 1 (exp (* 2 x))))))(+ temp temp)
  57.      ;(- 1 (/ 2 (+ 1 (exp (* 2 x)))))
  58.      ]
  59.     [else 1.0]))
  60. (define (mytanh2 x)
  61.   (cond
  62.     [(< x 0)
  63.      (- (mytanh2 (- x)))]
  64.     [(< x 1.29047841397589243466D-08)
  65.      x]
  66.     [(< x 19.06154746539849600897D+00)
  67.      (define temp (- 0.5 (/ 1 (+ 1 (exp (* 2 x))))))(+ temp temp)
  68.      ;(- 1 (/ 2 (+ 1 (exp (* 2 x)))))
  69.      ]
  70.     [else 1.0]))
  71. (define (mytanh3 x)
  72.   (cond
  73.     [(< x 0)
  74.      (- (mytanh3 (- x)))]
  75.     [(< x 1.29047841397589243466D-08)
  76.      x]
  77.     [(< x 0.54930614433405484570D+00)
  78.      (define p0 -0.16134119023996228053D+04)
  79.      (define p1 -0.99225929672236083313D+02)
  80.      (define p2 -0.96437492777225469787D+00)
  81.      (define q0 0.48402357071988688686D+04)
  82.      (define q1 0.22337720718962312926D+04)
  83.      (define q2 0.11274474380534949335D+03)
  84.      (define g (* x x))
  85.      (define R
  86.        (/ (* g (+ (* (+ (* p2 g) p1) g) p0))
  87.           (+ (* (+ (* (+ g q2) g) q1) g) q0)))
  88.      (+ x (* x R))]
  89.     [(< x 19.06154746539849600897D+00)
  90.      (define temp (- 0.5 (/ 1 (+ 1 (exp (* 2 x))))))(+ temp temp)
  91.      ;(- 1 (/ 2 (+ 1 (exp (* 2 x)))))
  92.      ]
  93.     [else 1.0]))
  94.  
  95. (define (tst x)
  96.   (define bfx (bf x))
  97.  
  98.   (println
  99.    (list (list (exp x) (exp (- x)))
  100.          (list '/ (- (exp x) (exp (- x))) (+ (exp x) (exp (- x))))
  101.          (list '/ (+ (exp (* 2 x)) 1) (- (exp (* 2 x)) 1))
  102.  
  103.          (tanh x)
  104.          (M:tanh x)
  105.          (/ (- (exp x) (exp (- x))) (+ (exp x) (exp (- x))))
  106.          (- (/ (exp x) (+ (exp x) (exp (- x))))
  107.             (/ (exp (- x)) (+ (exp x) (exp (- x)))))
  108.  
  109.          ))
  110.   (println
  111.    (list (list (bfexp bfx) (bfexp (bf- bfx)))
  112.          (list '/ (bf- (bfexp bfx) (bfexp (bf- bfx))) (bf+ (bfexp bfx) (bfexp (bf- bfx))))
  113.          (list '/ (bf+ (bfexp (bf* (bf 2) bfx)) (bf 1)) (bf- (bfexp (bf* (bf 2) bfx)) (bf 1)))
  114.  
  115.          (bftanh bfx)
  116.          (bf/ (bf- (bfexp bfx) (bfexp (bf- bfx))) (bf+ (bfexp bfx) (bfexp (bf- bfx))))
  117.  
  118.          )))
  119.  
  120. (define (test f #:seed [seed 123] #:nr  [nr 10000])
  121.   (random-seed seed)
  122.   (define H (make-hash))
  123.   (for* ([i (in-range nr)]
  124.          [x (in-value (* 20 (random)))])
  125.     (hash-update! H (- (f x) (bigfloat->real(bftanh (bf x)))) add1 0))
  126.   (define avg
  127.     (for/sum ([(k v)(in-hash H)])
  128.       (* v k)))
  129.   (define dev1
  130.     (for/sum ([(k v) (in-hash H)])
  131.       (* v (expt (- k avg) 2))))
  132.   (define dev2
  133.     (for/sum ([(k v) (in-hash H)])
  134.       (* v (expt k 2))))
  135.   (define S (sort (hash->list H) < #:key car))
  136.   (list (/ (hash-ref H 0.0 0) nr 0.01) avg dev1 dev2 (first S) (last S)))
  137.  
  138. (define (multi-test lst #:seed [seed 123] #:nr [nr 50000])
  139.   (map (λ (x) (cons x (test x #:seed seed #:nr nr)))
  140.        lst))
  141.  
  142. (define LST
  143.   (list tanh mytanh0 mytanh1 mytanh2 mytanh3
  144.         (λ (x) (/ (- (exp x) (exp (- x))) (+ (exp x) (exp (- x)))))
  145.         (λ (x) (- (/ (exp x) (+ (exp x) (exp (- x))))
  146.                   (/ (exp (- x)) (+ (exp x) (exp (- x))))))))
  147. (multi-test LST)
  148. (multi-test LST #:seed 582)
  149.  
  150. (collect-garbage)(collect-garbage)
  151. (time (test tanh #:seed 385796 #:nr 500000))
  152. (collect-garbage)(collect-garbage)
  153. (time (test mytanh0 #:seed 385796 #:nr 500000))
  154. (collect-garbage)(collect-garbage)
  155. (time (test mytanh3 #:seed 385796 #:nr 500000))
  156.  
  157. (require plot)
  158.  
  159. (plot (list (function tanh 0 20 #:color 'red)
  160.             (function mytanh0 0 20 #:color 'blue)
  161.             (function mytanh3 0 20 #:color 'green)))
  162.  
  163. (plot (list (function tanh    18.25 18.55 #:color 'red)
  164.             (function mytanh3 18.25 18.55 #:color 'green)))

=>