diff --git a/pkgs/math-pkgs/math-doc/math/scribblings/math-flonum.scrbl b/pkgs/math-pkgs/math-doc/math/scribblings/math-flonum.scrbl index dcbce61254..13ea2b3a12 100644 --- a/pkgs/math-pkgs/math-doc/math/scribblings/math-flonum.scrbl +++ b/pkgs/math-pkgs/math-doc/math/scribblings/math-flonum.scrbl @@ -286,6 +286,31 @@ identities hold: Further, choosing @racket[0.0] does not ensure that any additional identities hold. } +@defproc[(flbracketed-root [f (Flonum -> Flonum)] [a Flonum] [b Flonum]) Flonum]{ +Uses the @hyperlink["http://en.wikipedia.org/wiki/Brent%27s_method"]{Brent-Dekker method} to +find a floating-point root of @racket[f] (an @racket[x : Flonum] for which @racket[(f x)] is +very near a zero crossing) between @racket[a] and @racket[b]. +The values @racket[(f a)] and @racket[(f b)] must have opposite signs, but @racket[a] and @racket[b] +may be in any order. +@examples[#:eval untyped-eval + (define (f x) (+ 1.0 (* (+ x 3.0) (sqr (- x 1.0))))) + (define x0 (flbracketed-root f -4.0 2.0)) + (plot (list (x-axis) + (function f -4 2) + (function-label f x0)) + #:y-min -10) + (f (flprev x0)) + (f x0) + (flbracketed-root f -1.0 2.0)] +Caveats: +@(itemlist + @item{There is no guarantee that @racket[flbracketed-root] will find any @emph{particular} root. + Moreover, future updates to its implementation could make it find different ones.} + @item{There is currently no guarantee that it will find the @emph{closest} @racket[x] to an exact root.} + @item{It currently runs for at most 5000 iterations.}) +It usually requires far fewer iterations, especially if the initial bounds @racket[a] and @racket[b] are tight. +} + @defproc[(make-flexpt [x Real]) (Flonum -> Flonum)]{ Equivalent to @racket[(λ (y) (flexpt x y))] when @racket[x] is a flonum, but much more accurate for large @racket[y] when @racket[x] cannot be exactly represented diff --git a/pkgs/math-pkgs/math-lib/math/flonum.rkt b/pkgs/math-pkgs/math-lib/math/flonum.rkt index d603c0829d..11c0559ac8 100644 --- a/pkgs/math-pkgs/math-lib/math/flonum.rkt +++ b/pkgs/math-pkgs/math-lib/math/flonum.rkt @@ -12,6 +12,7 @@ "private/flonum/flonum-log1pmx.rkt" "private/flonum/flonum-polyfun.rkt" "private/flonum/flonum-error.rkt" + "private/flonum/flonum-brent-dekker.rkt" "private/flonum/expansion/expansion-base.rkt" "private/flonum/expansion/expansion-exp.rkt" "private/flonum/expansion/expansion-log.rkt" @@ -29,6 +30,7 @@ "private/flonum/flonum-log1pmx.rkt" "private/flonum/flonum-polyfun.rkt" "private/flonum/flonum-error.rkt" + "private/flonum/flonum-brent-dekker.rkt" "private/flonum/expansion/expansion-base.rkt" "private/flonum/expansion/expansion-exp.rkt" "private/flonum/expansion/expansion-log.rkt" diff --git a/pkgs/math-pkgs/math-lib/math/private/flonum/flonum-brent-dekker.rkt b/pkgs/math-pkgs/math-lib/math/private/flonum/flonum-brent-dekker.rkt new file mode 100644 index 0000000000..ec3fc463f0 --- /dev/null +++ b/pkgs/math-pkgs/math-lib/math/private/flonum/flonum-brent-dekker.rkt @@ -0,0 +1,106 @@ +#lang typed/racket/base + +(require "flonum-constants.rkt" + "flonum-functions.rkt") + +(provide flbracketed-root) + +(define debug-printf? #f) + +(define-syntax-rule (debugf e ...) + (when debug-printf? (printf e ...))) + +(: flbracketed-root* (-> (-> Flonum Flonum) Flonum Flonum Flonum Flonum Flonum)) +(define (flbracketed-root* f a fa b fb) + (let loop ([bisected? #t] [a a] [fa fa] [b b] [fb fb] [c a] [fc fa] [d 0.0] [n 0]) + (debugf "~v: ~v ~v~n" n a b) + (debugf "~v ~v~n~n" fa fb) + (define min-abs-ab (min (abs a) (abs b))) + (define ε (if (min-abs-ab . <= . +max-subnormal.0) +min.0 (* min-abs-ab epsilon.0))) + (cond + ;; If we got it right, return it + [(= fb 0.0) b] + ;; If a and b are too close, return b + [((abs (- a b)) . <= . ε) + (debugf "(abs (- a b)) <= ~v; bailing~n" ε) + b] + [(n . < . 5000) + (let*-values + ([(bisect? s) + (cond + ;; Weird rules for forcing bisection to make progress + [(or (and bisected? ((abs (- b c)) . < . (* 128.0 ε)) + (debugf "reason 1~n")) + (and (not bisected?) ((abs (- c d)) . < . (* 128.0 ε)) + (debugf "reason 2:~n"))) + (values #t 0.0)] + ;; Get an interpolated point + [else + (define fa-fb (- fa fb)) + (define fb-fc (- fb fc)) + (define fc-fa (- fc fa)) + (define da (* fa-fb (- fc-fa))) + (define db (* fb-fc (- fa-fb))) + (define dc (* fc-fa (- fb-fc))) + (cond [(or (= da 0.0) (= db 0.0) (= dc 0.0)) + ;; Secant method because quadratic method will fail + (values #f (- b (* fb (/ (- b a) (- fa-fb)))))] + [else + ;; Inverse quadratic interpolation method + (values #f (+ (/ (* a fb fc) da) + (/ (* b fc fa) db) + (/ (* c fa fb) dc)))])])] + ;; Run tests to determine whether it's a bad point + [(bisected? s) + (cond + [(or bisect? + (and (not (let ([s0 (/ (+ (* 3.0 a) b) 4.0)]) + (if (<= s0 b) (<= s0 s b) (<= b s s0)))) + (debugf "reason 3~n")) + (and bisected? ((abs (- s b)) . >= . (* 0.5 (abs (- b c)))) + (debugf "reason 4~n")) + (and (not bisected?) ((abs (- s b)) . >= . (* 0.5 (abs (- c d)))) + (debugf "reason 5~n"))) + ;; Bisect this time + (values #t (* 0.5 (+ a b)))] + [else + (values #f s)])] + [(d) c] + [(c fc) (values b fb)] + [(fs) (f s)] + ;; Replace the endpoint with the same sign as s + [(a fa b fb) (if ((* fa fs) . < . 0.0) + (values a fa s fs) + (values s fs b fb))] + ;; Make sure b is closer + [(a fa b fb) (if ((abs fa) . < . (abs fb)) + (values b fb a fa) + (values a fa b fb))]) + (loop bisected? a fa b fb c fc d (+ n 1)))] + [else b]))) + +(: flbracketed-root (-> (-> Flonum Flonum) Flonum Flonum Flonum)) +;; Find a root between a and b if f(a) and f(b) have opposite signs +(define (flbracketed-root f a b) + (define fa (f a)) + (define fb (f b)) + (cond [(= fa 0.0) a] + [(= fb 0.0) b] + [(or (flnan? a) (flnan? fa) (flnan? b) (flnan? fb)) +nan.0] + ;; Check signs + [((* fa fb) . >= . 0.0) + (debugf "(f ~v) = ~v and (f ~v) = ~v do not bracket a root~n" a fa b fb) + +nan.0] + [else + (let*-values + ([(a fa) (cond [(= a -inf.0) (values -max.0 (f -max.0))] + [(= a +inf.0) (values +max.0 (f +max.0))] + [else (values a fa)])] + [(b fb) (cond [(= b -inf.0) (values -max.0 (f -max.0))] + [(= b +inf.0) (values +max.0 (f +max.0))] + [else (values b fb)])] + ;; Make sure b is closer + [(a fa b fb) (if ((abs fa) . < . (abs fb)) + (values b fb a fa) + (values a fa b fb))]) + (flbracketed-root* f a fa b fb))]))