suzanne.soy/racket
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Added Brent-Dekker floating-point root finder

Browse Source
This commit is contained in:
Neil Toronto 2014-03-26 15:05:23 -06:00
parent 8069a9be54
commit 32e8ec76b0
3 changed files with 133 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

25
pkgs/math-pkgs/math-doc/math/scribblings/math-flonum.scrbl
View File

@ -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

2
pkgs/math-pkgs/math-lib/math/flonum.rkt
View File

@ -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"

106
pkgs/math-pkgs/math-lib/math/private/flonum/flonum-brent-dekker.rkt Normal file
View File

@ -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))]))