Rational{BigInt} problems

[Bill Hart]

I'm trying to rewrite my Nemo package to make it conform as much as possible to the conventions of the Julia language. 

On that score I finally relented and decided to stick with the Julia convention of using // to construct fractions, instead of /, and to stick with the (to me) strange definition of div for values of type BigInt{Rational} (historical note: it turns out this is the original definition Euclid used, despite the fact that modern algebraists do not define it that way in a "Euclidean domain", such as the field of rational numbers). I also decided to simply learn to get used to the inconsistent semantics of integer literals as either Int64's, Int128's or BigInt's depending on their size. I want Nemo to be as consistent with the rest of Julia, and its design intentions as possible, even though my application for the language is somehow quite different.

[[ Nemo is a computer algebra package for Julia, which allows constructing various Rings (e.g. polynomial rings, residue rings, finite fields, power series rings, p-adics, (and more recently number fields and order), recursively over each other). ]]

Of course I don't want to compromise on performance (at least, not much) or on neat code.

After reworking a substantial proportion of the Nemo codebase, I've resolved most of the problems, but Rational{BigInt} is currently causing me substantial difficulties which I've been unable to resolve.

In particular, I want Nemo to  make use of the // notation to construct fractions, as per the Julia convention. When using BigInt, this results in a Rational{BigInt} type.

I wonder if I could get some comment on the following issues:

1) The current Rational{BigInt} type appears to not be compatible with the GMP mpq_ type? It looks to me that effectively from a C interface perspective Rational{BigInt} is a type containing two pointers to mpz_structs. Whereas a GMP mpq_t contains two mpz_structs without the extra level of indirection. 

Is this because Julia does not support embedding C structs within its types, and Julia wants Rational{BigInt} to consist of a pair of Julia types, namely BigInt's? Or am I mistaken about this.

This is problematic because I can't call GMP for the mpq_t functionality that I need, at least not without creating an auxiliary mpq type in Julia, copying the contents and adding substantial gc overhead. Or I can reimplement the mpq functions in Julia, which ends up being extremely inefficient due to gc, etc.

2) The Rational{BigInt} type appears to be immutable, which means that it can't be passed as the output to a GMP function (or any C function, for that matter) anyway. This adds overhead for all the same reasons.

3) I can't find a way to add Rational{BigInt} to my FieldElem typeclass in Nemo, as it already belongs to a typeclass (Real or something). This means I have to have T <: Union(Rational{BigInt}, FieldElem) everywhere throughout my code, and it becomes extremely unwieldy. Are there any plans in the pipeline to ease this "posthoc inheritance" issue? It currently makes interfaces hard to extend.

4) Because GC ends up taking most of the time for some computations in Nemo, we came up with some unsafe functions mul!, addeq!, etc., which do what you think they do, to cut down on the number of objects allocated by reusing objects for output. But Rational{BigInt} defeats this optimisation by being immutable. We can't write a mul! for Rational{BigInt}!

This optimisation is not just desirable, it is absolutely fundamental to our entire approach. And it is being defeated by immutability!

5) Things are a little worse for me in that the type Flint type uses for bignum rationals in an fmpq. So I need to convert to and from Rational{BigInt} for every operation involving the Julia Rational{BigInt} type, e.g. doing something with a coefficient of a rational polynomial. The following is currently the best hack I have come up with to handle this (note fmpq is the Flint rational type). I display it here just to show how dire these problems become:

# Not really an Int, but Julia doesn't provide an analogous type which is sometimes

# an Int and other times a pointer to an mpz_struct

type fmpz

   data :: Int 

   fmpz() = new() # Sometimes we don't want flint to initialise, not Julia

end

function init(z::fmpz) 

   z.data = 0 # sometimes we do want to initialise in Julia

end

# Flint needs to clean up auxilliary data, such as mpz_t's it associates with 

# its fmpz's, but can't use a finalizer, in case the fmpz is mocked up (see below)

function clear(z::fmpz)

   ccall((:fmpz_clear, :libflint), Void, (Ptr{fmpz},), &z) 

end

# We want to pass a Julia BigInt to Flint, but they both use an mpz_t for integers

# that don't fit in an Int, and we definitely don't want to copy the data if we only

# want to read it. So we mock up an fmpz for flint. Note Julia must not call the

