How does setrounding work?

New to Julia
1

This is a session of my REPL

julia> setrounding(BigFloat, RoundDown)
MPFRRoundDown::MPFRRoundingMode = 3 

julia> function f()
        println(big"0.3")
end
f (generic function with 1 method)
 
julia> setrounding(BigFloat, RoundUp)
MPFRRoundUp::MPFRRoundingMode = 2
 
julia> f()
0.2999999999999999999999999999999999999999999999999999999999999999999999999999974
 
julia> function g()
        println(big"0.3")
end
g (generic function with 1 method)
 
julia> g()
0.3000000000000000000000000000000000000000000000000000000000000000000000000000017

Here I make the following conclusion: what determines the rounding behaviour of a function is the rounding mode that was in play at the time the function was created. Changing the rounding mode at a later time doesn’t change the function’s behavior.

Here is another session of the REPL:

julia> setrounding(BigFloat, RoundUp)
MPFRRoundUp::MPFRRoundingMode = 2
 
julia> function h()
        setrounding(BigFloat, RoundDown)
        if true
                println(big"0.3")
        end
end
h (generic function with 1 method)

julia> h()
0.3000000000000000000000000000000000000000000000000000000000000000000000000000017

Here I make the following conclusion: what determined the rounding behavior inside the if’s body was the “global setrounding”, and not the “setrounding” inside h’s body.

Because I am new to Julia, I don’t understand this behavior so well. I need some help to understand the behavior I am noticing here.

2

setrounding is indeed global, if you want it set temporarily, use the form that takes a function so that the original setting will be restored. See Integers and Floating-Point Numbers · The Julia Language.

Also, please quote your code.

2 Likes
3

Thank you for directing me to your post, I will pay attention those advices.

4

A reasonable conclusion given your data, but not entirely correct. Changing the rounding mode will affect the behavior of a function, no matter when it was defined:

julia> function big_pi()
         big(Ď€)
       end
big_pi (generic function with 1 method)

julia> big_pi()
3.141592653589793238462643383279502884197169399375105820974944592307816406286198

julia> setrounding(BigFloat, RoundUp)
MPFRRoundUp::MPFRRoundingMode = 2

julia> big_pi()
3.141592653589793238462643383279502884197169399375105820974944592307816406286233

(note how the last few digits have changed).

So why is your function different? I think what’s going on is that your function does not compute a new BigFloat. Instead, your function definition uses the big"" string macro to create a literal BigFloat, and that literal is part of the function definition. You can see this in the @code_lowered:

julia> @code_lowered g()
CodeInfo(
1 ─ %1 = Main.println(0.3000000000000000000000000000000000000000000000000000000000000000000000000000017)
└──      return %1

The literal 0.300....17 is part of your function definition. Changing the rounding mode later won’t affect that literal.

It’s the same as if you did:

function g()
  println(0.3)
end

Even if you later change how arithmetic on Float64 works (Julia allows you to do this!), nothing will change that 0.3 to some other literal.

When you re-define the function, the big"0.3" string macro is run again, and the rounding mode causes a different literal to be used in the body of your function.

2 Likes
5

Is setrounding an efficient way of controlling the roundiing behaviour? I am asking this because validated numerics is important to me , and I am writing code where I need rigourous bounds to the results of my numerical computations.

I have written a function that takes two BigFloat numbers as input, and rounds the result up

function round_up_sum(x, y)
	setrounding(BigFloat, RoundUp)
	return x + y
end

I am going to create similar functions: round_down_sum, round_up_times, round_down_times, and so on, so that I have control in the way I round my computations.

Is my implementation efficient? Is there a more efficient way of controlling the rounding behavior?

6

GitHub - JuliaIntervals/IntervalArithmetic.jl: Rigorous floating-point calculations using interval arithmetic in Julia might be a good place to start. It will be very efficient (although it might give slightly wide bounds due to using Floats)

2 Likes
7

I am going to read this package documentation, I feel it will help me a lot.

1 Like
8

For the sake of completeness StochasticArithmetic.jl also provides directed rounding for basic operations on FP numbers (currently only Float32 and Float64). Disclaimer: I’m the author of this package and it has not been very thoroughly tested, nor has it been used much since I wrote it. If what you need is interval arithmetic, then you should definitely use IntervalArithmetic.jl.
If, however, you rather need base building blocks, StochasticArithmetic.jl might contain interesting ideas. If you need increased precision, you’ll probably not be able to use it directly; but if you have to implement your own directed rounding, you might want to try its algorithms, which I tried to make as efficient as I could (at least for standard FP numbers; not sure whether that would generalize well to arbitrary precision).

Basically, directed roundings in StochasticArithmetic.jl are computed using an Error-Free Transformation to compute at the same time both a rounded result and the associated error. This first step is followed by a computation of either the predecessor or the successor depending on the error sign and the desired rounding mode.


For comparison, here is the sum of two vectors, rounded upwards using IntervalArithmetic and StochasticArithmetic:

using IntervalArithmetic
function f!(a, b)
    @inbounds for i in eachindex(a, b)
        # note that this also computes the downward rounding
        a[i] = (interval(a[i]) + interval(b[i])).hi
    end
end

using StochasticArithmetic
function g!(a, b)
    @inbounds for i in eachindex(a, b)
        # note that you'd have to add another instruction to compute
        # the downward rounding as well and make the comparison entirely
        # fair
        a[i] = +(UP, a[i], b[i])
    end
end
g(a, b) = +(UP, a, b)

This not an entirely fair comparison, since IntervalArithmetic always computes full intervals, with both the upward and downward rounding, but I think it does give an idea of how both approaches compare for Float64 numbers:

using BenchmarkTools
let n = 1_000
    x = 0.1*rand(n)
    y = 0.3*rand(n)

    x1 = copy(x); f!(x1, y)
    x2 = copy(x); g!(x2, y)
    @assert x1 == x2

    println("IntervalArithmetic");   @btime f!($x, $y)
    println("StochasticArithmetic"); @btime g!($x, $y)
end

IntervalArithmetic
  17.892 ÎĽs (0 allocations: 0 bytes)
StochasticArithmetic
  1.570 ÎĽs (0 allocations: 0 bytes)
3 Likes