From 53761c30f05bf9b95e9d80c7739230c31622e840 Mon Sep 17 00:00:00 2001 From: Pavel Panchekha Date: Mon, 4 Oct 2021 10:35:26 -0600 Subject: [PATCH 1/7] Fix #16 following discussion in #29 --- .../math/scribblings/math-number-theory.scrbl | 5 +++ .../math/private/number-theory/quadratic.rkt | 41 +++++++++++++++---- math-test/math/tests/number-theory-tests.rkt | 6 +++ 3 files changed, 43 insertions(+), 9 deletions(-) diff --git a/math-doc/math/scribblings/math-number-theory.scrbl b/math-doc/math/scribblings/math-number-theory.scrbl index b476568..30f7b04 100644 --- a/math-doc/math/scribblings/math-number-theory.scrbl +++ b/math-doc/math/scribblings/math-number-theory.scrbl @@ -801,6 +801,11 @@ 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 1e200 -2e200 1e200)] } @defproc[(quadratic-integer-solutions [a Real] [b Real] [c Real]) (Listof Integer)]{ diff --git a/math-lib/math/private/number-theory/quadratic.rkt b/math-lib/math/private/number-theory/quadratic.rkt index e8b5407..edaa47f 100644 --- a/math-lib/math/private/number-theory/quadratic.rkt +++ b/math-lib/math/private/number-theory/quadratic.rkt @@ -37,18 +37,41 @@ [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-values 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 + (sqrt (- (* b/2 b/2) (* a c)))] + [(= (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. + (* (flsqrt (+ (abs b/2) ac-sqrt)) (flsqrt (- (abs b/2) ac-sqrt)))] + [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 + ; Use the standard a/c swap trick to avoid cancellation + (cond + [(nan? sqrt-d) + (if (< b 0) + (list (/ (+ (- b/2) sqrt-d) a) (/ c (+ (- b/2) sqrt-d))) + (list (/ (- (- b/2) sqrt-d) a) (/ c (- (- b/2) sqrt-d))))] + [(zero? sqrt-d) + ; There is a danger here that ac-sqrt rounded itself, either up + ; or down, and these things are not actually equal. But this + ; will mostly affect the *number* of roots, not their actual + ; values, since in this case the sqrt-d + (list (- (/ b/2 a)))] + [else '()])) (: quadratic-integer-solutions : Real Real Real -> (Listof Integer)) (define (quadratic-integer-solutions a b c) diff --git a/math-test/math/tests/number-theory-tests.rkt b/math-test/math/tests/number-theory-tests.rkt index e4cdc17..a9184c9 100644 --- a/math-test/math/tests/number-theory-tests.rkt +++ b/math-test/math/tests/number-theory-tests.rkt @@ -9,6 +9,12 @@ (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)) +(check-equal? (quadratic-solutions 1e-8 1 1e-8) '(-99999999.99999999 -1.0000000000000002e-08)) +(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)) (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)) From 6ea27ba1f5fa061ee317b44927f47963eca8c5b1 Mon Sep 17 00:00:00 2001 From: Pavel Panchekha Date: Mon, 4 Oct 2021 11:56:51 -0600 Subject: [PATCH 2/7] Some bug fixes and improvements --- .../math/private/number-theory/quadratic.rkt | 58 +++++++++++++++---- math-test/math/tests/number-theory-tests.rkt | 4 ++ 2 files changed, 50 insertions(+), 12 deletions(-) diff --git a/math-lib/math/private/number-theory/quadratic.rkt b/math-lib/math/private/number-theory/quadratic.rkt index edaa47f..4aca0f7 100644 --- a/math-lib/math/private/number-theory/quadratic.rkt +++ b/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 @@ -43,33 +44,66 @@ (flsqrt (real->double-flonum c)))) (define b/2 (/ b 2)) - (define-values sqrt-d + (define sqrt-d (cond [(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 - (sqrt (- (* b/2 b/2) (* a c)))] + (flsqrt (real->double-flonum (- (* b/2 b/2) (* a c))))] [(= (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. - (* (flsqrt (+ (abs b/2) ac-sqrt)) (flsqrt (- (abs b/2) ac-sqrt)))] + ; + ; 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 = (a c - ac-sqrt^2) / 2 ac-sqrt + ; = a / 2 * c / ac-sqrt - 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))) + (/ (- (* (/ a ac-sqrt) c) ac-sqrt) 2) + (/ (- (* a (/ c ac-sqrt)) 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 - ; Use the standard a/c swap trick to avoid cancellation + ; return list of solutions to a a x^2 + b x + c = 0 (cond - [(nan? sqrt-d) + [(nan? sqrt-d) + ; Use the standard a/c swap trick to avoid cancellation (if (< b 0) (list (/ (+ (- b/2) sqrt-d) a) (/ c (+ (- b/2) sqrt-d))) (list (/ (- (- b/2) sqrt-d) a) (/ c (- (- b/2) sqrt-d))))] [(zero? sqrt-d) - ; There is a danger here that ac-sqrt rounded itself, either up - ; or down, and these things are not actually equal. But this - ; will mostly affect the *number* of roots, not their actual - ; values, since in this case the sqrt-d (list (- (/ b/2 a)))] [else '()])) diff --git a/math-test/math/tests/number-theory-tests.rkt b/math-test/math/tests/number-theory-tests.rkt index a9184c9..23e0d6b 100644 --- a/math-test/math/tests/number-theory-tests.rkt +++ b/math-test/math/tests/number-theory-tests.rkt @@ -15,6 +15,10 @@ (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. +(check-equal? (quadratic-solutions 1 (sqrt 8.0000000000000016) 2) '(1.0)) (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)) From 1ac181e0b3f15bf4fbd1dbef26cc17f81da5990c Mon Sep 17 00:00:00 2001 From: Pavel Panchekha Date: Mon, 4 Oct 2021 12:05:48 -0600 Subject: [PATCH 3/7] Woops! --- math-lib/math/private/number-theory/quadratic.rkt | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/math-lib/math/private/number-theory/quadratic.rkt b/math-lib/math/private/number-theory/quadratic.rkt index 4aca0f7..03f412a 100644 --- a/math-lib/math/private/number-theory/quadratic.rkt +++ b/math-lib/math/private/number-theory/quadratic.rkt @@ -86,8 +86,8 @@ ; same sign. (if (or (and (> (abs a) 1.0) (> ac-sqrt 1.0)) (and (< (abs a) 1.0) (< ac-sqrt 1.0))) - (/ (- (* (/ a ac-sqrt) c) ac-sqrt) 2) - (/ (- (* a (/ c ac-sqrt)) ac-sqrt) 2)) + (/ (- 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 From 42329f0923f1d5b66e7f7fcba281fecee88f7463 Mon Sep 17 00:00:00 2001 From: Pavel Panchekha Date: Mon, 4 Oct 2021 12:07:16 -0600 Subject: [PATCH 4/7] Fix test --- math-test/math/tests/number-theory-tests.rkt | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/math-test/math/tests/number-theory-tests.rkt b/math-test/math/tests/number-theory-tests.rkt index 23e0d6b..331d3b3 100644 --- a/math-test/math/tests/number-theory-tests.rkt +++ b/math-test/math/tests/number-theory-tests.rkt @@ -18,7 +18,8 @@ ; 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. -(check-equal? (quadratic-solutions 1 (sqrt 8.0000000000000016) 2) '(1.0)) +(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)) From 7e071cd25e8379b4693d6ef521f2a3d522f0fca4 Mon Sep 17 00:00:00 2001 From: Pavel Panchekha Date: Tue, 5 Oct 2021 11:48:08 -0600 Subject: [PATCH 5/7] Ensure that exact inputs produce exact outputs --- math-doc/math/scribblings/math-number-theory.scrbl | 4 ++++ math-lib/math/private/number-theory/quadratic.rkt | 5 ++++- 2 files changed, 8 insertions(+), 1 deletion(-) diff --git a/math-doc/math/scribblings/math-number-theory.scrbl b/math-doc/math/scribblings/math-number-theory.scrbl index 30f7b04..d265eb7 100644 --- a/math-doc/math/scribblings/math-number-theory.scrbl +++ b/math-doc/math/scribblings/math-number-theory.scrbl @@ -806,8 +806,12 @@ cancellation and overflow intelligently: @interaction[#:eval untyped-eval (quadratic-solutions 1e200 2e200 1e200) (quadratic-solutions 1e200 -2e200 1e200)] +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 diff --git a/math-lib/math/private/number-theory/quadratic.rkt b/math-lib/math/private/number-theory/quadratic.rkt index 03f412a..80cd400 100644 --- a/math-lib/math/private/number-theory/quadratic.rkt +++ b/math-lib/math/private/number-theory/quadratic.rkt @@ -49,7 +49,10 @@ [(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 - (flsqrt (real->double-flonum (- (* b/2 b/2) (* a c))))] + (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 From d06f3fb5f2a444dcd0dedd6341af896e717d48b3 Mon Sep 17 00:00:00 2001 From: Pavel Panchekha Date: Tue, 5 Oct 2021 11:48:31 -0600 Subject: [PATCH 6/7] Fix + update tests --- .../math/private/number-theory/quadratic.rkt | 19 ++++++++++--------- math-test/math/tests/number-theory-tests.rkt | 4 +++- 2 files changed, 13 insertions(+), 10 deletions(-) diff --git a/math-lib/math/private/number-theory/quadratic.rkt b/math-lib/math/private/number-theory/quadratic.rkt index 80cd400..21b780e 100644 --- a/math-lib/math/private/number-theory/quadratic.rkt +++ b/math-lib/math/private/number-theory/quadratic.rkt @@ -73,10 +73,10 @@ (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 = (a c - ac-sqrt^2) / 2 ac-sqrt - ; = a / 2 * c / ac-sqrt - ac-sqrt/2 + ; (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 @@ -102,13 +102,14 @@ ; return list of solutions to a a x^2 + b x + c = 0 (cond [(nan? sqrt-d) - ; Use the standard a/c swap trick to avoid cancellation - (if (< b 0) - (list (/ (+ (- b/2) sqrt-d) a) (/ c (+ (- b/2) sqrt-d))) - (list (/ (- (- b/2) sqrt-d) a) (/ c (- (- b/2) sqrt-d))))] + '()] [(zero? sqrt-d) (list (- (/ b/2 a)))] - [else '()])) + ; 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) diff --git a/math-test/math/tests/number-theory-tests.rkt b/math-test/math/tests/number-theory-tests.rkt index 331d3b3..7d57cd3 100644 --- a/math-test/math/tests/number-theory-tests.rkt +++ b/math-test/math/tests/number-theory-tests.rkt @@ -10,7 +10,9 @@ (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)) -(check-equal? (quadratic-solutions 1e-8 1 1e-8) '(-99999999.99999999 -1.0000000000000002e-08)) +(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)) From a33cbc7c62474c5e05842fa10099bc79f585cabb Mon Sep 17 00:00:00 2001 From: Pavel Panchekha Date: Tue, 5 Oct 2021 12:04:18 -0600 Subject: [PATCH 7/7] Tweak test to demonstrate underflow --- math-doc/math/scribblings/math-number-theory.scrbl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/math-doc/math/scribblings/math-number-theory.scrbl b/math-doc/math/scribblings/math-number-theory.scrbl index d265eb7..9c570cf 100644 --- a/math-doc/math/scribblings/math-number-theory.scrbl +++ b/math-doc/math/scribblings/math-number-theory.scrbl @@ -805,7 +805,7 @@ The implementation of @racket[quadratic-solutions] handles cancellation and overflow intelligently: @interaction[#:eval untyped-eval (quadratic-solutions 1e200 2e200 1e200) - (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)]