Option to use system libm.

[ntessore]

Thanks to the helpful pointer in JuliaLang/julia#8869, I added an option for USE_SYSTEM_LIBM=1. Since this makes numerical programs quite a bit faster on Macs, I would even suggest to make this the default behaviour, and add another --with-openlibm flag if necessary.

[staticfloat]

This looks really nice, thank you for taking the time to submit a patch! Just to be certain, do the tests (Base.runtests()) pass on your machine with this enabled? (Edited to match @tkelman's correction)

[ntessore]

I noticed that a test for math fails, namely gamma(1:25) == gamma(Float64[1:25])

It's related to rounding, but I'm not sure I understand what exactly the issue is. Here is the output:

gamma(1:25)             gamma(Float64[1:25])
-------------------------------------------------
1.0                     1.0                     
1.0                     1.0                     
2.0                     2.0                     
6.0                     6.0                     
24.0                    24.0                    
120.0                   120.0                   
720.0                   720.0                   
5040.0                  5040.0                  
40320.0                 40320.0                 
362880.0                362880.0                
3.6288e6                3.6288e6                
3.99168e7               3.991680000000003e7     
4.790016e8              4.790015999999976e8     
6.2270208e9             6.227020799999989e9     
8.71782912e10           8.717829119999985e10    
1.307674368e12          1.3076743679999983e12   
2.0922789888e13         2.0922789888000055e13   
3.55687428096e14        3.556874280960007e14    
6.402373705728e15       6.402373705727994e15    
1.21645100408832e17     1.2164510040883208e17   
2.43290200817664e18     2.43290200817664e18     
5.109094217170942e19    5.109094217170942e19    
1.1240007277776115e21   1.1240007277776115e21   
2.5852016738885247e22   2.5852016738885247e22   
6.204484017332391e23    6.204484017332391e23    

[staticfloat]

Can you post gamma(1:25) - gamma(Float64[1:25]) so we can see what's
going on?

On Sat, Nov 8, 2014 at 12:59 AM, Nicolas Tessore notifications@github.com
wrote:

I noticed that a test for math fails, namely gamma(1:25) ==
gamma(Float64[1:25])

It's related to rounding, but I'm not sure I understand what exactly the
issue is. Here is the output:

gamma(1:25) gamma(Float64[1:25])

1.0 1.0
1.0 1.0
2.0 2.0
6.0 6.0
24.0 24.0
120.0 120.0
720.0 720.0
5040.0 5040.0
40320.0 40320.0
362880.0 362880.0
3.6288e6 3.6288e6
3.99168e7 3.991680000000003e7
4.790016e8 4.790015999999976e8
6.2270208e9 6.227020799999989e9
8.71782912e10 8.717829119999985e10
1.307674368e12 1.3076743679999983e12
2.0922789888e13 2.0922789888000055e13
3.55687428096e14 3.556874280960007e14
6.402373705728e15 6.402373705727994e15
1.21645100408832e17 1.2164510040883208e17
2.43290200817664e18 2.43290200817664e18
5.109094217170942e19 5.109094217170942e19
1.1240007277776115e21 1.1240007277776115e21
2.5852016738885247e22 2.5852016738885247e22
6.204484017332391e23 6.204484017332391e23

—
Reply to this email directly or view it on GitHub
#125 (comment)
.

[tkelman]

I saw that same failure last time I tried USE_SYSTEM_LIBM on Ubuntu. IIRC the relative difference was O(eps)

[ntessore]

$ julia -e 'println(gamma(1:25) - gamma(Float64[1:25]))'
[0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,-2.9802322387695312e-8,2.384185791015625e-6,1.1444091796875e-5,0.000152587890625,0.001708984375,-0.0546875,-0.6875,6.0,-80.0,0.0,0.0,0.0,0.0,0.0]

[staticfloat]

Even the -80 is 14 orders of magnitude smaller than the actual numbers being compared. @ViralBShah is this the expected precision when dealing with system libm?

[staticfloat]

I'm merging this as-is for now, if we decide the precision is good enough, we can alter our tests and use system libm as default. In any case, thank you for your contribution!

[ntessore]

Checking what actually gets called in both cases:

julia> gamma(12)
3.99168e7

julia> gamma(12.)
3.991680000000003e7

julia> @which gamma(12)
gamma(n::Union(Int8,Int16,UInt16,UInt32,Int64,UInt64,UInt8,Int32)) at combinatorics.jl:45

julia> @which gamma(12.)
gamma(x::Float64) at special/gamma.jl:1

Checking combinatorics.jl:45

function gamma(n::Union(Int8,UInt8,Int16,UInt16,Int32,UInt32,Int64,UInt64))
    n < 0 && throw(DomainError())
    n == 0 && return Inf
    n <= 2 && return 1.0
    n > 20 && return gamma(float64(n))
    @inbounds return float64(_fact_table64[n-1])
end

And there we see that the factorial table is in integers, which round more nicely than the returned Float64s from libm.

[ViralBShah]

Cc: @alanedelman @stevengj @jiahao @andreasnoack

We should file this as an issue in base julia, if required.