Skip to content

Direct rounding for ival-pow and ival-pow2 #138

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 2 commits into from
Oct 17, 2025
Merged

Direct rounding for ival-pow and ival-pow2 #138

Show file tree
Hide file tree
Changes from 1 commit
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
2557bd2
Implements direct rounding for pow and pow2
obround Oct 17, 2025
c4c1986
Use mpfr_sqr
obround Oct 17, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 4 additions & 0 deletions mpfr.rkt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -142,6 +142,9 @@

(define mpfr-set-prec! set-mpfr-prec!)

(define (mpfr-pow2! out a rnd)
(mpfr-mul! out a a rnd))

(define (bfremainder x mod)
(define out (bf 0))
(mpfr-remainder! out x mod (bf-rounding-mode))
Expand Down Expand Up @@ -346,5 +349,6 @@
mpfr-log!
mpfr-const-pi!
mpfr-atan2!
mpfr-pow2!
mpfr-set-prec!
mpfr-set!)
6 changes: 6 additions & 0 deletions ops/core.rkt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -32,6 +32,12 @@
ival-union
ival-split
(rename-out [monotonic monotonic->ival] [comonotonic comonotonic->ival])
define*
immutable
monotonic-mpfr!
comonotonic-mpfr!
monotonic-mpfr
comonotonic-mpfr
ival-illegal
ival-pi
ival-e
Expand Down
68 changes: 34 additions & 34 deletions ops/pow.rkt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -13,15 +13,15 @@
[(and (< (mpfr-exp (ival-hi-val x)) 1) (not (bfinfinite? (ival-hi-val x)))) -1]
[else 0]))

(define (eppow a a! b b! a-class b-class)
(when (bfzero? a)
(set! a 0.bf)) ; Handle (-0)^(-1)
(define-values (val exact?) (bf-return-exact? bfexpt (list a b)))
(values val
(define (eppow! out a a! b b! a-class b-class rnd)
(define base (if (bfzero? a) 0.bf a)) ; Handle (-0)^(-1)
(mpfr-set-prec! out (bf-precision))
(define exact? (= 0 (mpfr-pow! out base b rnd)))
(values out
(or (and a! b! exact?)
(and a! (bf=? a 1.bf))
(and a! (bfzero? a) (not (= b-class 0)))
(and a! (bfinfinite? a) (not (= b-class 0)))
(and a! (bf=? base 1.bf))
(and a! (bfzero? base) (not (= b-class 0)))
(and a! (bfinfinite? base) (not (= b-class 0)))
(and b! (bfzero? b))
(and b! (bfinfinite? b) (not (= a-class 0))))))

Expand All @@ -39,24 +39,25 @@
(define y-class (classify-ival y))

(define (mk-pow a b c d)
(define out (new-ival (bf-precision)))
(define-values (lo lo!)
(rnd 'down
eppow
(endpoint-val a)
(endpoint-immovable? a)
(endpoint-val b)
(endpoint-immovable? b)
x-class
y-class))
(eppow! (ival-lo-val out)
(endpoint-val a)
(endpoint-immovable? a)
(endpoint-val b)
(endpoint-immovable? b)
x-class
y-class
'down))
(define-values (hi hi!)
(rnd 'up
eppow
(endpoint-val c)
(endpoint-immovable? c)
(endpoint-val d)
(endpoint-immovable? d)
x-class
y-class))
(eppow! (ival-hi-val out)
(endpoint-val c)
(endpoint-immovable? c)
(endpoint-val d)
(endpoint-immovable? d)
x-class
y-class
'up))

(define-values (real-lo! real-hi!)
(cond
Expand All @@ -69,17 +70,16 @@
;; Can't use >=, even though exp2-overflow-threshold is a
;; power of 2, because mpfr-exp is offset by 1 from the real
;; exponent, which matters when we add them.
(define log2-base (bf 0))
(define (log2-sum-exceeds-threshold? exponent-value base-value)
(mpfr-log2! log2-base base-value 'zero)
(> (+ (mpfr-exp exponent-value) (mpfr-exp log2-base)) (mpfr-exp exp2-overflow-threshold)))

(define must-overflow
(and (bfinfinite? hi)
(= (* x-class y-class) 1)
(> (+ (mpfr-exp bval) (mpfr-exp (rnd 'zero bflog2 aval)))
(mpfr-exp exp2-overflow-threshold))))
(and (bfinfinite? hi) (= (* x-class y-class) 1) (log2-sum-exceeds-threshold? bval aval)))

(define must-underflow
(and (bfzero? lo)
(= (* x-class y-class) -1)
(> (+ (mpfr-exp dval) (mpfr-exp (rnd 'zero bflog2 cval)))
(mpfr-exp exp2-overflow-threshold))))
(and (bfzero? lo) (= (* x-class y-class) -1) (log2-sum-exceeds-threshold? dval cval)))

(define real-lo! (or lo! must-underflow (and (bfzero? lo) a! b!)))
(define real-hi! (or hi! must-underflow must-overflow (and (bfinfinite? hi) c! d!)))
Expand Down Expand Up @@ -152,8 +152,8 @@
(let ([pospow (ival-pow-pos (ival-exact-fabs x) y)])
(ival-then (ival-assert ival-maybe) (ival-union (ival-neg pospow) pospow)))))

(define (ival-pow2 x)
((monotonic->ival (lambda (x) (bfmul x x))) (ival-pre-fabs x)))
(define* ival-pow2! (lambda (out x) ((monotonic-mpfr! mpfr-pow2!) out (ival-pre-fabs x))))
(define* ival-pow2 (immutable ival-pow2!))

(define (ival-pow x y)
(cond
Expand Down