# fmpz_clear function on a mocked up fmpz, hence the explicit init/clear functions

# above instead of a finalizer.

function fmpz_readonly(x::BigInt)

   temp = fmpz()

   ccall((:fmpz_init_set_readonly, :libflint), Void, (Ptr{fmpz}, Ptr{BigInt}), &temp, &x)

   return temp

end

   

# Convert back to a Julia BigInt

function BigInt(z::fmpz)

   r = BigInt()

   ccall((:fmpz_get_mpz, :libflint), Void, (Ptr{BigInt}, Ptr{fmpz}), &r, &z)

   return r

end

type fmpq

   num::Int # not really an Int, but Julia doesn't provide an analogous type

   den::Int # not really an Int, but Julia doesn't provide an analogous type

   fmpq() = new() # we only sometimes want to initialise the values in Julia itself

end

function init(z::fmpq)

   z.num = 0 # sometimes we do want to initialise in Julia

   z.den = 1

end

function clear(z::fmpq)

   ccall((:fmpq_clear, :libflint), Void, (Ptr{fmpq},), &z) # flint needs to clean up after use

end

# try to avoid copying all the data, since Flint and Julia both use mpz_t's for the

# components of their rational types, and we only want to *read* the values

function fmpq_readonly(x::Rational{BigInt})

   temp = fmpq()

   temp2 = fmpz()

   ccall((:fmpz_init_set_readonly, :libflint), Void, (Ptr{fmpz}, Ptr{BigInt}), &temp2, &x.num)

   temp.num = temp2.data

   ccall((:fmpz_init_set_readonly, :libflint), Void, (Ptr{fmpz}, Ptr{BigInt}), &temp2, &x.den)

   temp.den = temp2.data

   return temp

end

   

# We need to create an actual Julia Rational{BigInt} type to convert from Flint

# fmpq back to Julia. But copying the data is currently extremely problematic.

function Rational(z::fmpq)

   temp = fmpz()

   temp.data = z.num

   a = BigInt(temp)

   temp.data = z.den

   b = BigInt(temp)

   return Rational(a, b)

end

# This is an absolute mess. It should be a single call to a GMP mpq function or a

# simple efficient readonly conversion to flint types, a call to a flint function and

# a simple conversion back to Julia types (saving that final copy is probably not

# possible.

function gcd(a::Rational{BigInt}, b::Rational{BigInt})

   temp1 = fmpq_readonly(a)

   temp2 = fmpq_readonly(b)

   temp3 = fmpq()

   init(temp3)

   ccall((:fmpq_gcd, :libflint), Void, (Ptr{fmpq}, Ptr{fmpq}, Ptr{fmpq}), &temp3, &temp1, &temp2)

   r = Rational(temp3)

   clear(temp3)

   return r

end

This is incredibly painful and inefficient. And yes, I need all of that because of the different things I need to do with fmpz's and fmpq's. This is just a single example. The semantics needs to be the same as the rest of Nemo because of the recursive nature of the data structures in Nemo, so an individual hack is not useful, or we'd need to do the same thing right across Nemo.

I've thought about not using either Julia BigInt's or Rational{BigInt}'s in Nemo as the default bignum types.

But note that I cannot just make the BigInt type in Nemo to be an fmpz (the Flint integer type) since otherwise one has to type ZZ(123.....145) or fmpz(123......145) or something similar to enter BigNums instead of 123......145 (extremely painful for polynomials with large coefficients -- they can currently be cut-n-pasted from other systems, Magma, Sage, Pari, etc., modulo the use of // instead of / for constructing rationals of course), and any BigNum features of Julia would then require back and forth by the user to convert between the two types. I want to make Nemo as Julian as possible and to be able to interface with other packages as easily as possible!

I also cannot redefine the // operator to construct an fmpq instead of a Rational{BigInt} from BigInt's because apparently redefining operators breaks them for everyone else (i.e. all other packages that are subsequently loaded or currently loaded). I want Nemo users to maximally benefit from other packages available for Julia. (For example I tried redefining / for Int's and it broke the Julia profiler!!)

If Julia intends to use GMP mpz_t for BigInt, I **strongly** believe the Rational{BigInt} type should be binary compatible with the GMP mpq_t, and it should be a mutable type.

I realise the whole bignum thing in Julia needs reworking, and this is not intended to be a finalised implementation. But as things currently are, it is causing unnecessary difficulties for me. As it is, any potential solution I come up with causes me to lose something. I either loose:

1) The ability to use integer literals to enter BigInts,

