performance - 我怎样才能提高我的 Julia 程序的性能以获得出色的数字？

Question

我一直在玩一些Perl programs to calculate excellent numbers .尽管我的解决方案的运行时间是可以接受的，但我认为另一种语言，尤其是为数字设计的语言，可能会更快。 friend 建议Julia ，但我看到的表现太糟糕了，我一定是做错了什么。我已经浏览了 Performance Tips并没有看到我应该改进什么:

digits = int( ARGS[1] )

const k = div( digits, 2 )

for a = ( 10 ^ (k - 1) ) : ( 10 ^ (k) - 1 )
    front = a * (10 ^ k + a)

    root = floor( front ^ 0.5 )

    for b = ( root - 1 ): ( root + 1 )
        back = b * (b - 1);
        if back > front
            break
        end
        if log(10,b) > k
            continue
        end
        if front == back
            @printf "%d%d\n" a b
        end
    end
 end

我有一个等效的 C 程序，它比 Julia page 上标注的 2 倍快了一个数量级。 (尽管关于 Julia 速度的大多数 Stackoverflow 问题似乎都指出了该页面中存在缺陷的基准测试):

而我写的未经优化的纯 Perl 只需要一半的时间:

use v5.20;

my $digits = $ARGV[0] // 2;
die "Number of digits must be even and non-zero! You said [$digits]\n"
    unless( $digits > 0 and $digits % 2 == 0 and int($digits) eq $digits );

my $k  = ( $digits / 2 );

foreach my $n ( 10**($k-1) .. 10**($k) - 1 ) {
    my $front = $n*(10**$k + $n);
    my $root = int( sqrt( $front ) );

    foreach my $try ( $root - 2 .. $root + 2 ) {
        my $back = $try * ($try - 1);
        last if length($try) > $k;
        last if $back > $front;
        # say "\tn: $n back: $back try: $try front: $front";
        if( $back == $front ) {
            say "$n$try";
            last;
            }
        }
    }

我正在为 Mac OS X 使用预编译的 Julia，因为我无法编译源代码(但我没有尝试超越它的第一次)。我想这是其中的一部分。

另外，我看到任何 Julia 程序的启动时间大约为 0.7 秒(参见 Slow Julia Startup Time )，这意味着等效的编译 C 程序在 Julia 完成一次之前可以运行大约 200 次。随着运行时间的增加(更大的 digits 值)和启动时间意味着更少，我的 Julia 程序仍然很慢。

我还没有涉及非常大的数字(20 位以上的优秀数字)的部分，我没有意识到 Julia 处理这些问题的能力并不比大多数其他语言好。

这是我的 C 代码，它与我开始时有点不同。我生疏的、不优雅的 C 技能与我的 Perl 基本相同。

#include <math.h>
#include <stdio.h>
#include <stdlib.h>

int main( int argc, char *argv[] ) {
    long 
        k, digits,
        start, end,
        a, b,
        front, back,
        root
        ;

    digits = atoi( argv[1] );

    k = digits / 2;

    start = (long) pow(10, k - 1);
    end   = (long) pow(10, k);

    for( a = start; a < end; a++ ) {
        front = (long) a * ( pow(10,k) + a );
        root  = (long) floor( sqrt( front ) );

        for( b = root - 1; b <= root + 1; b++ ) {
            back = (long) b * ( b - 1 );
            if( back > front )  { break; }
            if( log10(b) > k ) { continue; }
            if( front == back ) {
                printf( "%ld%ld\n", a, b );
                }
            }
        }

    return 0;
    }

Answer 1

我已经针对以下代码对您的代码 ( brian.jl ) 进行了基准测试，该代码尝试对您的代码进行最少的更改并遵循 Julian 风格:

function excellent(digits)
    k = div(digits, 2)
    l = 10 ^ (k - 1)
    u = (10 ^ k) - 1

    for a in l:u
        front = a * (10 ^ k + a)
        root = isqrt(front)
        for b = (root - 1):(root + 1)
            back = b * (b - 1)
            back > front && break
            log10(b) > k && continue
            front == back && println(a,b)
        end
    end
end
excellent(int(ARGS[1]))

分出 u和 l是个人对可读性的偏好。作为基准，我机器上的 Julia 启动时间是:

$ time julia -e ''
real    0m0.248s
user    0m0.306s
sys 0m0.091s

因此，如果您从冷启动开始每次执行 Julia 所运行的计算大约为 0.3 秒，那么在此阶段 Julia 可能不是您的最佳选择。我通过了 16到脚本，并得到:

$ time julia brian.jl 16
1045751633986928
1140820035650625
3333333466666668

real    0m15.973s
user    0m15.691s
sys     0m0.586s

和

$ time julia iain.jl 16
1045751633986928
1140820035650625
3333333466666668

real    0m9.691s
user    0m9.839s
sys     0m0.155s

编写的此代码的一个限制是如果 digits>=20我们将超出 Int64 的存储量.出于性能原因，Julia 不会将整数类型自动提升为任意精度整数。我们可以通过将最后一行更改为:

digits = int(ARGS[1])
excellent(digits >= 20 ? BigInt(digits) : digits)

我们得到了 BigInt优秀的免费版本，这很好。暂时忽略这一点，在分析我的版本时，我发现大约 74% 的时间用于计算 log10 ，在 isqrt 上关注约 19% .我通过将最后一行替换为

excellent(4)  # Warm up to avoid effects of JIT
@profile excellent(int(ARGS[1]))
Profile.print()

现在，如果我们想涉足较小的算法更改，鉴于我们现在从分析器中了解到的信息，我们可以替换 log10行(只是检查有效位数)与 ndigits(b) > k && continue ，这给了我们

$ time julia iain.jl 16
1045751633986928
1140820035650625
3333333466666668

real    0m3.634s
user    0m3.785s
sys     0m0.153s

这将余额从 isqrt 更改为约 56%和 ~28% 来自 ndigits .进一步挖掘这 56%，大约一半用于执行 this line这似乎是一个非常明智的算法，因此任何改进都可能会改变比较的精神，因为它确实是一种完全不同的方法。使用 @code_native 调查机器码倾向于暗示没有其他什么太奇怪的事情发生了，尽管我没有深入研究这一点。

如果我允许自己进行一些更小的算法改进，我可以从 root+1 开始并且只做 ndigits检查一次，即

for a in l:u
    front = a * (10^k + a)
    root = isqrt(front)
    b = root + 1
    ndigits(b) > k && continue
    front == b*(b-1) && println(a,b)
    b = root
    front == b*(b-1) && println(a,b)
    b = root - 1
    front == b*(b-1) && println(a,b)
end

这让我明白

real    0m2.901s
user    0m3.050s
sys     0m0.154s

(我不相信需要后两个相等性检查，但我正在努力最小化差异!)。最后，我想我可以通过预计算来提高速度 10^k ，即 k10 = 10^k ，这似乎是在每次迭代时重新计算的。有了这个，我得到

real    0m2.518s
user    0m2.670s
sys     0m0.153s

与原始代码相比，这是一个相当不错的 20 倍改进。