|
20 | 20 | [4 4096] |
21 | 21 | [5 8192])) |
22 | 22 |
|
| 23 | +(define (prev-iter) |
| 24 | + (- (*sampling-iteration*) 1)) |
| 25 | + |
23 | 26 | (define (crosses-zero? x) |
24 | 27 | (not (equal? (mpfr-sign (ival-lo x)) (mpfr-sign (ival-hi x))))) |
25 | 28 |
|
|
93 | 96 | ; k = 2: logspan(x) |
94 | 97 | (define x (first srcs)) |
95 | 98 | (define y (second srcs)) |
96 | | - (list (list (logspan y) 0) ; bounds per x |
97 | | - (list (logspan x) 0))] ; bounds per y |
| 99 | + (list (cons (logspan y) 0) ; bounds per x |
| 100 | + (cons (logspan x) 0))] ; bounds per y |
98 | 101 |
|
99 | 102 | [(ival-div) |
100 | 103 | ; k = 1: logspan(y) |
101 | 104 | ; k = 2: logspan(x) + 2 * logspan(y) |
102 | 105 | (define x (first srcs)) |
103 | 106 | (define y (second srcs)) |
104 | | - (list (list (logspan y) 0) ; bounds per x |
105 | | - (list (+ (logspan x) (* 2 (logspan y))) 0))] ; bounds per y |
| 107 | + (list (cons (logspan y) 0) ; bounds per x |
| 108 | + (cons (+ (logspan x) (* 2 (logspan y))) 0))] ; bounds per y |
106 | 109 |
|
107 | 110 | [(ival-sqrt ival-cbrt) |
108 | 111 | ; sqrt: logspan(x)/2 - 1 |
109 | 112 | ; cbrt: logspan(x)*2/3 - 1 |
110 | 113 | (define x (first srcs)) |
111 | | - (list (list (quotient (logspan x) 2) 0))] |
| 114 | + (list (cons (quotient (logspan x) 2) 0))] |
112 | 115 |
|
113 | 116 | [(ival-add ival-sub) |
114 | 117 | ; k = 1: maxlog(x) - minlog(z) |
115 | 118 | ; k = 2: maxlog(y) - minlog(z) |
116 | 119 | (define x (first srcs)) |
117 | 120 | (define y (second srcs)) |
118 | 121 |
|
119 | | - (list (list (- (maxlog x) (minlog z)) |
| 122 | + (list (cons (- (maxlog x) (minlog z)) |
120 | 123 | (- (minlog x #:no-slack #t) (maxlog z #:no-slack #t))) ; bounds per x |
121 | | - (list (- (maxlog y) (minlog z)) |
| 124 | + (cons (- (maxlog y) (minlog z)) |
122 | 125 | (- (minlog y #:no-slack #t) (maxlog z #:no-slack #t))))] ; bounds per y |
123 | 126 |
|
124 | 127 | [(ival-pow) |
|
140 | 143 | (get-slack) |
141 | 144 | 0)) |
142 | 145 |
|
143 | | - (list (list (max (+ (maxlog y) (logspan x) (logspan z) x-slack) x-slack) |
| 146 | + (list (cons (max (+ (maxlog y) (logspan x) (logspan z) x-slack) x-slack) |
144 | 147 | (minlog y #:no-slack #t)) ; bounds per x |
145 | | - (list (max (+ (maxlog y) (max (abs (maxlog x)) (abs (minlog x))) (logspan z) y-slack) |
| 148 | + (cons (max (+ (maxlog y) (max (abs (maxlog x)) (abs (minlog x))) (logspan z) y-slack) |
146 | 149 | y-slack) |
147 | 150 | (minlog y #:no-slack #t)))] ; bounds per y |
148 | 151 |
|
149 | 152 | [(ival-exp ival-exp2) |
150 | 153 | ; maxlog(x) + logspan(z) |
151 | 154 | (define x (car srcs)) |
152 | | - (list (list (+ (maxlog x) (logspan z)) (minlog x #:no-slack #t)))] |
| 155 | + (list (cons (+ (maxlog x) (logspan z)) (minlog x #:no-slack #t)))] |
153 | 156 |
|
154 | 157 | [(ival-tan) |
155 | 158 | ; maxlog(x) + max(|minlog(z)|,|maxlog(z)|) + logspan(z) + 1 |
156 | 159 | (define x (first srcs)) |
157 | | - (list (list (+ (maxlog x) (max (abs (maxlog z)) (abs (minlog z))) (logspan z) 1) |
| 160 | + (list (cons (+ (maxlog x) (max (abs (maxlog z)) (abs (minlog z))) (logspan z) 1) |
158 | 161 | (+ (minlog x #:no-slack #t) |
159 | 162 | (min (abs (maxlog z #:no-slack #t)) (abs (minlog z #:no-slack #t))))))] |
160 | 163 |
|
161 | 164 | [(ival-sin) |
162 | 165 | ; maxlog(x) - minlog(z) |
163 | 166 | (define x (first srcs)) |
164 | | - (list (list (- (maxlog x) (minlog z)) (- (maxlog z #:no-slack #t))))] |
| 167 | + (list (cons (- (maxlog x) (minlog z)) (- (maxlog z #:no-slack #t))))] |
165 | 168 |
|
166 | 169 | [(ival-cos) |
167 | 170 | ; maxlog(x) - minlog(z) + min(maxlog(x), 0) |
168 | 171 | (define x (first srcs)) |
169 | | - (list (list (+ (- (maxlog x) (minlog z)) (min (maxlog x) 0)) (- (maxlog z #:no-slack #t))))] |
| 172 | + (list (cons (+ (- (maxlog x) (minlog z)) (min (maxlog x) 0)) (- (maxlog z #:no-slack #t))))] |
170 | 173 |
|
171 | 174 | [(ival-sinh) |
172 | 175 | ; maxlog(x) + logspan(z) - min(minlog(x), 0) |
173 | 176 | (define x (first srcs)) |
174 | | - (list (list (- (+ (maxlog x) (logspan z)) (min (minlog x) 0)) (max 0 (minlog x #:no-slack #t))))] |
| 177 | + (list (cons (- (+ (maxlog x) (logspan z)) (min (minlog x) 0)) (max 0 (minlog x #:no-slack #t))))] |
175 | 178 |
|
176 | 179 | [(ival-cosh) |
177 | 180 | ; maxlog(x) + logspan(z) + min(maxlog(x), 0) |
178 | 181 | (define x (first srcs)) |
179 | | - (list (list (+ (maxlog x) (logspan z) (min (maxlog x) 0)) (max 0 (minlog x #:no-slack #t))))] |
| 182 | + (list (cons (+ (maxlog x) (logspan z) (min (maxlog x) 0)) (max 0 (minlog x #:no-slack #t))))] |
180 | 183 |
|
181 | 184 | [(ival-log ival-log2 ival-log10) |
182 | 185 | ; log: logspan(x) - minlog(z) |
183 | 186 | ; log2: logspan(x) - minlog(z) + 1 |
184 | 187 | ; log10: logspan(x) - minlog(z) - 1 |
185 | 188 | (define x (first srcs)) |
186 | | - (list (list (+ (- (logspan x) (minlog z)) 1) (- (maxlog z #:no-slack #t))))] |
| 189 | + (list (cons (+ (- (logspan x) (minlog z)) 1) (- (maxlog z #:no-slack #t))))] |
187 | 190 |
|
188 | 191 | [(ival-asin) |
189 | | - ; maxlog(x) - log[1-x^2]/2 - minlog(z) |
190 | | - ; ^^^^^^^^^^^^ |
191 | | - ; condition of uncertainty |
192 | | - (define x (first srcs)) |
193 | | - (define slack |
194 | | - (if (>= (maxlog z) 1) ; Condition of uncertainty |
195 | | - (get-slack) ; assumes that log[1-x^2]/2 is equal to slack |
196 | | - 1)) |
| 192 | + ; Γ[asin] = | x / (sqrt(1-x^2) * arcsin(x))| |
| 193 | + ; ampl[asin] = maxlog(x) - log[1-x^2]/2 - minlog(z) |
| 194 | + ; ^^^^^^^^^^^^ |
| 195 | + ; unknown part |
197 | 196 |
|
198 | | - (list (list slack 0))] |
| 197 | + (define x (first srcs)) |
| 198 | + (if (>= (maxlog z) 1) |
| 199 | + (list (cons (get-slack) |
| 200 | + (get-slack (prev-iter)))) ; assumes that log[1-x^2]/2 is equal to slack |
| 201 | + (list (cons 1 0)))] |
199 | 202 |
|
200 | 203 | [(ival-acos) |
201 | | - ; maxlog(x) - log[1-x^2]/2 - minlog(z) |
202 | | - ; ^^^^^^^^^^^^ |
203 | | - ; condition of uncertainty |
204 | | - (define x (first srcs)) |
205 | | - (define slack |
206 | | - (if (>= (maxlog x) 0) ; Condition of uncertainty |
207 | | - (get-slack) ; assumes that log[1-x^2]/2 is equal to slack |
208 | | - 0)) |
| 204 | + ; Γ[acos] = |- x / (sqrt(1-x^2) * arccos(x))| |
| 205 | + ; ampl[acos] = maxlog(x) - log[1-x^2]/2 - minlog(z) |
| 206 | + ; ^^^^^^^^^^^^ |
| 207 | + ; unknown part |
209 | 208 |
|
210 | | - (list (list slack 0))] |
| 209 | + (define x (first srcs)) |
| 210 | + (if (>= (maxlog x) 0) |
| 211 | + (list (cons (get-slack) |
| 212 | + (get-slack (prev-iter)))) ; assumes that log[1-x^2]/2 is equal to slack |
| 213 | + (list (cons 0 0)))] |
211 | 214 |
|
212 | 215 | [(ival-atan) |
213 | | - ; logspan(x) - min(|minlog(x)|, |maxlog(x)|) - minlog(z) |
| 216 | + ; Γ[atan] = | x / ((1+x^2) * arctan(x))| |
| 217 | + ; ampl[atan] = logspan(x) - min(|minlog(x)|, |maxlog(x)|) - minlog(z) |
214 | 218 | (define x (first srcs)) |
215 | | - (list (list (- (logspan x) (min (abs (minlog x)) (abs (maxlog x))) (minlog z)) 0))] |
| 219 | + (list (cons (- (logspan x) (min (abs (minlog x)) (abs (maxlog x))) (minlog z)) |
| 220 | + (- (- (max (abs (minlog x #:no-slack #t)) (abs (maxlog x #:no-slack #t)))) |
| 221 | + (maxlog z #:no-slack #t) |
| 222 | + 2)))] |
216 | 223 |
|
217 | 224 | [(ival-fmod ival-remainder) |
218 | 225 | ; x mod y = x - y*q, where q is rnd_down(x/y) |
|
228 | 235 | (get-slack) ; y crosses zero |
229 | 236 | 0)) |
230 | 237 |
|
231 | | - (list (list (- (maxlog x) (minlog z)) 0) ; bounds per x |
232 | | - (list (max (+ (- (maxlog x) (minlog z)) slack) slack) 0))] ; bounds per y |
| 238 | + (list (cons (- (maxlog x) (minlog z)) 0) ; bounds per x |
| 239 | + (cons (max (+ (- (maxlog x) (minlog z)) slack) slack) 0))] ; bounds per y |
233 | 240 |
|
234 | 241 | ; Currently log1p has a very poor approximation |
235 | 242 | [(ival-log1p) |
|
240 | 247 | (define xhi (ival-hi x)) |
241 | 248 | (define xlo (ival-lo x)) |
242 | 249 |
|
243 | | - (define slack |
244 | | - (if (or (equal? (mpfr-sign xlo) -1) (equal? (mpfr-sign xhi) -1)) |
245 | | - (get-slack) ; if x in negative |
246 | | - 0)) |
247 | | - |
248 | | - (list (list (max (+ (- (maxlog x) (minlog z)) slack) slack) 0))] |
| 250 | + (list (if (or (equal? (mpfr-sign xlo) -1) (equal? (mpfr-sign xhi) -1)) |
| 251 | + (cons (max (+ (- (maxlog x) (minlog z)) (get-slack)) (get-slack)) 0) ; if x in negative |
| 252 | + (cons (- (maxlog x) (minlog z)) 0)))] |
249 | 253 |
|
250 | 254 | ; Currently expm1 has a very poor solution for negative values |
251 | 255 | [(ival-expm1) |
252 | 256 | ; log[Гexpm1] = log[x * e^x / expm1] <= max(1 + maxlog(x), 1 + maxlog(x) - minlog(z)) |
253 | 257 | (define x (first srcs)) |
254 | | - (list (list (max (+ 1 (maxlog x)) (+ 1 (- (maxlog x) (minlog z)))) 0))] |
| 258 | + (list (cons (max (+ 1 (maxlog x)) (+ 1 (- (maxlog x) (minlog z)))) 0))] |
255 | 259 |
|
256 | 260 | [(ival-atan2) |
257 | 261 | ; maxlog(x) + maxlog(y) - 2*max(minlog(x), minlog(y)) - minlog(z) |
258 | 262 | (define x (first srcs)) |
259 | 263 | (define y (second srcs)) |
260 | 264 |
|
261 | 265 | (make-list 2 |
262 | | - (list (- (+ (maxlog x) (maxlog y)) (* 2 (max (minlog x) (minlog y))) (minlog z)) 0))] |
| 266 | + (cons (- (+ (maxlog x) (maxlog y)) (* 2 (max (minlog x) (minlog y))) (minlog z)) |
| 267 | + (- (+ (minlog x #:no-slack #t) (minlog y #:no-slack #t)) |
| 268 | + (* 2 (min (maxlog x #:no-slack #t) (maxlog y #:no-slack #t))) |
| 269 | + (maxlog z #:no-slack #t))))] |
263 | 270 |
|
264 | 271 | [(ival-tanh) |
265 | 272 | ; logspan(z) + logspan(x) |
266 | 273 | (define x (first srcs)) |
267 | | - (list (list (+ (logspan z) (logspan x)) 0))] |
| 274 | + (list (cons (+ (logspan z) (logspan x)) 0))] |
268 | 275 |
|
269 | 276 | [(ival-atanh) |
270 | 277 | ; log[Гarctanh] = maxlog(x) - log[(1-x^2)] - minlog(z) = 1 if x < 0.5, otherwise slack |
271 | 278 | ; ^^^^^^^ |
272 | 279 | ; a possible uncertainty |
273 | 280 | (define x (first srcs)) |
274 | | - (list (list (if (>= (maxlog x) 1) |
275 | | - (get-slack) |
276 | | - 1) |
277 | | - 0))] |
| 281 | + (list (if (>= (maxlog x) 1) |
| 282 | + (cons (get-slack) (get-slack (prev-iter))) |
| 283 | + (cons 1 0)))] |
278 | 284 |
|
279 | 285 | [(ival-acosh) |
280 | 286 | ; log[Гacosh] = log[x / (sqrt(x-1) * sqrt(x+1) * acosh)] <= -minlog(z) + slack |
281 | 287 | (define z-exp (minlog z)) |
282 | | - (define slack |
283 | | - (if (< z-exp 2) ; when acosh(x) < 1 |
284 | | - (get-slack) |
285 | | - 0)) |
286 | | - |
287 | | - (list (list (max (- slack z-exp) slack) 0))] |
| 288 | + (list (if (< z-exp 2) ; when acosh(x) < 1 |
| 289 | + (cons (max (- (get-slack) z-exp) (get-slack)) 0) |
| 290 | + (cons 0 0)))] |
288 | 291 |
|
289 | 292 | [(ival-pow2) |
290 | 293 | ; same as multiplication |
291 | 294 | (define x (first srcs)) |
292 | | - (list (list (+ (logspan x) 1) 0))] |
| 295 | + (list (cons (+ (logspan x) 1) 0))] |
293 | 296 |
|
294 | 297 | ; TODO |
295 | | - [(ival-erfc ival-erf ival-lgamma ival-tgamma ival-asinh ival-logb) (list (list (get-slack) 0))] |
| 298 | + [(ival-erfc ival-erf ival-lgamma ival-tgamma ival-asinh ival-logb) |
| 299 | + (list (cons (get-slack) (get-slack (prev-iter))))] |
296 | 300 | ; TODO |
297 | | - [(ival-ceil ival-floor ival-rint ival-round ival-trunc) (list (list (get-slack) 0))] |
| 301 | + [(ival-ceil ival-floor ival-rint ival-round ival-trunc) |
| 302 | + (list (cons (get-slack) (get-slack (prev-iter))))] |
298 | 303 |
|
299 | | - [else (map (λ (_) (list 0 0)) srcs)])) ; exponents for arguments |
| 304 | + [else (map (λ (_) (cons 0 0)) srcs)])) ; exponents for arguments |
0 commit comments