#### /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
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))
```