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?

CodePudding user response：

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.

CodePudding user response：

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)