From 14bbd662e9c6b13aaf348d2878147be5b420fdb0 Mon Sep 17 00:00:00 2001
From: Neil Toronto <neil.toronto@gmail.com>
Date: Tue, 22 Apr 2014 10:46:19 -0600
Subject: [PATCH] Added flfma/error and flfma

---
 .../math/private/flonum/flonum-error.rkt      | 17 +++++++++--
 .../math/private/flonum/flonum-functions.rkt  |  6 ++--
 .../private/flonum/flonum-more-functions.rkt  | 29 +++++++++++++++++--
 3 files changed, 45 insertions(+), 7 deletions(-)

diff --git a/pkgs/math-pkgs/math-lib/math/private/flonum/flonum-error.rkt b/pkgs/math-pkgs/math-lib/math/private/flonum/flonum-error.rkt
index 7933074e11..699cf557ca 100644
--- a/pkgs/math-pkgs/math-lib/math/private/flonum/flonum-error.rkt
+++ b/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)
@@ -84,7 +86,7 @@
   ;; Returns a*b+c and its rounding error
   (define-syntax-rule (fast-flfma/error a-expr b-expr c-expr)
     (let*-values ([(y2 y1)  (fast-fl*/error a-expr b-expr)]
-                  [(h0 h1)   (fast-fl+/error c-expr y1)]
+                  [(h0 h1)  (fast-fl+/error c-expr y1)]
                   [(h3 h2)  (fast-fl+/error h0 y2)])
       (values h3 (fl+ h2 h1))))
   
@@ -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
diff --git a/pkgs/math-pkgs/math-lib/math/private/flonum/flonum-functions.rkt b/pkgs/math-pkgs/math-lib/math/private/flonum/flonum-functions.rkt
index 73357d663a..3d179535a4 100644
--- a/pkgs/math-pkgs/math-lib/math/private/flonum/flonum-functions.rkt
+++ b/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])
diff --git a/pkgs/math-pkgs/math-lib/math/private/flonum/flonum-more-functions.rkt b/pkgs/math-pkgs/math-lib/math/private/flonum/flonum-more-functions.rkt
index 5cc3fdf91b..88a42f1dd5 100644
--- a/pkgs/math-pkgs/math-lib/math/private/flonum/flonum-more-functions.rkt
+++ b/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