特定值的幂运算符问题

Question

我正在尝试做一个简单的函数来检查阶乘和斯特林近似之间的差异：

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)

当i=16和i=24时英镑的结果是0！。所以，我检查了两个值的幂，结果是 0:

julia> 16^16
0

julia> 24^24
0

我在 R 中编写了相同的代码，没有任何问题。我做错了什么或者我对 Julia 不了解，我可能应该了解什么？

Answer 1

根据 Integers and Floating-Point Numbers 的 Julia 文档，Julia 整数似乎是 32 位或 64 位，具体取决于您的系统。你的求幂溢出了你的值，即使它们是 64 位的。

Julia 看起来支持 Arbitrary Precision Arithmetic，您需要用它来存储较大的结果值。

根据 Overflow Section，写 big(n) 使得 n 任意精度。

Answer 2

虽然问题已经在另一个 post 上得到了回答，但还有一件事值得一提。

Julia 是为数不多的允许您定义自己的原始类型的语言之一 - 因此您仍然可以使用快速固定精度的数字来处理巨大的值。有一个包 BitIntegers：

BitIntegers.@define_integers 512

现在您可以：

julia>  Int512(2)^500
3273390607896141870013189696827599152216642046043064789483291368096133796404674554883270092325904157150886684127560071009217256545885393053328527589376

通常，即使是大的定点算术数，您也会获得更好的性能。例如：

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

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

Answer 3

有一个简单的解决方案可以解决您的问题，它不涉及使用 BigInt 或任何专门的数字类型，而且速度要快得多。只需稍微调整您的数学表达式即可。

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

好的，我们来测试一下：

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

让我们尝试计时：

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

编辑：好像是

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

可能是更好的解决方案，因为数值更接近原始值。

特定值的幂运算符问题

Problem with power operator for specific values

math

julia