PageRenderTime 66ms CodeModel.GetById 9ms app.highlight 54ms RepoModel.GetById 1ms app.codeStats 0ms

/pkgs/math-pkgs/math-lib/math/private/flonum/flonum-error.rkt

http://github.com/plt/racket
Unknown | 190 lines | 161 code | 29 blank | 0 comment | 0 complexity | 39fd07ef94e43ac0a8ee5ed9a6ea0083 MD5 | raw file
  1#lang racket/base
  2
  3(provide flsplit
  4         fast-mono-fl+/error
  5         fast-mono-fl-/error
  6         fast-fl+/error
  7         fast-fl-/error
  8         fast-fl*/error
  9         fast-flsqr/error
 10         fast-fl//error
 11         fl+/error
 12         fl-/error
 13         fl*/error
 14         flsqr/error
 15         fl//error)
 16
 17(module untyped-defs racket/base
 18  (require (for-syntax racket/base)
 19           "flonum-functions.rkt")
 20  
 21  (provide (all-defined-out))
 22  
 23  ;(: flsplit (Flonum -> (Values Flonum Flonum)))
 24  ;; Splits a flonum into a two flonums `hi' and `lo' with 26 bits precision each, such that
 25  ;; |hi| >= |lo| and hi + lo = a. (The extra sign bit accounts for the missing bit.)
 26  ;; This function returns (values +nan.0 +nan.0) for |a| >= 1.3393857490036326e+300.
 27  (define-syntax-rule (flsplit a-expr)
 28    (let ([a a-expr])
 29      (let* ([c   (fl* a (fl+ 1.0 (flexpt 2.0 27.0)))]
 30             [x2  (fl- c (fl- c a))])
 31        (values x2 (fl- a x2)))))
 32  
 33  ;; =================================================================================================
 34  ;; Fast monotone addition and subtraction
 35  
 36  ;(: fast-mono-fl+/error (Flonum Flonum -> (Values Flonum Flonum)))
 37  ;; Returns a+b and its rounding error
 38  ;; Assumes |a| >= |b|
 39  (define-syntax-rule (fast-mono-fl+/error a-expr b-expr)
 40    (let ([a a-expr] [b b-expr])
 41      (let ([x2  (+ a b)])
 42        (values x2 (- b (- x2 a))))))
 43  
 44  ;(: fast-mono-fl-/error (Flonum Flonum -> (Values Flonum Flonum)))
 45  ;; Returns a+b and its rounding error
 46  ;; Assumes |a| >= |b|
 47  (define-syntax-rule (fast-mono-fl-/error a-expr b-expr)
 48    (let ([a a-expr] [b b-expr])
 49      (let ([x2  (- a b)])
 50        (values x2 (- (- a x2) b)))))
 51  
 52  ;; =================================================================================================
 53  ;; Fast arithmetic that returns rounding error
 54  
 55  ;(: fast-fl+/error (Flonum Flonum -> (Values Flonum Flonum)))
 56  ;; Returns a+b and its rounding error
 57  (define-syntax-rule (fast-fl+/error a-expr b-expr)
 58    (let ([a a-expr] [b b-expr])
 59      (let* ([x2  (fl+ a b)]
 60             [v   (fl- x2 a)])
 61        (values x2 (fl+ (fl- a (fl- x2 v)) (fl- b v))))))
 62  
 63  ;(: fast-fl-/error (Flonum Flonum -> (Values Flonum Flonum)))
 64  ;; Returns a-b and its rounding error
 65  (define-syntax-rule (fast-fl-/error a-expr b-expr)
 66    (let ([a a-expr] [b b-expr])
 67      (let* ([x2  (fl- a b)]
 68             [v   (fl- x2 a)])
 69        (values x2 (fl- (fl- a (fl- x2 v)) (fl+ b v))))))
 70  
 71  ;(: fast-fl*/error (Flonum Flonum -> (Values Flonum Flonum)))
 72  ;; Returns a*b and its rounding error
 73  (define-syntax-rule (fast-fl*/error a-expr b-expr)
 74    (let ([a a-expr] [b b-expr])
 75      (let*-values ([(x2)  (fl* a b)]
 76                    [(a2 a1)  (flsplit a)]
 77                    [(b2 b1)  (flsplit b)])
 78        (values x2 (- (fl- (fl- (fl- (fl- x2 (fl* a2 b2))
 79                                     (fl* a1 b2))
 80                                (fl* a2 b1))
 81                           (fl* a1 b1)))))))
 82  
 83  ;(: fast-flfma/error (Flonum Flonum Flonum -> (Values Flonum Flonum)))
 84  ;; Returns a*b+c and its rounding error
 85  (define-syntax-rule (fast-flfma/error a-expr b-expr c-expr)
 86    (let*-values ([(y2 y1)  (fast-fl*/error a-expr b-expr)]
 87                  [(h0 h1)   (fast-fl+/error c-expr y1)]
 88                  [(h3 h2)  (fast-fl+/error h0 y2)])
 89      (values h3 (fl+ h2 h1))))
 90  
 91  #;; If we had hardware fused multiply-add:
 92  (define-syntax-rule (fast-fl*/error a-expr b-expr)
 93    (let ([a a-expr] [b b-expr])
 94      (let ([x2  (fl* a b)])
 95        (values x2 (flfma a b (- x2))))))
 96
 97  ;(: fast-flsqr/error (Flonum -> (Values Flonum Flonum)))
 98  ;; Returns a*a and its rounding error
 99  (define-syntax-rule (fast-flsqr/error a-expr)
100    (let ([a a-expr])
101      (let*-values ([(x2)  (fl* a a)]
102                    [(a2 a1)  (flsplit a)])
103        (values x2 (- (fl- (fl- (fl- x2 (fl* a2 a2))
104                                (fl* 2.0 (fl* a2 a1)))
105                           (fl* a1 a1)))))))
106  
107  ;(: fast-fl//error (Flonum Flonum -> (Values Flonum Flonum)))
108  ;; Returns a/b and its rounding error
109  (define-syntax-rule (fast-fl//error a-expr b-expr)
110    (let ([a a-expr] [b b-expr])
111      (let*-values ([(x2)  (fl/ a b)]
112                    [(w2 w1)  (fast-fl*/error x2 b)])
113        (fast-mono-fl+/error x2 (fl/ (fl- (fl- a w2) w1) b)))))
114  
115  )  ; module
116
117(require (submod "." untyped-defs))
118
119(module typed-defs typed/racket/base
120  (require racket/performance-hint
121           (submod ".." untyped-defs)
122           "flonum-functions.rkt"
123           "utils.rkt")
124  
125  (provide (all-defined-out))
126  
127  ;; =================================================================================================
128  ;; Function versions of the above that are well-defined for the largest domain, and return 0.0 as
129  ;; the second argument whenever the first isn't rational
130  
131  (begin-encourage-inline
132    
133    (: fl+/error (Flonum Flonum -> (Values Flonum Flonum)))
134    (define (fl+/error a b)
135      (let-values ([(x2 x1)  (fast-fl+/error a b)])
136        (values x2 (if (flrational? x2) x1 0.0))))
137    
138    (: fl-/error (Flonum Flonum -> (Values Flonum Flonum)))
139    (define (fl-/error a b)
140      (let-values ([(x2 x1)  (fast-fl-/error a b)])
141        (values x2 (if (flrational? x2) x1 0.0))))
142    
143    (: fl*/error (Flonum Flonum -> (Values Flonum Flonum)))
144    (define (fl*/error a b)
145      (let ([x2  (fl* a b)])
146        (values x2 (if (and (flrational? x2) (not (flsubnormal? x2)))
147                       (let*-values ([(da db)  (values (near-pow2 a) (near-pow2 b))]
148                                     [(d)  (fl* da db)]
149                                     [(d?)  (and (d . fl> . 0.0) (d . fl< . +inf.0))]
150                                     [(a2 a1)  (flsplit (fl/ a da))]
151                                     [(b2 b1)  (flsplit (fl/ b db))]
152                                     [(x2)  (if d? (fl/ x2 d) (fl/ (fl/ x2 da) db))]
153                                     [(x1)  (- (fl- (fl- (fl- (fl- x2 (fl* a2 b2))
154                                                              (fl* a1 b2))
155                                                         (fl* a2 b1))
156                                                    (fl* a1 b1)))])
157                         (if d? (fl* x1 d) (fl* (fl* x1 da) db)))
158                       0.0))))
159    
160    (: flsqr/error (Flonum -> (Values Flonum Flonum)))
161    (define (flsqr/error a)
162      (let ([x2  (fl* a a)])
163        (values x2 (if (and (flrational? x2) (not (flsubnormal? x2)))
164                       (let*-values ([(d)  (near-pow2 a)]
165                                     [(d^2)  (fl* d d)]
166                                     [(d^2?)  (and (d^2 . fl> . 0.0) (d^2 . fl< . +inf.0))]
167                                     [(a2 a1)  (flsplit (fl/ a d))]
168                                     [(x2)  (if d^2? (fl/ x2 d^2) (fl/ (fl/ x2 d) d))]
169                                     [(x1)  (- (fl- (fl- (fl- x2 (fl* a2 a2))
170                                                         (fl* 2.0 (fl* a1 a2)))
171                                                    (fl* a1 a1)))])
172                         (if d^2? (fl* x1 d^2) (fl* (fl* x1 d) d)))
173                       0.0))))
174    
175    (: fl//error (Flonum Flonum -> (Values Flonum Flonum)))
176    (define (fl//error a b)
177      (let ([x2  (fl/ a b)])
178        (values x2 (if (and (flrational? x2) (flrational? b))
179                       (let* ([d  (near-pow2/div a b)]
180                              [a  (fl/ a d)]
181                              [b  (fl/ b d)])
182                         (let-values ([(w2 w1)  (fl*/error x2 b)])
183                           (fl/ (fl- (fl- a w2) w1) b)))
184                       0.0))))
185    
186    )  ; begin-encourage-inline
187  
188  )  ; module
189
190(require (submod "." typed-defs))