'Problem with power operator for specific values

I am trying to do a simple function to check the differences between factorial and Stirling's approximation:

using DataFrames

n = 24
df_res = DataFrame(i = BigInt[],
                    f = BigInt[],
                    st = BigInt[])

for i in 1:n
    fac = factorial(big(i))
    sterling = i^i*exp(-i)*sqrt(i)*sqrt(2*pi)
    res = DataFrame(i = [i],
                    f = [fac],
                    st = [sterling])
    df_res = [df_res;res]
end

first(df_res, 24)

The result for sterling when i= 16 and i= 24 is 0!. So, I checked power for both values and the result is 0:

julia> 16^16
0

julia> 24^24
0

I did the same code in R, and there are no issues. What am I doing wrong or what I don't know about Julia and I probably should?

mathjulia


Solution 1:[1]

It appears that Julia integers are either 32-bit or 64-bit, depending on your system, according to the Julia documentation for Integers and Floating-Point Numbers. Your exponentiation is overflowing your values, even if they're 64 bits.

Julia looks like it supports Arbitrary Precision Arithmetic, which you'll need to store the large resultant values.

According to the Overflow Section, writing big(n) makes n arbitrary precision.

Solution 2:[2]

While the question has been answered at the another post one more thing is worth saying.

Julia is one of very few languages that allows you to define your very own primitive types - so you can be still with fast fixed precision numbers yet handle huge values. There is a package BitIntegers for that:

BitIntegers.@define_integers 512

Now you can do:

julia>  Int512(2)^500
3273390607896141870013189696827599152216642046043064789483291368096133796404674554883270092325904157150886684127560071009217256545885393053328527589376

Usually you will get better performance for for even big fixed point arithmetic numbers. For an example:

julia> @btime Int512(2)^500;
  174.687 ns (0 allocations: 0 bytes)

julia> @btime big(2)^500;
  259.643 ns (9 allocations: 248 bytes)

Solution 3:[3]

There is a simple solution to your problem that does not involve using BigInt or any specialized number types, and which is much faster. Simply tweak your mathematical expression slightly.

foo(i) = i^i*exp(-i)*sqrt(i)*sqrt(2*pi)  # this is your function
bar(i) = (i / exp(1))^i * sqrt(i) * sqrt(2*pi)  # here's a better way

OK, let's test it:

1.7.2> foo(16)
0.0  # oops. not what we wanted

1.7.2> foo(big(16)) # works
2.081411441522312838373895982304611417026205959453251524254923609974529540404514e+13 

1.7.2> bar(16)  # also works
2.0814114415223137e13

Let's try timing it:

1.7.2> using BenchmarkTools

1.7.2> @btime foo(n) setup=(n=16)
  18.136 ns (0 allocations: 0 bytes)
0.0

1.7.2> @btime foo(n) setup=(n=big(16))
  4.457 ?s (25 allocations: 1.00 KiB) # horribly slow
2.081411441522312838373895982304611417026205959453251524254923609974529540404514e+13

1.7.2> @btime bar(n) setup=(n=16)
  99.682 ns (0 allocations: 0 bytes) # pretty fast
2.0814114415223137e13

Edit: It seems like

baz(i) = float(i)^i * exp(-i) * sqrt(i) * sqrt(2*pi)

might be an even better solution, since the numerical values are closer to the original.

Sources

This article follows the attribution requirements of Stack Overflow and is licensed under CC BY-SA 3.0.

Source: Stack Overflow

Solution Source
Solution 1 rgettman
Solution 2 Przemyslaw Szufel
Solution 3

Related Questions