2) The ability to use the // operator to construct Rational{BigInt}'s; or

3) I end up with a lot of extremely ugly, inefficient code that does lots of copying.

I'm really trying hard to stick with Julian conventions as much as possible, but the current implementation makes this incredibly hard.

I welcome any comments/suggestions.

One other suggestion I could make is that modules by default  use the definition of operators from Base, unless they are redefined inside that module or a module that is imported by that module. This would include the REPL user. 

This would allow Nemo to redefine // for BigInt. Users of Nemo would have it redefined for them, including someone using Nemo from the REPL, but no one else would have the operators redefined. In particular, if a module internally uses //, it would use the Base definition, even if the user was using Nemo.

Somehow, I thought this was the case. But I broke the profiler doing this.

I realise most discussion occurs on GitHub tickets. But if possible I would like to first discuss the totality of the issues, not individual issues which isolated, could get closed as invalid or wontfix. Of course I am happy to open individual tickets later once I understand what the current thinking is.

Bill.

P.S: I am prepared to help implement something better for BigInt and Rational{BigInt}. One of my areas of expertise is efficient bignum implementation.

[Andreas Noack]

See https://github.com/JuliaLang/julia/pull/8699#issuecomment-66459451 and https://github.com/JuliaLang/julia/pull/9270#issuecomment-66230923

As you can see, there is some resistance against updating big numbers in place. However, even though the big numbers are immutable, the fields are not, so you should be able to define mutating versions by accessing the fields. The new GC will also handle big numbers better and it will probably be merged soon.

[Bill Hart]

I guess I don't mind BigInt's being immutable, if they are fast for integers that fit into a single word, ** so long as I can still mutate the contents of an mpz_t ** if that happens to be what they contain. 

There is no way to do efficient bignums if large bignums need to be copied every single operation. 

No more than five minutes after I posted the above (after sweating this the whole day) I realised I might be able to do better with the current setup by mutating the BigInt's inside the Rational{BigInt}'s. You thought of this much quicker than I did!

This at least eases some of the issues I am having with 5 in my list above.

What we really need is for a Rational{BigInt} to be an immutable type containing the equivalent of either two Int's or a pointer to an mpq_t (so that the GMP functions and other standard conversion routines of the many other libraries out there that implement rationals, can be used on it), the internals of which are mutable.

What I would be worried about is that a Rational{BigInt} just becomes an immutable pair of BigInt's (themselves either an Int or pointer to mpz_struct), which is certainly not compatible with GMP.

I really hope that the people working on bignums have thought of this. I realise this is hard, and requires compiler support, but it is a standard problem that occurs often with C interfacing.

[Bill Hart]

Looking at those tickets, I think the proposal is to have Immutable, which might have mutable internals, and Atom, which is indivisible and immutable, treated precisely like an Int or other metal backed type.

I think this definitely sounds great, with BigInt and Rational being of the former type. The compiler won't be able to assume that the internal value of an BigInt didn't change if it contains a large integer value, but operations for which no array allocation is required for the internals of a BigInt will be able to be implemented efficiently by placing them on the stack (and simply copying them around by value).

That seems to be precisely is wanted, assuming Rationals are implemented as an immutable struct containing an immutable pointer to a mutable mpq_struct (so that it can be passed to GMP and other libraries that support this interface). Whatever space exists inside the Rational type can be used to represent small rationals if they fit. If the rational value is not big enough to need a heap allocation, none is made, the whole thing is on the stack and there is no GC needed.

If I'm misunderstanding the proposal there, please let me know. Big thumbs up if not.

Bill.

On 16 January 2015 at 16:11, Bill Hart <goodwi...@googlemail.com> wrote:

[Avik Sengupta]

On your point 1,

Rational{BigInt} is  simply the generic Rational type defined as follows:

immutable Rational{T<:Integer} <: Real

what you get is therefore many operations defined on the rational type for many different kinds, but with no assumptions about the internal type parameter, other than it being some kind of an integer. In particular, there is no direct relation to GMP.

If it is easy/fast to convert the two mpz_t values inside Rational{BigInt} to an mpq_t, then you may want to redefine certain operations for Rational{BigInt}.  It that is not the case however, you may be better of defining a completely separate "BigRational <: Real" .

