(module math.flonum.more-functions (fpsqrt1pm1 fpsinh fpcosh fptanh fpasinh fpacosh fpatanh make-fpexpt fpexpt+ fpexpt1p fpfma) (import scheme chicken.type chicken.flonum (only chicken.base define-values let-values assert include) math.flonum.functions math.flonum.constants math.flonum.exp math.flonum.log math.flonum.error math.flonum.fpvector math.flonum.utils) (include "math-types.scm") ;; --------------------------------------------------------------------------------------------------- ;; sqrt(1+x)-1 (: fpsqrt1pm1 (float -> float)) (define (fpsqrt1pm1 x) (cond [(fp> (fpabs x) 0.75) (fp- (fpsqrt (fp+ 1.0 x)) 1.0)] [else (fpexpm1 (fp* 0.5 (fplog1p x)))])) ;; --------------------------------------------------------------------------------------------------- ;; Hyperbolic sine (: fpsinh (float -> float)) (define (fpsinh x) (cond [(fp< x 0.0) ;; Odd function (- (fpsinh (- x)))] [(fp< x (fpexpt 2.0 -26.0)) ;; sinh(x) ~ x x] [(fp< x 18.5) ;; sinh(x) = (exp(2*x) - 1) / (2*exp(x)) (let ((y (fpexpm1 x))) (fp* 0.5 (fp+ y (fp/ y (fp+ y 1.0)))))] [(fp< x (fplog +max.0)) ;; sinh(x) ~ exp(x) / 2 (fp* 0.5 (fpexp x))] [else ;; sinh(x) ~ exp(x) / 2 = (exp(x/2) / 2) * exp(x/2) (let ((y (fpexp (fp* 0.5 x)))) (fp* (fp* 0.5 y) y))])) ;; --------------------------------------------------------------------------------------------------- ;; Hyperbolic cosine (: fpcosh (float -> float)) (define (fpcosh x) ;; cosh(x) = cosh(-x) (let ([x (fpabs x)]) (cond [(fp< x (fpexpt 2.0 -26.0)) ;; cosh(x) ~ 1 1.0] [(fp< x (fp* 0.5 (fplog 2.0))) ;; cosh(x) = 1 + (exp(x) - 1)^2 / (2*exp(x)) (let ((y (fpexpm1 x))) (fp+ 1.0 (fp/ (fp* y y) (fp* 2.0 (fp+ 1.0 y)))))] [(fp< x 18.5) ;; cosh(x) = (exp(x) + 1/exp(x)) / 2 (let ((y (fpexp x))) (fp+ (fp* 0.5 y) (fp/ 0.5 y)))] [(fp< x (fplog +max.0)) ;; cosh(x) ~ exp(x) / 2 (fp* 0.5 (fpexp x))] [else ;; cosh(x) ~ exp(x) / 2 = (exp(x/2) / 2) * exp(x/2) (let ((y (fpexp (fp* 0.5 x)))) (fp* (fp* 0.5 y) y))]))) ;; --------------------------------------------------------------------------------------------------- ;; Hyperbolic tangent (: fptanh (float -> float)) (define (fptanh x) (cond [(fp< x 0.0) ;; tanh(x) = -tanh(-x) (- (fptanh (- x)))] [(fp< x 1e-16) ;; tanh(x) ~ x + x^2 (fp* x (fp+ 1.0 x))] [(fp< x 0.5) ;; tanh(x) = (exp(2*x) - 1) / (exp(2*x) + 1) (let ((y (fpexpm1 (fp* -2.0 x)))) (- (fp/ y (fp+ 2.0 y))))] [(fp< x 19.5) ;; tanh(x) = (exp(2*x) - 1) / (exp(2*x) + 1) (let ((y (fpexp (fp* 2.0 x)))) (fp/ (fp- y 1.0) (fp+ y 1.0)))] [(fp<= x +inf.0) ;; tanh(x) ~ 1 1.0] [else +nan.0])) ;; --------------------------------------------------------------------------------------------------- ;; Inverse hyperbolic sine (: fpasinh (float -> float)) (define (fpasinh x) (cond [(fp< x 0.0) (- (fpasinh (- x)))] [(fp< x 2e-8) x] [(fp< x 0.00018) ;; Taylor series order 3 (fp* x (fp+ 1.0 (fp* (fp* #i-1/6 x) x)))] [(fp< x 1.0) ;; Standard definition, rearranged to preserve digits (fplog1p (fp+ x (fpsqrt1pm1 (fp* x x))))] [(fp< x 3e3) ;; Standard definition (fplog (fp+ x (fpsqrt (fp+ (fp* x x) 1.0))))] [(fp< x 1e307) ;; Laurent series in 1/x at 0+ order from -1 to 1 (fp+ (fplog (fp* x 2.0)) (fp/ 1.0 (fp* (fp* 4.0 x) x)))] [else ;; Laurent series, rearranged to not overfpow (fp+ (fplog x) (fplog 2.0))])) ;; --------------------------------------------------------------------------------------------------- ;; Inverse hyperbolic cosine (: fpacosh (float -> float)) (define (fpacosh x) (cond [(fp< x 1.0) +nan.0] [(fp< x 1.5) ;; Standard definition, rearranged to preserve digits when x is near 1.0 (let ((y (fp- x 1.0))) (fplog1p (fp+ y (fpsqrt (fp+ (fp* y y) (fp* 2.0 y))))))] [(fp< x 1e8) ;; Standard definition (fplog (fp+ x (fpsqrt (fp- (fp* x x) 1.0))))] [(fp< x 1e307) ;; Laurent series in 1/x at 0+ order from -1 to 0 (fplog (fp* x 2.0))] [else ;; Laurent series, rearranged to avoid overfpow (fp+ (fplog x) (fplog 2.0))])) ;; --------------------------------------------------------------------------------------------------- ;; Inverse hyperbolic tangent (: flatanh (float -> float)) (define (fpatanh x) (cond [(fp< x 0.0) (- (fpatanh (- x)))] [(fp< x 1e-8) x] [(fp< x 0.00015) ;; Taylor series order 2 (fp+ x (fp* (fp* (fp* #i1/3 x) x) x))] [(fp< x 0.5) ;; Standard definition, rearranged to preserve digits when x is near 0.0 (fp* 0.5 (fp- (fplog1p x) (fplog1p (- x))))] [(fp< x 1.0) ;; Standard definition (fp* 0.5 (fplog (fp/ (fp+ 1.0 x) (fp- 1.0 x))))] [(fp= x 1.0) +inf.0] [else +nan.0])) ;; --------------------------------------------------------------------------------------------------- ;; Exponential with high-precision bases (: make-fpexpt (exact-rational -> (float -> float))) (define (make-fpexpt b) (assert (positive? b) 'make-fpexpt "not a positive ratnum:" b) (define b-hi (exact->inexact b)) (define b-lo (exact->inexact (- (/ (inexact->exact b-hi) b) 1))) (cond [(fp= b-lo 0.0) (lambda (x) (fpexpt b-hi x))] [else (lambda (x) (fp/ (fpexpt b-hi x) (fpexp (fp* x (fplog1p b-lo)))))])) (: fpexpt+ (float float float -> float)) (define (fpexpt+ a b y) (define-values (x-hi x-lo) (fast-fp+/error a b)) (fp/ (fpexpt x-hi y) (fpexp (fp* y (fplog1p (- (/ x-lo x-hi))))))) (: fpexpt1p (float float -> float)) (define (fpexpt1p x y) (cond [(and (> x -0.5) (< x +inf.0)) (let-values ([(a-hi a-lo) (fast-fp+/error 1.0 x)]) (fp/ (fpexpt a-hi y) (fpexp (fp* y (fplog1p (- (/ a-lo a-hi)))))))] [else (fpexpt (+ 1.0 x) y)])) ;; --------------------------------------------------------------------------------------------------- ;; Fused multiply-add (: slow-fpfma (float float float -> float)) (define (slow-fpfma a b c) (define n (near-pow2 (max (fpsqrt (abs a)) (fpsqrt (abs b))))) (define 1/n (/ 1.0 n)) (* n n (fast-fpfma (* a 1/n) (* b 1/n) (* c 1/n 1/n)))) (: fast-flfma (float float float -> float)) (define (fast-fpfma a b c) (let-values ([(d-hi d-lo) (fast-fpfma/error a b c)]) (+ d-hi d-lo))) (: flfma (float float float -> float)) (define (fpfma a b c) (let ([d (fast-fpfma a b c)]) (if (fprational? d) d (slow-fpfma a b c)))))