[ANN] MPFR_wrap.jl

Community
,
1

Hello everyone,

I would like to announce MPFR_wrap.jl, a set of wrappers around MPFR, a multiple-precision floating-point library for in-place computation. MPFR_wrap.jl does not introduce new types and makes instead use of Julia’s BigFloats.

Example

Most operations defined in Base on BigFloats however allocate temporary variables. Consider the computation of sqrt(a+b): a new BigFloat is allocated to hold the result of a+b. Using MPFR_wrap, one can do

add!(a, a, b) # a ← a + b
sqrt!(a, a) # a ← sqrt(a)

provided the content of a can be overwritten.

Installation

The package is not registered yet: if you want to give it a try, at the package prompt (“]”), type

(@v1.4) pkg> add https://github.com/mzaffalon/MPFR_wrap.jl

This is my first public module: feedback is welcome.

Background

This module was born out of curiosity to implement the computation of pi using iterative methods as explained in this review by David H. Bailey and to compare with the timings of table 3.

The in-place computation of pi to 100_000 digits takes 130ms on my laptop, compared to 350ms using the operations defined in Base; a substantial difference between the two methods can be seen when one consider allocations: 400 kB compared to 12 MB.

Disclaimer: I am not an expert in the computation of pi.

6 Likes
2

Thanks for doing this!

I’m actually happy to see that the non-in-place versions that I get by just using BigFloat with a nice syntax are losing less than a factor of 3 in performance; I thought it might be much more than that.

3 Likes
3

This may have something to do with this commit,

Using that allows us to use our fast memory-pool gc and avoid adding finalizers and use the slow malloc/free functions. (and in some cases, some of it might even end up on the stack!)

although I haven’t tested it.

4

Now all you need is a macro that can convert code written with normal syntax like a += b into in-place code.

1 Like
5

I have been wondering if c ← a + b and a ← sqrt(a) are a better idea, but this is above my pay grade.

With dot notation, I don’t know how to do b ← a - b.

1 Like
6

It would probably be possible to overload the broadcast machinery here, so that a .+= b would do inplace addition. That should even be possible with sqrt, although I am not 100% sure how robust such a solution would be overall.

1 Like
7

It’s fun to be a pirate! :pirate_flag: :laughing:

julia> using MPFR_wrap

julia> Base.copyto!(a::BigFloat, b::Base.Broadcast.Broadcasted{<:Any,<:Any,typeof(+)}) = add!(a, b.args...)

julia> big"1.2" .+= big"2.3"
3.5
1 Like
8

The problem is if you would write

c = a
a += b

since now c gets modified (since it is the same object as a).

2 Likes
9

This is quite nice!

How about

add!(a, b) # same as add!(a, a, b) 
sqrt!(a) # same as sqrt!(a, a)

?

1 Like
10

Taking the broadcasting idea a bit further:

using MPFR_wrap
using Base.Broadcast: Broadcasted

mutating!(::typeof(-), dest, x) = neg!(dest, x)
mutating!(::typeof(abs), dest, x) = abs!(dest, x)
for f in (
    :sqrt, :cbrt,
    :log, :log2, :log10, :log1p, :exp, :exp2, :exp10, :expm1,
    :cos, :sin, :tan, :sec, :csc, :cot, :acos, :asin, :atan,
    :cosh, :sinh, :tanh, :sech, :csch, :coth, :acosh, :asinh, :atanh,
)
    @eval mutating!(::typeof($f), dest, x) = $(Symbol(f, :!))(dest, x)
end

mutating!(::typeof(+), dest, x, y) = add!(dest, x, y)
mutating!(::typeof(-), dest, x, y) = sub!(dest, x, y)
mutating!(::typeof(*), dest, x, y) = mul!(dest, x, y)
mutating!(::typeof(/), dest, x, y) = div!(dest, x, y)
mutating!(::typeof(^), dest, x, y) = pow!(dest, x, y)

function Base.copyto!(dest::BigFloat, b::Broadcasted)
    dest_modified = false
    args = (
       if i isa Broadcasted
           copyto!(dest_modified ? zero(dest) : (dest_modified=true; dest), i)
       else
           i 
       end for i in b.args
    )
    return mutating!(b.f, dest, args...)
end

This should work for nested functions as well:

julia> big"0." .= sqrt.(big"9.") ./ cbrt.(big"8.")
1.5
2 Likes
11

I am mimicking MPFR’s API, which exposes add!(c, a, b), sqrt!(c, a).

Are you suggesting add!(a, b) to avoid typing twice the same variable?

12

Yes. I belive I’ve seen this pattern a few times, that there are two versions, foo!(b, a) and foo!(a), equal to foo!(a, a).

1 Like