Regards
-
Avik

[Bill Hart]

On 16 January 2015 at 18:37, Avik Sengupta <avik.s...@gmail.com> wrote:

On your point 1,

Rational{BigInt} is  simply the generic Rational type defined as follows:

immutable Rational{T<:Integer} <: Real

what you get is therefore many operations defined on the rational type for many different kinds, but with no assumptions about the internal type parameter, other than it being some kind of an integer. In particular, there is no direct relation to GMP.

Yes, but this is why the // operator should not yield a Rational{BigInt}. You can't get efficient big rationals this way. It's not just a problem that you can't interface with GMP. It's a problem for any other library that is designed to interface with GMP mpq_t types.

The Rational{T} thing is a cool idea, but it doesn't work in practice for composite types, especially ones where people *really* care about performance.

 

If it is easy/fast to convert the two mpz_t values inside Rational{BigInt} to an mpq_t, then you may want to redefine certain operations for Rational{BigInt}.  It that is not the case however, you may be better of defining a completely separate "BigRational <: Real" .

I don't think the latter is a solution for me, since constructing rationals with the // operator then doesn't construct the right kind of rational. 

Doing conversions to mpq_t is problematic because of the cost of conversion. At best you have to copy 96 bytes for a round trip for a two operand function with one return value and allocate multiple temporaries. That's a big overhead in some cases. (I think it may even be worse than this in practice due to deep copies and is certainly terrible at present for other reasons.)

Bill.

[Bill Hart]

Just to make this concrete, here is a C program written using GMP mpq's:

#include "mpir.h"

int main(void)

{

   mpq_t a, b;

   long i;

   mpq_init(a);

   mpq_init(b);

   mpq_set_ui(a, 1, 12345678912347);

   for (i = 0; i < 1000000; i++) {

      mpq_set_ui(b, i, 12345678912347);

      mpq_add(a, a, b);

   }

   gmp_printf("%Qd\n", a);

   mpq_clear(a);

   mpq_clear(b);

   return 0;

}

Here is how long it takes to run:

wbhart@hilbert:~$ time ./testrat

17241362069/425713065943

real    0m0.454s

user    0m0.448s

sys     0m0.008s

Let's do something similar in Julia:

julia> function doit()

          a = 1//BigInt(12345678912347)

          for i in 0:999999

             b = i//BigInt(12345678912347)

             a += b

          end

          return a

       end

doit (generic function with 1 method)

Here is the timing for Julia:

julia> @time doit()

elapsed time: 16.362996735 seconds (800972352 bytes allocated, 35.30% gc time)

17241362069//425713065943

julia> 16.36/0.454

36.0352422907489

[Bill Hart]

Here, for comparison, is a rather more direct approach to GMP, not that the comparison can be made entirely fair:

type MyRat

    data1::Ptr{Void}

    alloc1::Int32

    size1::Int32

    data2::Ptr{Void}

    alloc2::Int32

    size2::Int32

    function MyRat()

       r = new()

       ccall((:__gmpq_init, :libgmp), Void, (Ptr{MyRat},), &r)

       finalizer(r, _mpq_clear_fn)

       return r

    end

    function MyRat(a::Int, b::Int)

        r = new()

        ccall((:__gmpq_init, :libgmp), Void, (Ptr{MyRat},), &r)

        ccall((:__gmpq_set_ui, :libgmp), Void, 

              (Ptr{MyRat}, Int, Int), &r, a, b)

        finalizer(r, _mpq_clear_fn)

        return r

    end

end

function +(a::MyRat, b::MyRat)

    r = MyRat()

    ccall((:__gmpq_add, :libgmp), Void,

          (Ptr{MyRat}, Ptr{MyRat}, Ptr{MyRat}), &r, &a, &b)

    return r

end

function doit()

    a = MyRat(1, 12345678912347)

    for i in 0:999999

        b = MyRat(i, 12345678912347)

        a += b

    end

    return a

end

julia> @time doit()

elapsed time: 3.214663515 seconds (120015168 bytes allocated, 38.73% gc time)

julia> 3.215/0.454

7.081497797356827

[Andreas Noack]

Could you try with an inplace add! as well? Next, could you try the timings with this branch https://github.com/JuliaLang/julia/pull/8699?

