Added flfma/error and flfma

pkgs/math-pkgs/math-lib/math/private/flonum/flonum-error.rkt

@@ -8,11 +8,13 @@
          fast-fl*/error
          fast-flsqr/error
          fast-fl//error
+         fast-flfma/error
          fl+/error
          fl-/error
          fl*/error
          flsqr/error
-         fl//error)
+         fl//error
+         flfma/error)
 
 (module untyped-defs racket/base
   (require (for-syntax racket/base)
@@ -183,6 +185,17 @@
                            (fl/ (fl- (fl- a w2) w1) b)))
                        0.0))))
     
+    (: flfma/error (-> Flonum Flonum Flonum (Values Flonum Flonum)))
+    (define (flfma/error a b c)
+      (define-values (x2 x1) (fast-flfma/error a b c))
+      (cond [(flrational? (+ x2 x1))  (values x2 x1)]
+            [else
+             (define n (near-pow2 (max (flsqrt (abs a)) (flsqrt (abs b)))))
+             (define 1/n (/ 1.0 n))
+             (define n^2 (* n n))
+             (let-values ([(x2 x1)  (fast-flfma/error (* a 1/n) (* b 1/n) (* c 1/n 1/n))])
+               (values (* n^2 x2) (* n^2 x1)))]))
+    
     )  ; begin-encourage-inline
   
   )  ; module

pkgs/math-pkgs/math-lib/math/private/flonum/flonum-functions.rkt

@@ -38,15 +38,15 @@
   
   (define flrational?
     (λ: ([x : Flonum])
-      (and (x . fl> . -inf.0) (x . fl< . +inf.0))))
+      (fl< (flabs x) +inf.0)))
   
   (define flinfinite?
     (λ: ([x : Flonum])
-      (or (x . fl= . -inf.0) (x . fl= . +inf.0))))
+      (fl= (flabs x) +inf.0)))
   
   (define flnan?
     (λ: ([x : Flonum])
-      (not (and (x . fl>= . -inf.0) (x . fl<= . +inf.0)))))
+      (not (fl<= (flabs x) +inf.0))))
   
   (define flinteger?
     (λ: ([x : Flonum])

pkgs/math-pkgs/math-lib/math/private/flonum/flonum-more-functions.rkt

@@ -6,12 +6,14 @@
          "flonum-exp.rkt"
          "flonum-log.rkt"
          "flonum-error.rkt"
-         "flvector.rkt")
+         "flvector.rkt"
+         "utils.rkt")
 
 (provide flsqrt1pm1
          flsinh flcosh fltanh
          flasinh flacosh flatanh
-         make-flexpt flexpt+ flexpt1p)
+         make-flexpt flexpt+ flexpt1p
+         flfma)
 
 ;; ---------------------------------------------------------------------------------------------------
 ;; sqrt(1+x)-1
@@ -188,3 +190,26 @@
           [else  (flexpt (+ 1.0 x) y)]))
   
   )  ; begin-encourage-inline
+
+;; ---------------------------------------------------------------------------------------------------
+;; Fused multiply-add
+
+(: slow-flfma (-> Flonum Flonum Flonum Flonum))
+(define (slow-flfma a b c)
+  (define n (near-pow2 (max (flsqrt (abs a)) (flsqrt (abs b)))))
+  (define 1/n (/ 1.0 n))
+  (* n n (fast-flfma (* a 1/n) (* b 1/n) (* c 1/n 1/n))))
+
+(begin-encourage-inline
+  
+  (: fast-flfma (-> Flonum Flonum Flonum Flonum))
+  (define (fast-flfma a b c)
+    (let-values ([(d-hi d-lo)  (fast-flfma/error a b c)])
+      (+ d-hi d-lo)))
+
+  (: flfma (-> Flonum Flonum Flonum Flonum))
+  (define (flfma a b c)
+    (let ([d  (fast-flfma a b c)])
+      (if (flrational? d) d (slow-flfma a b c))))
+  
+  )  ; begin-encourage-inline