[RFC] MPFR and MPC on Base

[andrioni]

Hello,

as it was discussed in #2564, I am doing this pull request with a tentative integration of both MPFR and MPC on Base. It is currently working and passes all tests (I used the same tests I wrote when they were in packages), so it's easy for anyone to try it.

Known issues:

• make check for MPFR and MPC runs their test suites, lengthening considerably the build process, as make check is always called, and I have yet to discover how to bypass this.
• There is no written documentation currently, nor examples other than the tests.
• Support for different rounding modes.
• Renaming MPFRFloat to BigFloat.
• Renaming precision-related functions
• convert methods for BigComplex
• Migration to a direct representation of the structs involved, instead of a bunch of memory allocated as an array
  • BigInt
  • BigFloat
• mpfr_hypot
• Better basic operations rules between BigTypes and primitives to avoid using promotion to deal with it
• Set 256 as the default precision.
• Add the needed -dev libraries to .travis.yml
• Reimplement sum for MPFRFloat in a way that works, as it needs an array of pointers.
• MPFR has no make distclean option, it seems.

[ViralBShah]

@staticfloat Do you think we can support these in travis before we merge this?

[stevengj]

The types should probably be called BigFloat (replacing the current type of the same name) and BigComplex, rather than MPFRFloat and MPCComplex. (Which libraries are used to implement a given functionality should be mostly invisible to the user.)