[Bill Hart]

Here is a quick implementation with an inplace operator:

type MyRat

    data1::Ptr{Void}

    alloc1::Int32

    size1::Int32

    data2::Ptr{Void}

    alloc2::Int32

    size2::Int32

    function MyRat()

       r = new()

       ccall((:__gmpq_init, :libgmp), Void, (Ptr{MyRat},), &r)

       finalizer(r, _mpq_clear_fn)

       return r

    end

    function MyRat(a::Int, b::Int)

        r = new()

        ccall((:__gmpq_init, :libgmp), Void, (Ptr{MyRat},), &r)

        ccall((:__gmpq_set_ui, :libgmp), Void, 

              (Ptr{MyRat}, Int, Int), &r, a, b)

        finalizer(r, _mpq_clear_fn)

        return r

    end

end

function _mpq_clear_fn(a::MyRat)

    ccall((:__gmpq_init, :libgmp), Void, (Ptr{MyRat},), &a)

end

function +(a::MyRat, b::MyRat)

    r = MyRat()

    ccall((:__gmpq_add, :libgmp), Void,

          (Ptr{MyRat}, Ptr{MyRat}, Ptr{MyRat}), &r, &a, &b)

    return r

end

function add!(a::MyRat, b::MyRat, c::MyRat)

    ccall((:__gmpq_add, :libgmp), Void,

          (Ptr{MyRat}, Ptr{MyRat}, Ptr{MyRat}), &a, &b, &c)

end

function doit()

    a = MyRat(1, 12345678912347)

    for i in 0:999999

        b = MyRat(i, 12345678912347)

        # a += b

        add!(a, a, b)

    end

    return a

end

julia> @time doit()

elapsed time: 1.232903956 seconds (65186360 bytes allocated, 34.14% gc time)

julia> 1.233/0.454

2.7158590308370045

It's probably not really the time we should be comparing with, because to expect Julia to use delayed evaluation for arithmetic is not realistic I think. Still, it makes a fairer comparison with the C code.

I will try the gc branch you suggested, but I have just returned from overseas, so won't get to it immediately. I am very much looking forward to trying that, however!

Bill.

[Bill Hart]

Sorry, I clearly had a mistake in that code. I meant:

julia> function _mpq_clear_fn(a::MyRat)

           ccall((:__gmpq_clear, :libgmp), Void, (Ptr{MyRat},), &a)

       end

_mpq_clear_fn (generic function with 1 method)

julia> @time doit()

elapsed time: 1.649238662 seconds (56000152 bytes allocated, 38.89% gc time)

julia> 1.65/0.454

3.634361233480176

[Steven G. Johnson]

I think the reason for using Rational{BigInt} instead of a custom mpq-based RationalBigint type was simply economy of code.  This way, given a BigInt type we got a rational type (and all the interoperability with other types) "for free" without writing any additional code.

That being said, I think it would be straightforward to use the mpq code with Rational{BigInt} values.

(The Rational{BigInt} type is almost compatible with mpq_t already -- the only problem is that BigInt is a mutable type.  If BigInt were immutable, since it is an exact mirror of the mpz_t struct, the Rational type would be a mirror of the mpq_t struct.)

Basically, you would just define something like:

function +(x::Rational{BigInt}, y::Rational{BigInt})
    # unpack x and y into Julia mirrors x_ and y_ of the GMP __mpq_struct:
    x_ = mpq_struct(x)
    y_ = mpq_struct(y)
    z_ = mpq_struct() # uninitialized mpq for return value of mpq_add
    ccall((:mpq_add,:libgmp), Void, (Ptr{mpq_struct}, Ptr{mpq_struct}, Ptr{mpq_struct}), &z_, &x_, &y_)
    return convert(Rational{BigInt}, z_)
end

The Julia analogue of __mpq_struct to use in the above code would look like:

immutable mpz_struct # just an immutable version of BigInt
   alloc::Cint
   size::Cint
   d::Ptr{Base.GMP.Limb}
end

convert(::Type{mpz_struct}, i::BigInt) = mpz_struct(i.alloc, i.size, i.d)

type mpq_struct
     num::mpz_struct
     den::mpz_struct
     mpq_struct(num, den) = new(num, den)
     mpq_struct(q::Rational{BigInt}) = mpq_struct(num, den)
     mpq_struct() = new() # uninitialized
