diff --git a/math-doc/math/scribblings/math-number-theory.scrbl b/math-doc/math/scribblings/math-number-theory.scrbl index b476568..9c570cf 100644 --- a/math-doc/math/scribblings/math-number-theory.scrbl +++ b/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 diff --git a/math-lib/math/private/number-theory/quadratic.rkt b/math-lib/math/private/number-theory/quadratic.rkt index e8b5407..21b780e 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 @@ -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) diff --git a/math-test/math/tests/number-theory-tests.rkt b/math-test/math/tests/number-theory-tests.rkt index e4cdc17..7d57cd3 100644 --- a/math-test/math/tests/number-theory-tests.rkt +++ b/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))