Correspondingly, I would use get_bigfloat_precision etcetera to get/set the precision. The short names prec and prec2 should probably be ``get_bigcomplex_precision(x::BigComplex)(and just dompc_get_prec2`), or have a generic `get_precision` that works for all floating-point types.

You don't need iscomplex(::MPCComplex) = true since true is the default for any subtype of Number.

Why isn't there a convert(::Type{Complex{Float64}}, x::MPCComplex) function etcetera?

[andrioni]

I will do the replacements and renamings right now.

The convert functions are lacking mainly because I was finishing them in MPFR (as MPC is just a library around a pair of MPFR floats). Implementing them will be quite direct now.

I'll add these issues to the checklist on the top, and update it as I do them, thanks!

[andrioni]

Also, one question: I'm currently using 53 bits as the default precision (i.e. a Float64), but BigFloat currently uses the GMP default, 64 bits. Should I up the default to match it (for compatibility), or should I use a higher number (maybe 256), to give a nice precision by default for those who don't want to mess with it?

[stevengj]

I would default to 256-ish (~octuple) precision.

It would be much nicer to implement BigComplex as Complex{BigFloat}, so that you share all of the Complex{T} methods. You can still call MPC for the operations as needed. In fact, if you do it right (make BigFloat an immutable wrapper around a single mpfr_t pointer), Complex{BigFloat} should even be bit-compatible with mpc_t, so you won't even need to do a conversion (although making an mpc_t on the fly as needed should be cheap).

I'm confused: why is MPFRFloat.mpfr a Vector{Int32}, instead of being an opaque mpfr_t pointer?

[StefanKarpinski]

I'm confused: why is MPFRFloat.mpfr a Vector{Int32}, instead of being an opaque mpfr_t pointer?

I don't know if it's the case here, but in cases where a fixed amount of memory is used to store some C blob, you can make the thing serializable by using a Julia array for the storage. E.g. that's how compiled PCRE regex patterns work: https://github.com/JuliaLang/julia/blob/master/base/regex.jl#L10. Not sure if that's what's going on here though.

[stevengj]

@StefanKarpinski, that doesn't seem to be what's going on here. mpfr_t in C is an array (pointer) of one __mpfr_struct, which consists of three integers plus a pointer. MPFR.jl uses a Vector{Int32} of length 5 instead, where I guess 5 comes from the number of Int32s you need to store 3 32-bit int variables and a 64-bit pointer?? One problem is that 5 seems wrong to me: with struct padding, it seems like more space is required on a 64-bit architecture, and moreover one of the 3 integer fields is a long which is 64 bits on 64-bit Unix systems. In any case, the naive serialize function won't work for this because you can't just serialize the value of the pointer, you have to serialize the contents of the pointed-to array. So, it needs its own serialize method anyway.

The problem here is that MPFR's interface, now that I look at it, requires the caller to know the layout (or at least the size) of the __mpfr_struct, since it doesn't provide a way to allocate an opaque mpfr_t pointer.

There are two possibilities that occur to me. One is to just make BigFloat a type with some padding that is at least as big as __mpfr_struct, i.e. type BigFloat; pad0::Int32; pad1::Int32; .......; end. The other (better, I think) is to try and match the actual size by copying the original struct definition from mpfr.h

type BigFloat
     prec::Cint # problem: may be a Cshort in some MPFR compilations
     sign::Cint
     exp::Clong
     d::Ptr{Void} # pointer to mp_limb_t, which can vary in size depending on how GMP is compiled
end

However, since this will be part of the Julia build, it should be straightforward to extract the correct types from the header files to insert into the .jl file (e.g. you can run a little C program to print the sizeof values, and there are even tricks to get sizeof that only require running the compiler so that it works when cross-compiling).

(It shouldn't be immutable since the contents are, in fact, mutable by the C API.) This would avoid the double allocation (where you now allocate a type whose contents are a pointer to an array that contains another struct). Then you would pass a Ptr{BigFloat} (via &) to MPFR functions and it should work, assuming the struct types match. Getting Complex{BigFloat} to work with MPC is a more tricky, however, since I don't think Complex{BigFloat} will then have the same layout as a struct of two __mpfr_structs. (It would if BigFloat were immutable, but BigFloat must be mutable.)

[stevengj]

Also, it's not completely clear to me that we need MPC. If we just do Complex{BigFloat}, the generic Complex{T} methods seem like they provide pretty much all of what MPC provides.

(There doesn't seem to be anything too special in the MPC implementations; they just convert everything into the corresponding operations on the real and imaginary parts. They might be more clever about minimizing the number of temporary values compared to Complex{T}. It would be interesting to have a benchmark or two.)

[ViralBShah]

@stevengj That is useful to know. In that case, let's get MPFR in shape. It would be nice to leverage Complex{BigFloat} as much as possible and avoid adding another library.

[stevengj]

PS. It looks like you need to wrap the mpfr_hypot function. (This is necessary for abs(::Complex{BigFloat}) to work.)

[staticfloat]

Everything looks good. Once this lands, we can make some very simple changes to .travis.yml. No packaging required on my end! :D

I suppose we could also include the travis changes in this pull request..... that might actually be the best plan, as then we get to test the travis integration with the pull request itself.

[stevengj]

Another TODO item: BigFloat really shouldn't use promotion rules to perform e.g. *(::BigFloat,::Float64), since converting the second argument to a BigFloat before multiplying is heinously inefficient compared to calling mpfr_mul_d.

This is also a flaw of our BigInt implementation, by the way, since GMP provides special functions to add a BigInt to a Clong etcetera that we aren't exploiting (see #2839).

[stevengj]

I notice that the BigInt implementation in bigint.jl pulls the same ugly stunt: it allocates an Array(Int32, 4) to hold an __mpz_struct, which should really be a

type BigInt
    alloc::Cint
    size::Cint
    d::Ptr{Void}
 end

to match the C struct and avoid the double allocation. Again, I think it should be mutable (since the entries are mutable in the C API).

[andrioni]

I wrote a comment here, but you all basically said everything I was going to say and more, so thanks for the help.

I was trying to make the suggestion work with BigInt as preliminary test, but I got a bunch of stack overflows randomly, and then I wasn't able to compile Julia anymore because it wasn't able to read libreadline.a and some other files, so I'm doing a fresh install of Julia and updating gcc/glibc on this computer to work on it again.

[andrioni]

@staticfloat @ViralBShah the build I did for this pull request is already with a modified .travis.yml to use the system MPFR and MPC (they are needed for newer GCC versions, so they are bundled with it in almost every Linux distro)

[staticfloat]

hahaha, so I ended up putting USE_SYSTEM_MPFR in there twice! Awesome. Nice work, @andrioni. I would only ask you to include libmpfr-dev in the apt-get install line, as many users use .travis.yml as a template for an easy Ubuntu install, and if they don't have mpfr, they will need that line too.

[andrioni]

BigInt now uses the approach proposed by @stevengj, and the .travis.yml now installs the dependencies correctly.

[andrioni]

__GMP_MP_SIZE_T_INT only gets set on some Cray architectures, so _MPFR_PREC_FORMAT gets set to 3 and so both the default precision and the exponent variables are long, if my interpretation of mpfr.h is right.

(by the way, kudos to the one who implemented Cint and Clong types.)

[stevengj]

pointer_from_objref won't work because it returns a jl_value_t*, which has some header information before the member data starts. @vtjnash, will #2818 fix this (so that pointer_from_objref(x::Foo) returns a pointer that can be passed as a C struct corresponding to Foo)?

[vtjnash]

Actually the purpose of #2818 is to help eliminate uses of pointer_from_objref. In this case, it would make the following work: ptrarr = [x.mpfr for x in arr] (although the previous suggestion would also be made to work, but be more difficult from a garbage collection and usage standpoint)

[stevengj]

@andrioni, note that you shouldn't need to have a separate mpfr_struct type. Just put the members directly in BigFloat:

type BigFloat{N} <: FloatingPoint
    prec::Clong
    sign::Cint
    exp::Clong
    d::Ptr{Void}
end

Similarly, BigInt shouldn't need a separate mpz_struct.

[andrioni]

I was having some issues with constructors before, but then again, it seems a recompile after make distclean solved all my stack overflows. I'm updating it again right now.

[andrioni]

Program received signal EXC_BAD_ACCESS, Could not access memory.
Reason: KERN_PROTECTION_FAILURE at address: 0x00007fff5f3ffff8
0x0000000101f5dc48 in ?? ()

Ok, I got more stack overflows. Here's the code I tried to use. I'll keep improving this using the auxiliary mp*_struct types, and when I find a solution for these stack overflows (basically, the constructor seems to call itself endlessly if you define the type and some functions for it, but not if there aren't extra functions).

[stevengj]

@JeffBezanson, do you see why this version of the BigInt() constructor would lead to a stack overflow?

[StefanKarpinski]

Passing &b to gmpz_init where b is the Julia BigInt object seems like a bad idea. That pointer is going to be a jl_value_t* and include the box (type tag, etc.) in front of the actual value. I'm not sure what the best way to do this is, but you need to get a pointer to the data of the value, not the value itself. Could be as simple as adding 1 to the pointer, but that feels kind of hacky. Also, @vtjnash's "recurse types" branch might change this by making Julia reference point at the beginning of the data of Julia values instead of the box.

[andrioni]

Why it works then when I use the auxiliary mpz_struct, shouldn't it be boxed as well?
(I had a similar problem for the mpc_struct, as using two mpfr_struct as members directly doesn't work, but hand 'unrolling' them (i.e. putting the mpfr_struct fields directly, as in f17210d242597cc64b1037776a5318e15d0647ee) works fine.)

[StefanKarpinski]

If that works, then maybe I'm wrong about this. The & syntax may provide a pointer to the data rather than the box. In that case the issue must be something else.

[andrioni]

Can someone comment on the rounding modes? Currently I'm using constants in the MPFR module, should I rename/export them?

Other than that, I think it is finally reasonably finished. I've also rebased it (up to 9a72fab), so the merge should be simple.

[StefanKarpinski]

See #2976 – we don't currently expose rounding modes, but we should have a uniform interface across IEEE 754 and other types of floats for rounding modes. So I would say hold off exposing if for now and leave it as part of #2976.

[JeffBezanson]

Why are we now back to using separate BigInt and mpz_struct types? Not a big deal, I guess I can change it myself.
Also, unfortunately there is a new merge conflict.

[andrioni]

Ahem, I unfortunately destroyed the changes on the repo trying to merge it on my lab computer, I'll restore it as soon as I get back home.

[andrioni]

Which should be in twenty minutes.

[andrioni]

Restored now. (that's why one should never use rebase to rewrite history nor try to manually push changes on a Friday afternoon :))

[JeffBezanson]

Awesome! Hoping to merge this today.

[andrioni]

Feel free to assign me to any issues related to BigFloat or BigInt, I'm willing to support them and help as much as I can :)

[JeffBezanson]

Thanks, and thanks again for doing this. After a couple small tweaks it's looking good!

[ViralBShah]

@andrioni It would be great if you could add the documentation for BigFloat and even review the BigInt documentation. For the rest of the issues that are not addressed here (make check, sum, etc.), I think we should open new issues so that these are not lost.

Fantastic work.