end

convert(::Type{Rational{BigInt}}, q::mpq_struct) = Rational{BigInt}(convert(BigInt, q.num), convert(BigInt, q.den))

convert(::Type{BigInt}, z::mpz_struct) = .... not sure how to do this nicely...

The last step is a problem since there is currently no "internal" constructor that allows us to construct a BigInt from (alloc,size,d).  That would be trivial to add to Julia, of course, but it would be nicer to just make BigInt immutable and then we'd get mpz-compatibility of Rational{BigInt} for free.

At that point, it would just be a matter of defining all the methods to call mpq for the corresponding Julia functions.

[Steven G. Johnson]

I've opened an issue:

https://github.com/JuliaLang/julia/issues/9826

[Bill Hart]

I wondered whether the immutable types were basically represented as structs. I'd experimented with doing that, but found it only worked in certain circumstances, so abandoned it.

Ultimately I think you really want bignum arithmetic inlined for small values. You could of course store small values in the mpz_struct, but it is more efficient to store small values in a single word with a bit that selects whether the value is an immediate integer or a pointer to an mpz_struct.

This is particularly important for arrays of BigInts. The cache performance increase due to lower memory usage is huge. That's especially important for polynomials and matrices of (small) BigInts.

If the BigInt type were implemented this way, the Rational{BigInt} type would then coincidentally be equivalent to the Flint rational type, but not the GMP mpq_t type. 

On the other hand, I just realised that if the mpq_struct is allocated in one go and the Rational{BigInt} type just contains pointers to each of the consecutive mpz_structs inside that mpq_struct, then you do get GMP compatibility.

I'm not sure if this can be made to work in general.

I'm not sure why you would want an immutable mpz_struct in your example. GMP will modify the mpz_struct when it is the output of a function, so you don't want the compiler assuming it hasn't changed. It may need to reallocate. It seems that all your immutable mpz_struct does is create the requirement for an additional copy of everything.

Bill.

[Bill Hart]

@Steven, I don't think you can make BigInt an immutable struct. GMP needs to modify the contents, including reallocating and changing the size of the integer.

For example,

mpz_add(r, b, c)

could change any of the size, data or alloc fields of r.

It is the pointer to the mpz_struct you want to be immutable, so that you can coopt it for small integers which are not backed by an actual mpz_struct. That ought to be immutable because you don't want the extra indirection to if that you would have if it were an ordinary type

[Steven G. Johnson]

On Saturday, January 17, 2015 at 7:05:09 PM UTC-5, Bill Hart wrote:

@Steven, I don't think you can make BigInt an immutable struct. GMP needs to modify the contents, including reallocating and changing the size of the integer.

You can still make it immutable.   e.g. you could get a return value via:

r = Array(BigInt, 1)
ccall(:mpz_add, ...., r, &b, &c)
return r[1]

