Skip to content

Commit 4274cfb

Browse files
authored
Merge pull request #115 from herbie-fp/fix-110
Document that `ival-fmod` and `ival-remainder` fail refinement
2 parents 6b4492a + d10bece commit 4274cfb

File tree

3 files changed

+135
-95
lines changed

3 files changed

+135
-95
lines changed

ops/fmod.rkt

Lines changed: 42 additions & 44 deletions
Original file line numberDiff line numberDiff line change
@@ -10,46 +10,48 @@
1010
0.bf
1111
(bfmul a b)))
1212

13+
;; WARNING: double rounding issues below, it does not refine
14+
;;
15+
;;
16+
;; (define x
17+
;; (parameterize ([bf-precision 88])
18+
;; (ival (bf "2.023002359077489336695397166e258"))))
19+
;;
20+
;; (define y
21+
;; (parameterize ([bf-precision 91])
22+
;; (ival 0.bf (bf "6.4878140443992047719110337054e233"))))
23+
;;
24+
;; (define yhi
25+
;; (parameterize ([bf-precision 91])
26+
;; (ival 0.bf (bf "6.4878140443992047719110335995e233"))))
27+
;;
28+
;;
29+
;; Here yhi refines y, so you'd think fmod(x, yhi) refines fmod(x, y)
30+
;; but in fact it doesn't! Due to double-rounding, in fmox(x, y),
31+
;; we don't think there's an intersection along the top edge, but in
32+
;; fmod(x, yhi) we falsely do, which leads to a bigger interval.
33+
;;
34+
;; For this reason ival-fmod and ival-remainder skip refinement tests.
35+
;;
36+
;; Ideally we'd fix this and get even tighter bounds, but maybe it
37+
;; doesn't matter?
38+
39+
;; Assumes both `x` and `y` are entirely positive
1340
(define (ival-fmod-pos x y err? err)
14-
;; Assumes both `x` and `y` are entirely positive
15-
(define precision (max (bf-precision) (ival-max-prec x) (ival-max-prec y)))
16-
(define a
17-
(parameterize ([bf-precision precision])
18-
(rnd 'down bftruncate (bfdiv (ival-lo-val x) (ival-hi-val y)))))
19-
(define b
20-
(parameterize ([bf-precision precision])
21-
(rnd 'up bftruncate (bfdiv (ival-hi-val x) (ival-hi-val y)))))
41+
(define a (rnd 'down bftruncate (bfdiv (ival-lo-val x) (ival-hi-val y))))
42+
(define b (rnd 'up bftruncate (bfdiv (ival-hi-val x) (ival-hi-val y))))
2243
(cond
2344
[(bf=? a b) ; No intersection along `y.hi` edge
24-
(define c
25-
(parameterize ([bf-precision precision])
26-
(rnd 'down bftruncate (bfdiv (ival-hi-val x) (ival-hi-val y)))))
27-
(define d
28-
(parameterize ([bf-precision precision])
29-
(rnd 'up bftruncate (bfdiv (ival-hi-val x) (ival-lo-val y)))))
45+
(define c (rnd 'down bftruncate (bfdiv (ival-hi-val x) (ival-hi-val y))))
46+
(define d (rnd 'up bftruncate (bfdiv (ival-hi-val x) (ival-lo-val y))))
3047
(cond
3148
[(bf=? c d) ; No intersection along `x.hi` either; use top-left/bottom-right point
32-
(define lo
33-
(rnd 'down
34-
bfsub
35-
(ival-lo-val x)
36-
(parameterize ([bf-precision precision])
37-
(rnd 'up bfmul* c (ival-hi-val y)))))
38-
(define hi
39-
(rnd 'up
40-
bfsub
41-
(ival-hi-val x)
42-
(parameterize ([bf-precision precision])
43-
(rnd 'down bfmul* c (ival-lo-val y)))))
49+
(define lo (rnd 'down bffmod (ival-lo-val x) (ival-hi-val y)))
50+
(define hi (rnd 'up bffmod (ival-hi-val x) (ival-lo-val y)))
4451
(ival (endpoint lo #f) (endpoint hi #f) err? err)]
4552
[else
4653
(ival (endpoint 0.bf #f)
47-
(endpoint (rnd 'up
48-
bfmax2
49-
(parameterize ([bf-precision precision])
50-
(bfdiv (ival-hi-val x) (rnd 'down bfadd c 1.bf)))
51-
0.bf)
52-
#f)
54+
(endpoint (rnd 'up bfdiv (ival-hi-val x) (rnd 'down bfadd c 1.bf)) #f)
5355
err?
5456
err)])]
5557
[else (ival (endpoint 0.bf #f) (endpoint (ival-hi-val y) #f) err? err)]))
@@ -63,8 +65,8 @@
6365
(or (ival-err x) (ival-err y) (and (bf=? (ival-lo-val y) 0.bf) (bf=? (ival-hi-val y) 0.bf))))
6466
(define y* (ival-exact-fabs y))
6567
(cond
66-
[(bflte? (ival-hi-val x) 0.bf) (ival-neg (ival-fmod-pos (ival-exact-neg x) y* err? err))]
67-
[(bfgte? (ival-lo-val x) 0.bf) (ival-fmod-pos x y* err? err)]
68+
[(= (mpfr-sign (ival-hi-val x)) -1) (ival-neg (ival-fmod-pos (ival-exact-neg x) y* err? err))]
69+
[(= (mpfr-sign (ival-lo-val x)) 1) (ival-fmod-pos x y* err? err)]
6870
[else
6971
(define-values (neg pos) (split-ival x 0.bf))
7072
(ival-union (ival-fmod-pos pos y* err? err)
@@ -80,14 +82,9 @@
8082
(define d (rnd 'up bfround (bfdiv (ival-hi-val x) (ival-lo-val y))))
8183
(cond
8284
[(bf=? c d) ; No intersection along `x.hi` either; use top-left/bottom-right point
83-
(define y* (rnd 'up bfdiv (ival-hi-val y) 2.bf))
84-
(ival
85-
(endpoint
86-
(rnd 'down bfmax2 (bfsub (ival-lo-val x) (rnd 'up bfmul c (ival-hi-val y))) (bfneg y*))
87-
#f)
88-
(endpoint (rnd 'up bfmin2 (bfsub (ival-hi-val x) (rnd 'down bfmul c (ival-lo-val y))) y*) #f)
89-
err?
90-
err)]
85+
(define lo (rnd 'down bfremainder (ival-lo-val x) (ival-hi-val y)))
86+
(define hi (rnd 'up bfremainder (ival-hi-val x) (ival-lo-val y)))
87+
(ival (endpoint lo #f) (endpoint hi #f) err? err)]
9188
[else
9289
;; NOPE! need to subtract half.bf one way, add it another!
9390
(define y*-hi (rnd 'up bfdiv (bfdiv (ival-hi-val x) (rnd 'down bfadd c half.bf)) 2.bf))
@@ -111,8 +108,9 @@
111108
(or (ival-err x) (ival-err y) (and (bf=? (ival-lo-val y) 0.bf) (bf=? (ival-hi-val y) 0.bf))))
112109
(define y* (ival-exact-fabs y))
113110
(cond
114-
[(bflte? (ival-hi-val x) 0.bf) (ival-neg (ival-remainder-pos (ival-exact-neg x) y* err? err))]
115-
[(bfgte? (ival-lo-val x) 0.bf) (ival-remainder-pos x y* err? err)]
111+
[(= (mpfr-sign (ival-hi-val x)) -1)
112+
(ival-neg (ival-remainder-pos (ival-exact-neg x) y* err? err))]
113+
[(= (mpfr-sign (ival-lo-val x)) 1) (ival-remainder-pos x y* err? err)]
116114
[else
117115
(define-values (neg pos) (split-ival x 0.bf))
118116
(ival-union (ival-remainder-pos pos y* err? err)

scribblings/core.scrbl

Lines changed: 67 additions & 29 deletions
Original file line numberDiff line numberDiff line change
@@ -57,24 +57,38 @@ functions, which return different intervals at different
5757

5858
@section{Interval Operations}
5959

60-
Rival provides a large set of interval operations. All of these
61-
operations are sound, meaning that the output interval always contains
62-
all valid outputs from points in the input intervals.
63-
64-
Most operations are also weakly complete, meaning that the endpoints
65-
of the output interval come from some point in the input intervals
66-
(rounded outwards). Not all operations are weakly complete, however.
60+
Rival aims to ensure three properties of all helper functions:
61+
62+
@itemlist[
63+
@item{
64+
@italic{Soundness} means output intervals contain any
65+
output on inputs drawn from the input intervals.
66+
IEEE-1788 refers to this as the output interval being @italic{valid}.
67+
}
68+
69+
@item{
70+
@italic{Refinement} means, moreover, that narrower input intervals
71+
lead to narrower output intervals. Rival's movability flags make this
72+
a somewhat more complicated property than typical.
73+
}
74+
75+
@item{
76+
@italic{Weak completeness} means, moreover, that Rival returns
77+
the narrowest possible valid interval. IEEE-1788 refers
78+
to this as the output interval being @italic{tight}.
79+
}
80+
]
6781

68-
These operations have their output precision determined by
69-
@racket[bf-precision].
82+
Weak completeness (tightness) is the strongest possible property,
83+
while soundness (validity) is the weakest, with refinement somewhere
84+
in between.
7085

7186
@deftogether[(
7287
@defproc[(ival-add [a ival?] [b ival?]) ival?]
7388
@defproc[(ival-sub [a ival?] [b ival?]) ival?]
7489
@defproc[(ival-neg [a ival?]) ival?]
7590
@defproc[(ival-mul [a ival?] [b ival?]) ival?]
7691
@defproc[(ival-div [a ival?] [b ival?]) ival?]
77-
@defproc[(ival-fma [a ival?] [b ival?] [c ival?]) ival?]
7892
@defproc[(ival-fabs [a ival?]) ival?]
7993
@defproc[(ival-sqrt [a ival?]) ival?]
8094
@defproc[(ival-cbrt [a ival?]) ival?]
@@ -86,15 +100,10 @@ These operations have their output precision determined by
86100
@defproc[(ival-log2 [a ival?]) ival?]
87101
@defproc[(ival-log10 [a ival?]) ival?]
88102
@defproc[(ival-log1p [a ival?]) ival?]
89-
@defproc[(ival-log1b [a ival?]) ival?]
90-
@defproc[(ival-pow [a ival?] [b ival?]) ival?]
91-
@defproc[(ival-sin [a ival?]) ival?]
92-
@defproc[(ival-cos [a ival?]) ival?]
93-
@defproc[(ival-tan [a ival?]) ival?]
103+
@defproc[(ival-logb [a ival?]) ival?]
94104
@defproc[(ival-asin [a ival?]) ival?]
95105
@defproc[(ival-acos [a ival?]) ival?]
96106
@defproc[(ival-atan [a ival?]) ival?]
97-
@defproc[(ival-atan2 [a ival?] [b ival?]) ival?]
98107
@defproc[(ival-sinh [a ival?]) ival?]
99108
@defproc[(ival-cosh [a ival?]) ival?]
100109
@defproc[(ival-tanh [a ival?]) ival?]
@@ -103,10 +112,6 @@ These operations have their output precision determined by
103112
@defproc[(ival-atanh [a ival?]) ival?]
104113
@defproc[(ival-erf [a ival?]) ival?]
105114
@defproc[(ival-erfc [a ival?]) ival?]
106-
@defproc[(ival-tgamma [a ival?]) ival?]
107-
@defproc[(ival-lgamma [a ival?]) ival?]
108-
@defproc[(ival-fmod [a ival?] [b ival?]) ival?]
109-
@defproc[(ival-remainder [a ival?] [b ival?]) ival?]
110115
@defproc[(ival-rint [a ival?]) ival?]
111116
@defproc[(ival-round [a ival?]) ival?]
112117
@defproc[(ival-ceil [a ival?]) ival?]
@@ -118,15 +123,48 @@ These operations have their output precision determined by
118123
@defproc[(ival-fdim [a ival?] [b ival?]) ival?]
119124
)]{
120125
These are all interval functions with arguments in the order
121-
corresponding to the same-name @code{math.h} functions. Barring
122-
bugs, all are sound. Most are weakly complete, though some more
123-
complex functions aren't, including @racket[ival-pow],
124-
@racket[ival-fma], @racket[ival-fmod], and @racket[ival-atan2]. Even
125-
these fuctions still make a best-effort attempt to produce
126-
relatively narrow intervals. For example, @racket[ival-fma] is
127-
implemented via the formula @code{(fma a b c) = (+ (* a b) c)},
128-
which that it accumulates multiple rounding errors. The result is
129-
therefore not maximally tight, but typically still pretty close.
126+
corresponding to the same-name @code{math.h} functions. The
127+
precision of the output can be set with @racket[bf-precision].
128+
All of these functions are weakly complete, returning the tightest
129+
possible intervals for the strongest possible guarantees.
130+
}
131+
132+
@deftogether[(
133+
@defproc[(ival-fma [a ival?] [b ival?] [c ival?]) ival?]
134+
@defproc[(ival-pow [a ival?] [b ival?]) ival?]
135+
@defproc[(ival-sin [a ival?]) ival?]
136+
@defproc[(ival-cos [a ival?]) ival?]
137+
@defproc[(ival-tan [a ival?]) ival?]
138+
@defproc[(ival-atan2 [a ival?] [b ival?]) ival?]
139+
)]{
140+
These interval functions, like the previous set, are analogous to
141+
the same-name @code{math.h} functions and set their precision with
142+
@racket[bf-precision]. However, these functions are more complex and
143+
do not guarantee weak completeness. We do, however, have high
144+
confidence that they satisfy the refinement property.
145+
}
146+
147+
@deftogether[(
148+
@defproc[(ival-fmod [a ival?] [b ival?]) ival?]
149+
@defproc[(ival-remainder [a ival?] [b ival?]) ival?]
150+
)]{
151+
Like the others, these interval functions take arguments and return
152+
values analogous to the same-name @code{math.h} functions and
153+
produce output with @racket[bf-precision] precision. However,
154+
these functions do not guarantee refinement in all cases due to
155+
several subtle double-rounding cases.
156+
}
157+
158+
@deftogether[(
159+
@defproc[(ival-tgamma [a ival?]) ival?]
160+
@defproc[(ival-lgamma [a ival?]) ival?]
161+
)]{
162+
These two interval functions (which take arguments and return
163+
values analogous to the same-name @code{math.h} functions and
164+
produce output with @racket[bf-precision] precision) are extremely
165+
slow, and we have only moderate confidence that these functions
166+
satisfy soundness in all cases. We do not recommended using these
167+
functions in typical use cases or at high precision.
130168

131169
@history[#:changed "1.7" @elem{Added @racket[ival-tgamma] and @racket[ival-lgamma]}]
132170
}

test.rkt

Lines changed: 26 additions & 22 deletions
Original file line numberDiff line numberDiff line change
@@ -274,6 +274,9 @@
274274
(define num-tests 1000)
275275
(define num-witnesses 10)
276276

277+
;; These fail refinement due to double-rounding sending us down the wrong branch
278+
(define non-refine-tests (list ival-fmod ival-remainder))
279+
;; These are super slow due to trying to find the minimum
277280
(define slow-tests (list ival-lgamma ival-tgamma))
278281
(define num-slow-tests 25)
279282

@@ -308,28 +311,29 @@
308311
(with-check-info (['intervals is] ['points xs] ['precs (list out-prec in-precs)])
309312
(check ival-contains? iy y)))
310313

311-
(with-check-info
312-
(['intervals is] ['points xs] ['iy iy] ['y y] ['precs (list out-prec in-precs)])
313-
(for ([k (in-naturals)]
314-
[i is]
315-
[x xs])
316-
(define-values (ilo ihi) (ival-split i x))
317-
(when (and ilo ihi)
318-
(define iylo
319-
(parameterize ([bf-precision out-prec])
320-
(apply ival-fn (list-set is k ilo))))
321-
(define iyhi
322-
(parameterize ([bf-precision out-prec])
323-
(apply ival-fn (list-set is k ihi))))
324-
(with-check-info (['split-argument k] ['ilo ilo] ['ihi ihi] ['iylo iylo] ['iyhi iyhi])
325-
(check-ival-equals? iy
326-
(parameterize ([bf-precision out-prec])
327-
(ival-union iylo iyhi))))))
328-
(when (or (ival-lo-fixed? iy) (ival-hi-fixed? iy))
329-
(define iy*
330-
(parameterize ([bf-precision 128])
331-
(apply ival-fn is)))
332-
(check ival-refines? iy iy*))))
314+
(unless (set-member? non-refine-tests ival-fn)
315+
(with-check-info
316+
(['intervals is] ['points xs] ['iy iy] ['y y] ['precs (list out-prec in-precs)])
317+
(for ([k (in-naturals)]
318+
[i is]
319+
[x xs])
320+
(define-values (ilo ihi) (ival-split i x))
321+
(when (and ilo ihi)
322+
(define iylo
323+
(parameterize ([bf-precision out-prec])
324+
(apply ival-fn (list-set is k ilo))))
325+
(define iyhi
326+
(parameterize ([bf-precision out-prec])
327+
(apply ival-fn (list-set is k ihi))))
328+
(with-check-info (['split-argument k] ['ilo ilo] ['ihi ihi] ['iylo iylo] ['iyhi iyhi])
329+
(check-ival-equals? iy
330+
(parameterize ([bf-precision out-prec])
331+
(ival-union iylo iyhi))))))
332+
(when (or (ival-lo-fixed? iy) (ival-hi-fixed? iy))
333+
(define iy*
334+
(parameterize ([bf-precision 128])
335+
(apply ival-fn is)))
336+
(check ival-refines? iy iy*)))))
333337

334338
(define (run-tests)
335339
(check ival-contains? (ival-bool #f) #f)

0 commit comments

Comments
 (0)