fllog+, fllog-quotient & fllog-hypot

document (and improve) fllog+ and fllog-quotient, and add fllog-hypot

fllog+
* change boundaries for double-double calculation
`(fflog+ -1.4293498385921057e-012 +2.5111776760269606)`

fllog-quotient
* correct handling if x is +nan.0
* double-double calculation if x/y is near 1
* avoid unnecessary multiplication with s
`(fflog-quotient 500000000001. +500000000000.)`
`(fflog-quotient +nan.0 +2.)`

fllog-hypot
* avoid over/underflow
* accuracy where (hypot x y) is near 1
`(fllog-hypot 1e-321 +1e-320)`
`(fllog-hypot 1.1084129893159664 +1.619387121491612e-005)`

Author: bdeket

math-doc/math/scribblings/math-flonum.scrbl

@@ -121,6 +121,13 @@ For factorial-like functions that return sensible values for non-integers, see
 @racket[gamma] and @racket[beta].
 }
 
+@deftogether[(@defproc[(fllog+ [a Flonum] [b Flonum]) Flonum]
+              @defproc[(fllog-quotient [a Flonum] [b Flonum]) Flonum]
+              @defproc[(fllog-hypot [a Flonum] [b Flonum]) Flonum])]{
+Like @racket[(fllog (+ a b))] , @racket[(fllog (/ a b))] and @racket[(fllog (flhypot a b))]
+but avoiding over/underflow and accurate when the argument to @racket[fllog] is near 1
+(within 1 @tech{ulp}).}
+
 @deftogether[(@defproc[(fllog-factorial [n Flonum]) Flonum]
               @defproc[(fllog-binomial [n Flonum] [k Flonum]) Flonum]
               @defproc[(fllog-permutations [n Flonum] [k Flonum]) Flonum]

math-lib/math/private/flonum/flonum-log.rkt

@@ -14,7 +14,8 @@
          lg1+ lg+ lg1- lg- lgsum
          fllog-quotient
          fllog2
-         fllogb)
+         fllogb
+         fllog-hypot)
 
 (begin-encourage-inline
 
@@ -99,9 +100,10 @@
 
   (: fllog+ (Flonum Flonum -> Flonum))
   ;; Computes log(a+b) in a way that is accurate for a+b near 1.0
+  ;; random testing shows 99.9% <0.5ulp and 100%< 1ulp
   (define (fllog+ a b)
     (define a+b (+ a b))
-    (cond [((flabs (- a+b 1.0)) . < . (fllog 2.0))
+    (cond [(< 0.3 a+b 2.7)
            ;; a+b is too close to 1.0, so compute in higher precision
            (define-values (a+b a+b-lo) (fast-fl+/error a b))
            (- (fllog a+b) (fllog1p (- (/ a+b-lo a+b))))]
@@ -167,16 +169,21 @@
                           (+ max-log-x (fllog1p s)))))))))))
 
 (: fllog-quotient (Flonum Flonum -> Flonum))
-;; Computes (fllog (/ x y)) in a way that reduces error and avoids under-/overflow
+;; Computes (fllog (/ x y)) in a way that reduces error, avoids under-/overflow
+;; and is more accurate for x ≈ y.  random testing shows 96% <0.5 ulp 99.5%<1ulp
+;; and 100% < 1.3 (when x very big and y very small or reversed)
 (define (fllog-quotient x y)
   (let ([x  (flabs x)]
         [y  (flabs y)]
         [s  (fl/ (flsgn x) (flsgn y))])
     (cond [(s . fl> . 0.0)
            (define z (fl/ x y))
-           (cond [(and (z . fl> . +max-subnormal.0) (z . fl< . +inf.0))  (fllog (fl* s z))]
+           (cond [(< 0.5 z 2.5)
+                  (define-values (z z_) (fl//error x y))
+                  (fl- (fllog z) (fllog1p (- (/ z_ z))))]
+                 [(and (z . fl> . +max-subnormal.0) (z . fl< . +inf.0))  (fllog z)]
                  [else  (fl+ (fllog x) (- (fllog y)))])]
-          [(s . fl= . 0.0)  -inf.0]
+          [(and (s . fl= . 0.0) (not (flnan? x)))  -inf.0]
           [else  +nan.0])))
 
 (define log-max.0 (fllog +max.0))
@@ -264,3 +271,24 @@
                 ;; Oh noes! We had overflows or underflows!
                 ;; Not a lot we can do without introducing more error, so just return y
                 y])]))
+
+(: fllog-hypot (-> Flonum Flonum Flonum))
+;; computes (log (hypot x y)) avoiding over/underflow and more accurate when (hypot x y) ≈ 1.0
+;; random testing shows 99.9% < 0.5ulp and 100% < 1ulp
+(define (fllog-hypot X Y)
+  (define x (flabs X))
+  (define y (flabs Y))
+  (define m (flhypot x y))
+  (cond
+    [(or (<= m .23e-307)
+         (and (flinfinite? m) (flrational? x) (flrational? y)))
+     (fl* 0.5 (lg+ (fl* 2. (fllog x)) (fl* 2. (fllog y))))]
+    [(<= 0.15 m 2.7)
+     (define-values (x² x²_) (fl*/error x x))
+     (define-values (y² y²_) (fl*/error y y))
+     (define-values (x+y x+y_) (fl2+ x² x²_ y² y²_))
+     (fl* 0.5 (fl- (fllog x+y) (fllog1p (- (/ x+y_ x+y)))))]
+    [(<= 0.025 m 25.)
+     (fl* 0.5 (fllog (+ (flexpt x 2.) (flexpt y 2.))))]
+    [else
+     (fllog m)]))

math-test/math/tests/function-tests.rkt

@@ -101,3 +101,13 @@
 (check-= (Fresnel-C 20) 0.4999873349723443881870062136976602164476 (ε))
 (check-= (Fresnel-C 50) 0.4999991894307279679558101639817919070024 (ε))
 
+;; ---------------------------------------------------------------------------------------------------
+;; fllog+ fllog-quotient fllog-hypot
+(check-= (fllog+ 1. 1e-20) 1e-20 (* 10 (flulp 1e-20)))
+(check-= (fllog-quotient +min.0 2.0) -745.1332191019411 1.e-10)
+(check-= (fllog-quotient 1.8244133740471687e+267 +8.64804592730887e-063) 758.2970057934524 1.)
+(check-= (fllog-quotient 1.8244133740471687e+267 +8.64804592730887e-063) 758.2970057934524 1.)
+(check-= (fllog-quotient 500000000001. +500000000000.) 2e-12 1e-20)
+(check-= (fllog-hypot #i3/5 #i4/5) 2.22e-17 1e-18)
+(check-= (fllog-hypot +min.0 +min.0) -744.0934983311013 1e-10)
+(check-= (fllog-hypot +max.0 +max.0) 710.1292864836639 1e-10)