Fuss-Catalan数的Julia记忆代码

Question

从Fuss-Catalan系列C{4}_n (请参阅在线整数序列百科全书OEIS A002293)中，我想使用memoization计算第n项。我在下面编写的代码可以工作，但在我的笔记本电脑上运行n=200时需要大约43秒。有没有办法进一步加快速度呢？

numterms = 20
C4 = Array{BigInt}(numterms+1) # memoization dictionary
fill!(C4,-1) # -1 implies not yet computed
C4[1] = 1 # Base case for n=0, C[n+1] provides nth term

function catalan4(n,C)
    C[n+1] == -1 || return C[n+1]
    sum1 = convert(BigInt,0)
    for i in 1:n
       sum2 = convert(BigInt,0)
       for j in 1:(n-i+1)
            sum3 = convert(BigInt,0)
            for k in 1:(n-i-j+2)
                sum3+= catalan4(k-1,C)*catalan4(n-i-j-k+2,C)
            end
            sum2 +=  catalan4(j-1,C)*sum3          
       end
       sum1 +=  catalan4(i-1,C)*sum2 
    end
    C[n+1] = sum1    
    return sum1
end 

for i in 1:numterms
    println(i,"\t",catalan4(i,C4))
end    

这提供了预期的效果：

1   1
2   4
3   22
4   140
5   969
6   7084
7   53820
8   420732
9   3362260
10  27343888
11  225568798
12  1882933364
13  15875338990
14  134993766600
15  1156393243320
16  9969937491420
17  86445222719724
18  753310723010608
19  6594154339031800
20  57956002331347120

谢谢!

Answer 1

所以，在Stefan上面的回答之后，我尝试了几种方法。为了查看这个立方体是否真的是Julia的BigInt，我在上面代码的C版本中尝试了GMP mpz_t整数。实际上，n=200的速度大约快了10倍。然而，我有时需要n=4096，而GMP和C并不能真正帮助我减少计算时间。然后我注意到内部循环也在被反复地重新计算。所以我记住了两个内部和，如下所示，现在代码给了我n=200，用Julia和BigInts的计算时间从43s减少到0.06秒！

numterms = 200
C4 = Array{BigInt}(numterms+1) # memoization dictionary
fill!(C4,-1) # -1 implies not yet computed
C4[1] = 1 # Base case for n=0, C[n+1] provides nth term

# memoize partial sums as well
sum2store = similar(C4)
sum3store = similar(C4)
fill!(sum2store, -1)
fill!(sum3store, -1)

function catalan4(n,C)
    C[n+1] == -1 || return C[n+1]
    sum1 = convert(BigInt,0)
    for i in 1:n
        if sum2store[n-i+1] == -1
            sum2 = convert(BigInt,0)
            for j in 1:(n-i+1)
                if sum3store[n-i-j+2] == -1
                    sum3 = convert(BigInt,0)
                    for k in 1:(n-i-j+2)
                        sum3+= catalan4(k-1,C)*catalan4(n-i-j-k+2,C)
                    end
                    sum3store[n-i-j+2] = sum3
                end
                sum2 +=  catalan4(j-1,C)*sum3store[n-i-j+2]
            end
            sum2store[n-i+1] = sum2
        end
        sum1 +=  catalan4(i-1,C)*sum2store[n-i+1]
    end
    C[n+1] = sum1
    return sum1
end

现在，我不需要Julia的BigInt速度很快--我很高兴我可以使用它。

Answer 2

我怀疑糟糕的BigInt性能是这里的罪魁祸首。您可能想尝试一下Nemo，它使用Flint library来处理任意精度的整数，据我所知效率要高得多。如果你想保持标准的Julia (Nemo是基于Julia的，但在某些语言语义上与Julia有所不同，它是一个计算机代数系统)，并且你的数字被限制为小于2^127，那么你可以尝试使用Int128 -这些可以是堆栈分配的，并保存在寄存器中，不像BigInts，它必须是堆分配的，并且LLVM不知道如何推理(它不能将新BigInts的创建转换为现有LLVM的突变)。创建使用两个/四个Int128值进行运算的自定义Int256和Int512类型可能并不太难，特别是当您只需要支持加法和乘法时。