Reimplement quadratic-solutions to avoid cancellation and overflow

math-doc/math/scribblings/math-number-theory.scrbl

@@ -801,8 +801,17 @@ Returns a list of all real solutions to the equation @math-style{a x^2 + b x +c
                       (quadratic-solutions 1 0 -1)
                       (quadratic-solutions 1 2 1)
                       (quadratic-solutions 1 0 1)]  
+The implementation of @racket[quadratic-solutions] handles
+cancellation and overflow intelligently:
+  @interaction[#:eval untyped-eval
+                      (quadratic-solutions 1e200 2e200 1e200)
+                      (quadratic-solutions 1e-200 -2e-200 1e-200)]  
+For exact inputs, if the output can be exactly represented, it will be:
+  @interaction[#:eval untyped-eval
+                      (quadratic-solutions -1 2/3 1/3)]  
 }
 
+
 @defproc[(quadratic-integer-solutions [a Real] [b Real] [c Real]) (Listof Integer)]{
 Returns a list of all integer solutions to the equation @math-style{a x^2 + b x +c = 0}.
   @interaction[#:eval untyped-eval

math-lib/math/private/number-theory/quadratic.rkt

@@ -1,8 +1,9 @@
 #lang typed/racket/base
 
 (require racket/list
-         (only-in racket/math conjugate)
-         "types.rkt")
+         (only-in racket/math conjugate sgn nan?)
+         "types.rkt"
+         "../flonum/flonum-functions.rkt")
 
 (provide complex-quadratic-solutions  ; complex coefficients
          quadratic-solutions          ; real coefficients
@@ -37,18 +38,78 @@
               [q    (/ (+ b (* sign sqrt-d)) -2)])
          (list (/ q a) (/ c q)))])))
 
-
 (: quadratic-solutions : Real Real Real -> (Listof Real))
 (define (quadratic-solutions a b c)
-  ; return list of solutions to a a x^2 + b x + c = 0
-  (let ([d (- (* b b) (* 4 a c))])
+  (define ac-sqrt (* (flsqrt (real->double-flonum a))
+                     (flsqrt (real->double-flonum c))))
+  (define b/2 (/ b 2))
+
+  (define sqrt-d
     (cond
-      [(< d 0) '()]
-      [(= d 0) (list (/ b (* -2 a)))]
-      [else
-       (let ([sqrt-d (sqrt d)])
-         (list (/ (- (- b) sqrt-d) (* 2 a))
-               (/ (+ (- b) sqrt-d) (* 2 a))))])))
+     [(and (exact? a) (exact? b) (exact? c))
+      ; If a/b/c are exact, we want to keep as much of the computation
+      ; as possible exact, so the sqrt goes at the end
+      (define sqrt-d-number (sqrt (- (* b/2 b/2) (* a c))))
+      (if (real? sqrt-d-number)
+          sqrt-d-number
+          +nan.0)]
+     [(= (sgn a) (sgn c))
+      ; In this case we know that ac is positive so ac-sqrt has the
+      ; right sign. In this case we use difference of squares, which
+      ; allows us to do the sqrt operations without overflowing.
+      ;
+      ; So the plan is sqrt(|b/2| + ac-sqrt) sqrt(|b/2| - ac-sqrt)
+      (let* ([term1 (- (abs b/2) ac-sqrt)]
+             [term2 (+ (abs b/2) ac-sqrt)])
+             ; But there is a danger here that ac-sqrt rounded itself, either
+             ; up or down, and so the second term evaluates to zero, when it
+             ; should not. In this case, ac-sqrt has an error of ac-sqrt *
+             ; epsilon, and for the subtraction to yield zero we must have
+             ; |b/2| ~~ ac-sqrt, so the overall magnitude of the error is b
+             ; sqrt(epsilon). That's not good enough, so in this case we
+             ; expand the series out by one tick. That brings the error to a
+             ; small number of ulps, which is good enough.
+        (* (flsqrt
+            (if (= term1 0)
+                ; In this case we just need the error term of ac-sqrt. We solve:
+                ;
+                ;   (ac-sqrt - err)^2 = a c
+                ;   ac-sqrt^2 - 2 ac-sqrt err + err^2 = a c
+                ;   err = (ac-sqrt^2 - a c) / 2 ac-sqrt
+                ;       = (ac-sqrt - a c / ac-sqrt) / 2
+                ;
+                ; In this derivation I'm dropping the err^2 term
+                ; because it's too small to matter. The key is to
+                ; avoid overflow in the final formula. The most
+                ; important thing to remember here is that for c's and
+                ; ac-sqrt's exponents to have the opposite signs (and
+                ; thus for the final formula to result in overflow) we
+                ; must have a's and c's exponents to have opposite
+                ; signs, so a's and ac-sqrt's exponents must have the
+                ; same sign.
+                (if (or (and (> (abs a) 1.0) (> ac-sqrt 1.0))
+                        (and (< (abs a) 1.0) (< ac-sqrt 1.0)))
+                    (/ (- ac-sqrt (* (/ a ac-sqrt) c)) 2)
+                    (/ (- ac-sqrt (* a (/ c ac-sqrt))) 2))
+                term1))
+           ; This one has no worries about cancellation because it's
+           ; the sum of two positive values
+           (flsqrt term2)))]
+     [else
+      ; In this case hypot is perfect.
+      (flhypot (real->double-flonum b/2) ac-sqrt)]))
+
+  ; return list of solutions to a a x^2 + b x + c = 0
+  (cond
+   [(nan? sqrt-d) 
+    '()]
+   [(zero? sqrt-d)
+    (list (- (/ b/2 a)))]
+   ; Use the standard a/c swap trick to avoid cancellation
+   [(< b 0)
+    (list (/ c (+ (- b/2) sqrt-d)) (/ (+ (- b/2) sqrt-d) a))]
+   [else
+    (list (/ (- (- b/2) sqrt-d) a) (/ c (- (- b/2) sqrt-d)))]))
 
 (: quadratic-integer-solutions : Real Real Real -> (Listof Integer))
 (define (quadratic-integer-solutions a b c)

math-test/math/tests/number-theory-tests.rkt

@@ -9,6 +9,19 @@
 (check-equal? (quadratic-solutions 1 0 -4) '(-2 2))
 (check-equal? (quadratic-solutions 1 0 +4) '())
 (check-equal? (quadratic-solutions 1 0 0)  '(0))
+(check-equal? (quadratic-solutions 1 1e8 1) '(-1e8 -1e-8))
+(let ([roots (sort (quadratic-solutions 1e-8 1 1e-8) <)])
+  (check-= (first roots) -99999999.99999999 1e-7)
+  (check-= (second roots) -1.0000000000000002e-08 1e-23))
+(check-equal? (quadratic-solutions 1e200 2e200 1e200) '(-1.0))
+(check-equal? (quadratic-solutions 1e200 -2e200 1e200) '(1.0))
+(check-equal? (quadratic-solutions 1e-200 2e-200 1e-200) '(-1.0))
+(check-equal? (quadratic-solutions 1e-200 -2e-200 1e-200) '(1.0))
+; In this case, b^2 > 4 ac so there should be two solutions, but
+; `sqrt` loses the last bit, so a naive computation will yield
+; one root.
+(define sqrt-8-rounded-up (sqrt 8.0000000000000016))
+(check-equal? (length (quadratic-solutions 1 sqrt-8-rounded-up 2)) 2)
 (check-equal? (quadratic-integer-solutions 1 0 -4) '(-2 2))
 (check-equal? (quadratic-integer-solutions 1 0 +4) '())
 (check-equal? (quadratic-integer-solutions 1 0 0)  '(0))