[為什麼學習Julia]Julia性能測試之「冰雹數」
本專欄由"MATLAB"更名為"MATLAB or Julia", 注意這裡的"or"並不是"live or die"裡面做選擇的意思, 而是數學裡面的並集的意思.
然後我搜索了一下關鍵詞"MATLAB", 發現我的專欄由原先的第二, 下降到了第七, 第六名是"1關注, 0文章"的專欄, 知乎的排名規則很迷啊!!!
搜索了一下關鍵詞"Julia", 發現我的專欄竟然排到了第一, 儘管我的專欄沒有寫過一篇Julia主題的文章,知乎的排名規則很迷啊!!!
為了對得起這個Julia排名, 我決定寫一篇關於Julia的文章. 先聲明啊, 我是從11月28日開始學Julia的, 寫的Julia程序肯定是比較幼稚的, 希望熟悉Julia的讀者包涵與指正.
想想我學Julia的初衷是因為看到了一個Julia的官方性能測試報告, 我被它的高速與代碼簡潔所吸引, 然後就開始了學習Julia之旅. 既有C/C++的性能, 又有MATLAB或Python的代碼簡潔度, 這完全是我夢寐以求的編程語言啊!
性能測試的鏈接:
Julia Micro-Benchmarks
Julia排名第二, 而MATLAB的排名是中間靠後的, MATLAB的這個排名符合我的使用經驗.
然後一對比Julia與MATLAB的測試代碼, 我被震驚到了, 兩者的代碼相似度極高!
也就是說, 如果我會使用MATLAB, 那麼我能快速學習使用Julia.
成本低(學習成本), 收益很高(因為運行速度快), 有什麼理由不學Julia呢?
以上是我學習Julia的理由.
-----------------------正題分割線-----------------------------------------------------
作為一個做事嚴謹的人(疑心病患者晚期)對Julia官網的benchmark的還是有懷疑的, 懷疑官方特意挑選了對Julia有利的測試例子, 又或者故意將其他語言的速度拖慢.
於是, 我想要對Julia的性能進行親手動手測試.
測試例子怎麼選呢?
我首先想到就是Project Euler第14題.
Problem 14 - Project EulerLongest Collatz sequence
The following iterative sequence is defined for the set of positive integers:n → n/2 (n is even)n → 3n + 1 (n is odd)Using the rule above and starting with 13, we generate the following sequence:13 → 40 → 20 → 10 → 5 → 16 → 8 → 4 → 2 → 1It can be seen that this sequence (starting at 13 and finishing at 1) contains 10 terms. Although it has not been proved yet (Collatz Problem), it is thought that all starting numbers finish at 1.
Which starting number, under one million, produces the longest chain?NOTE: Once the chain starts the terms are allowed to go above one million.
這題我思考了很久, 代碼也優化了很多次, 還在知乎上進行了討論:
Matlab做project euler 14怎樣提速？
@Falccm 的代碼差不多是MATLAB做這題的性能極限了.(我試圖提速他的這個代碼, 每次都無法成功)
於是, 我將他的代碼作為MATLAB的benchmark.
Falccm的MATLAB代碼:
function [f,lmax] = pe14(n)
[l(n), l(1), lmax] = deal(0, 1, 0);
for k = 2:n
a = k; lk = 0;
while a >= k
lk = lk + 1;
if a == floor(a/2)*2, a = a*.5;
else a = a*3 + 1; end
end
lk = lk + l(a);
l(k) = lk;
if lk > lmax, lmax = lk; f = k; end
end
測試一下速度(注意, 這裡用的是1e8, 而不是原題要求的1e6, 原因是後者運算時間大概為0.06 秒, 太短了, 每次運行的相對浮動較大):
在MATLAB的R2018b上運行.
>> tic;
[f,lmax] = pe14(1e8)
toc;
f =
63728127
lmax =
950
時間已過 4.134375 秒。
我將Falccm的MATLAB直接改寫為Julia語言:
function euler14(N::Int)
len_star = 1
n_star = 1
dp = zeros(Int16, N)
dp[1] = 1
for n=2:N
t = n
cnt = 0
while t >= n
if t % 2 == 0
t = div(t, 2)
else
t = 3*t+1
end
cnt += 1
end
dp[n] = dp[t] + cnt
if dp[n] > len_star
len_star = dp[n]
n_star = n
end
end
n_star, len_star
end
測試(在Julia 1.0.1上運行):
julia> @time euler14(10^8)
1.920117 seconds (10 allocations: 190.735 MiB, 1.67% gc time)
(63728127, 950)
julia> using BenchmarkTools
julia> @benchmark euler14(10^8)
BenchmarkTools.Trial:
memory estimate: 190.74 MiB
allocs estimate: 5
--------------
minimum time: 1.880 s (0.62% GC)
median time: 1.915 s (3.53% GC)
mean time: 1.915 s (3.00% GC)
maximum time: 1.951 s (4.76% GC)
--------------
samples: 3
evals/sample: 1
Julia的速度是MATLAB的4.134375/1.920117 = 2.15倍.
代碼複雜度差不多, 卻獲得了速度的提升.
然後我使用位運算, 以及其他一些技巧, 進一步提速.
function euler0014(N::Int)
len_star = 1
n_star = 1
dp = zeros(Int16, N)
dp[1] = 1
dp[2] = 2
for n=3:2:N
t = n
t += (t + 1) >> 1
cnt = 2
while t >= n
if t & 1 == 0
t >>= 1
cnt += 1
else
t += (t + 1) >> 1
cnt += 2
end
end
dp[n] = dp[t] + cnt
dp[n+1] = dp[div(n+1, 2)] + 1
if dp[n] > len_star
len_star = dp[n]
n_star = n
end
end
n_star, len_star
end
測試(在Julia 1.0.1上運行):
julia> @time euler0014(10^8)
1.438625 seconds (10 allocations: 190.735 MiB, 2.77% gc time)
(63728127, 950)
julia> @benchmark euler0014(10^8)
BenchmarkTools.Trial:
memory estimate: 190.74 MiB
allocs estimate: 5
--------------
minimum time: 1.375 s (0.83% GC)
median time: 1.438 s (2.92% GC)
mean time: 1.432 s (3.19% GC)
maximum time: 1.477 s (4.75% GC)
--------------
samples: 4
evals/sample: 1
Julia的速度是MATLAB的4.134375/1.438625 = 2.87倍.
千萬別問我為什麼不對MATLAB版的程序進行位運算提速, 讀者有興趣動手試一下就明白了, MATLAB的位運算速度感人.
最後, 上大殺器: C++.
使用的是上面知乎討論帖 @rsa 的C++代碼, 我添加了計時功能, 並且修改了2處bug(變數"l"的初始值由0改為了1, "f[1]"的初始值由0改為了1):
#include<cstdio>
#include<time.h>
const int n = 100000000;
int f[n + 1];
int main() {
clock_t tic, toc;
tic = clock();
int a = 1, l = 1;
f[1] = 1;
for (int i = 2;i <= n;i++) {
long long x = i;
int t = 0;
while (x >= i) {
if (x & 1)x += x + 1 >> 1, t += 2;
else x >>= 1, t++;
}
f[i] = t += f[x];
if (t>l)a = i, l = t;
}
toc = clock();
printf("%d, %d
", a, l);
printf("Use Time:%f seconds
", 1.0 * (toc - tic) / CLOCKS_PER_SEC);
getchar();
return 0;
}
演算法和我優化的Julia代碼差不多, 都考慮使用了位運算, 奇數時, 直接進行兩步(因為奇數變換後, 肯定是偶數, 可以再除以2).
運行結果(VS2017社區版, release模式下):
發現這個例子裡面Julia的運算速度比C++的還要快!
為了嚴謹, 我將我優化的Julia演算法直接修改為C++代碼:
#include<cstdio>
#include<time.h>
const int n = 100000000;
int f[n + 10];
int main() {
freopen("output.txt", "w", stdout);
clock_t tic, toc;
int a = 1, l = 1;
f[1] = 1;
f[2] = 2;
tic = clock();
for (int i = 3;i <= n;i+=2) {
long long x = i;
x += (x + 1) >> 1;
int t = 2;
while (x >= i) {
if (x & 1) {
x += (x + 1) >> 1;
t += 2;
}
else {
x >>= 1;
t++;
}
}
f[i] = t + f[x];
f[(i + 1)] = f[(i + 1) / 2] + 1;
if (f[i] > l) {
a = i;
l = f[i];
}
}
toc = clock();
printf("%d, %d
", a, l);
printf("Use Time:%f seconds
", 1.0 * (toc - tic) / CLOCKS_PER_SEC);
return 0;
}
運行結果為:
63728127, 950
Use Time:1.869000 seconds
運行時間總結(單位是秒):
- MATLAB: 4.134375
- Julia(沒有位運算): 1.920117
- Julia(位運算): 1.438625
- C++(rsa的代碼): 2.18
- C++(我寫的代碼): 1.869
為什麼Julia會比C++快呢?
我也搞不太懂, Julia我是初學, C++也是個半吊子, 不好意思. 是因為自帶了並行計算? 好像說不通, 因為這個演算法的循環當前值是依賴於歷史值的, 屬於迭代演算法, 並不能並行計算.
知道答案的讀者可以在評論區留言.
經過一些網友的測試, 發現他們用其他C++的編譯器運行速度是要快於Julia的, 所以, 我這運行速度慢, 應該是VS2017的鍋, 或者是我的配置不太對.
總結:
Julia的優勢是既短又快.
本文通過一個簡單的例子說明了這點.
本來是想讓Julia和MATLAB比較速度的, 結果發現比C++還快, 出乎的我的意料.
---2019年1月1日更新------------------------------------------------
使用宏命令"@code_warntype", 會發現變數"len_star"的類型並不"stable"
因此, 將代碼改為如下:
function euler0014(N::Int)
len_star = Int16(1)
n_star = 1
dp = Vector{Int16}(undef, N)
dp[1] = 1
dp[2] = 2
for n=3:2:N
t = n
t += (t + 1) >> 1
cnt = Int16(2)
while t >= n
if t & 1 == 0
t >>= 1
cnt += Int16(1)
else
t += (t + 1) >> 1
cnt += Int16(2)
end
end
dp[n] = dp[t] + cnt
dp[n+1] = dp[div(n+1, 2)] + Int16(1)
if dp[n] > len_star
len_star = dp[n]
n_star = n
end
end
n_star, len_star
end
主要變化為: 將一些字面整數轉化為Int16類型的, 因為在64為操作系統裡面, 字面整數默認為Int64類型.
改完以後, 在使用"@code_warntype"檢查一下, 發現類型全部都是"stable"了.
性能測試:
julia> @benchmark euler0014(10^8)
BenchmarkTools.Trial:
memory estimate: 190.73 MiB
allocs estimate: 2
--------------
minimum time: 1.201 s (0.98% GC)
median time: 1.232 s (1.11% GC)
mean time: 1.236 s (3.06% GC)
maximum time: 1.273 s (6.39% GC)
--------------
samples: 5
evals/sample: 1
添加宏命令"@inbounds":
function euler0014(N::Int)
len_star = Int16(1)
n_star = 1
dp = Vector{Int16}(undef, N)
dp[1] = 1
dp[2] = 2
for n=3:2:N
t = n
t += (t + 1) >> 1
cnt = Int16(2)
while t >= n
if t & 1 == 0
t >>= 1
cnt += Int16(1)
else
t += (t + 1) >> 1
cnt += Int16(2)
end
end
@inbounds dp[n] = dp[t] + cnt
@inbounds dp[n+1] = dp[div(n+1, 2)] + Int16(1)
@inbounds if dp[n] > len_star
len_star = dp[n]
n_star = n
end
end
n_star, len_star
end
性能測試:
julia> @benchmark euler0014(10^8)
BenchmarkTools.Trial:
memory estimate: 190.73 MiB
allocs estimate: 2
--------------
minimum time: 1.112 s (1.06% GC)
median time: 1.123 s (1.16% GC)
mean time: 1.142 s (3.33% GC)
maximum time: 1.185 s (7.09% GC)
--------------
samples: 5
evals/sample: 1
性能總結:
運行時間總結(單位是秒):
- MATLAB: 4.134375
- Julia(沒有位運算): 1.920117
- Julia(位運算): 1.438625
- Julia(類型stable) 1.236
- Julia(@inbounds) 1.142
- C++(rsa的代碼): 2.18
- C++(我寫的代碼): 1.869
求與C++的對比, 我用的是VS2017, 速度慢, 我估計這次能真的超越C++了.
推薦閱讀：