although it would be nice to have a better way to pass "output" arguments by reference in ccall (see https://github.com/JuliaLang/julia/issues/2322).

[Jameson Nash]

I have a PR coming soon for that issue which may be highly relevant.

[Bill Hart]

Again, you make a copy of the struct when you return it. I really don't think this will be efficient.

The way we do it in flint is to have a cache of mpz_struct's allocated on the heap. They are passed around by pointer within flint. That means only a single word to copy when returning them from a function (just the pointer). And as I mentioned, only a single word is required for small integers, since we store them as immediate words in the pointer, if they fit.

I mean something like this:

type mpz_t

   alloc::Int32

   size::Int32

   data::Ptr{UInt}

   function mpz_t()

      r = new()

      ccall((:__gmpz_init, :libgmp), Void, (Ptr{mpz_t},), &r)

      return r

   end

end

immutable BigInt

   ptr::mpz_t # or a mangled immediate integer, if the least significant bit of the "pointer" is set

   BigInt() = new(mpz_t())

   BigInt(a::Int) = new(reinterpret(mpz_t, a))

end

Rational{BigInt} is not immediately compatible with mpq_t if you do this. But no one said that the two pointers it would contain needed to point to two nonconsecutive mpz_structs. If they are consecutive then the two together are compatible with GMP's mpq_t. Though one never needs to set the second pointer, as it is never needed.

I almost get away with:

type mpq_t

   alloc1::Int32

   size1::Int32

   data1::Ptr{UInt}

   alloc2::Int32

   size2::Int32

   data2::Ptr{UInt}

   function mpq_t()

      r = new()

      ccall((:__gmpq_init, :libgmp), Void, (Ptr{mpq_t},), &r)

      return r

   end

end

immutable Rational{T}

   ptr1::T # or a pair of immediate integers representing a small rational value

   ptr2::T # ...  Whenever this is actually an mpz_t, it points to the middle of an mpq_struct that ptr1 points to the start of

                        # though it never actually needs to be set in practice. 

   Rational(a::T, b::T) = new(a, b)

   Rational(a::mpq_t) = new(reinterpret(T, a))

   Rational(a::Int, b::Int) = new(reinterpret(T, a), reinterpret(T, b))

end

But the reinterpret requires T to be a bitstype, so I can't quite do the following:

a = BigInt()

b = Rational{BigInt}(mpq_t())

If I change it to this for now:

immutable Rational{T}

   ptr1::T # or a pair of immediate integers representing a small rational value

   ptr2::T # ...  Whenever this is actually an mpz_t, it points to the middle of an mpq_struct that ptr1 points to the start of

                        # though it never actually needs to be set in practice. 

   Rational(a::T, b::T) = new(a, b)

   Rational(a::mpq_t) = new(a)

   # Rational(a::Int, b::Int) = new(reinterpret(T, a), reinterpret(T, b)) # doesn't work

end

then the following compiles:

a = BigInt()

b = Rational{mpq_t}(mpq_t())

Bill.

   

[Bill Hart]

Having said all this, I would be very interested in comparing a very clever implementation along the lines of what you've suggested, Steven, with one along the lines of what I have suggested. I can see a couple of opportunities for optimisation that your method might have that mine might not, and vice versa.

If it is a simple tradeoff, then the simplest approach should probably win out.

The disadvantage of my method is slightly higher overhead for things that are just over 1 word in size. However, in most applications, it is the very smallest sizes you want the greatest optimisation for. There are exceptions. Relatively sparse multivariate polynomial arithmetic seems to like the case of three words optimised, so it can accumulate a sum of products of one word integers without overflow.  

My approach certainly relies on incremental gc. Otherwise the gc works against the caching.

On the other hand, if gc is the bottleneck for small operations, then your approach might win hands down.

Bill. 

[Bill Hart]

I notice the gc branch got merged with master, so I decided to try it out. Unfortunately I've been unable to build Julia since doing git pull.

I'm using the very latest llvm-svn (so I can use the Cxx package, which I'm using for a project I'm working on). I patched the lldb makefile to get around one build failure I'm familiar with, but now I have another:

codegen.cpp: In function âvoid jl_init_codegen()â:

codegen.cpp:5117:72: error: âconst class llvm::TargetSubtargetInfoâ has no member named âgetDataLayoutâ

     engine_module->setDataLayout(jl_TargetMachine->getSubtargetImpl()->getDataLayout());

Does anyone know how to sort this one out?

[Bill Hart]

I simply removed the line of code that was causing Julia to fail to build, for now.

Anyhow, the times from earlier in the thread went down as follows with the new gc:

16.36s => 6.28s

3.215s => 1.47s

1.65s => 0.85s

Certainly a big improvement. The generic 6.28s is still nowhere near the 0.45s of C though. So I still maintain we need to change the rational bignum format to mpq_t (at least).

[Jameson Nash]

If you're going to live on the llvm-svn branch, you're probably going to need to learn to deal with build failures like this quite frequently. This was broken (by an upstream change to llvm) about 5 hours before you sent your email. Fortunately, they also just branched 3.6.0, so you can switch to that if you don't want to avoid dealing with upstream breakage.

[Bill Hart]

Thanks for looking into it. 

Alas, I'm not living on that branch by choice. Cxx seems to need functionality there, for now. And locally we are investigating Cxx quite seriously for a Julia project we are working on, which wraps a C++ project.

It sounds like things might settle down soon anyway, as Cxx can surely pin to LLVM 3.6.

[Tim Holy]

I think you're the one who has to pin LLVM to 3.6, in your Make.user file.

--Tim

[Bill Hart]

Thanks for reminding me of this. 

What I really meant is that hopefully Cxx won't continue to require llvm-svn.