Browse files

Use ULPs to measure error, not relative error.

  • Loading branch information...
pavpanchekha committed Mar 20, 2014
1 parent e08db35 commit d28f281ee0ea852c7cf591344313c228e8f226a8
Showing with 25 additions and 19 deletions.
  1. +5 −5 casio/analyze-local-error.rkt
  2. +10 −2 casio/common.rkt
  3. +10 −12 casio/points.rkt
@@ -1,6 +1,6 @@
#lang racket
(require racket/match)
(require math/flonum)
(require math/bigfloat)
(require casio/common)
(require casio/points)
@@ -36,14 +36,14 @@
(λ (c loc)
(let* ([exact (repeat (bf c))] [approx (repeat ((*precision*) c))]
[error (repeat (relative-error (->flonum exact) approx))])
[error (repeat (1+ (flulp-error (->flonum c) (->flonum c))))])
(annotation c exact approx error error loc)))
(λ (v loc)
(let* ([var (for/list ([vmap varmap]) (cdr (assoc v vmap)))]
[exact (map bf var)] [approx (map (*precision*) var)]
[error (map relative-error (map ->flonum exact) approx)])
[error (map (compose 1+ flulp-error) (map ->flonum exact) (map ->flonum approx))])
(annotation v exact approx error error loc)))
@@ -61,9 +61,9 @@
(map (curry apply approx-op) approx-inputs)]
(map relative-error (map ->flonum exact-ans) semiapprox-ans)]
(map (compose flulp-error) (map ->flonum exact-ans) (map ->flonum semiapprox-ans))]
(map relative-error (map ->flonum exact-ans) approx-ans)])
(map (compose 1+ flulp-error) (map ->flonum exact-ans) (map ->flonum approx-ans))])
(annotation expr exact-ans approx-ans local-error cumulative-error loc)))))
(define (interesting-error? l)
@@ -4,7 +4,7 @@
(require data/order)
(provide reap println ->flonum *precision* cotan square ordinary-float?
list= list< enumerate take-up-to *debug* debug debug-reset pipe)
list= list< enumerate take-up-to *debug* debug debug-reset pipe 1+)
; Precision for approximate evaluation
(define *precision* (make-parameter real->double-flonum))
@@ -54,7 +54,12 @@
(define (->flonum x)
[(real? x) ((*precision*) x)]
[(bigfloat? x) ((*precision*) (bigfloat->flonum x))]))
[(bigfloat? x) ((*precision*) (bigfloat->flonum x))]
[(complex? x)
(if (= (imag-part x) 0)
(->flonum (real-part x))
[else (error "Invalid number" x)]))
; Functions used by our benchmarks
(define (cotan x)
@@ -70,6 +75,9 @@
(or (= x1 x2)
(and (nan? x1) (nan? x2))))
(define (1+ x)
(+ 1 x))
(define (list= l1 l2)
(and l1 l2 (andmap =-or-nan? l1 l2)))
@@ -1,12 +1,12 @@
#lang racket
(require math/flonum)
(require math/bigfloat)
(require casio/common)
(require casio/programs)
(provide prepare-points *points* *exacts*
errors errors-compare errors-difference
relative-error make-exacts)
errors errors-compare errors-difference make-exacts)
(define *points* (make-parameter '()))
(define *exacts* (make-parameter '()))
@@ -81,21 +81,19 @@
[exacts* (filter-exacts pts exacts)])
(values pts* exacts*)))
(define (relative-error approx exact)
(if (and (real? approx) (real? exact))
(abs (/ (- exact approx) (max (abs exact) (abs approx))))
(define (errors prog points exacts)
(let ([fn (eval-prog prog mode:fl)])
(for/list ([point points] [exact exacts])
(relative-error (fn point) exact))))
(let ([out (fn point)])
(if (real? out)
(+ 1 (flulp-error out (->flonum exact)))
(define errors-compare-cache (make-hasheq))
(define (reasonable-error? x)
; TODO : Why do we need the 100% error case?
(not (or (= x 1.0) (infinite? x) (nan? x))))
(not (or (infinite? x) (nan? x))))
(define (errors-compare errors1 errors2)
(map (λ (x) (cond [(< x 0) '<] [(> x 0) '>] [#t '=]))
@@ -115,9 +113,9 @@
; but that lead to stupid behavior;
; a least-significant-bit error versus no error
; would be very heavily weighed
[(= error1 0) (log (/ 1e-16 error2))]
[(= error2 0) (log (/ error1 1e-16))]
[#t (log (/ error1 error2))])]
[(= error1 0) (log error2)]
[(= error2 0) (log error1)]
[#t (/ (log (/ error1 error2)) (log 2))])]
[(or (and (reasonable-error? error1) (not (reasonable-error? error2))))
[(or (and (not (reasonable-error? error1)) (reasonable-error? error2)))

0 comments on commit d28f281

Please sign in to comment.