我突然想比较一下C语言和Julia的速度。于是在AI的辅助下写了这段代码。
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#define PI 3.14159265358979323846
#define restenergy 0.511e6
#define h1 756
#define h 3
#define circum 1360.4
#define alpha_c 1.56e-5
#define centerenergy 6e9
#define v1 3.6395053705870048e+06
#define r 0.1801760487622639
#define phi_s 2.1575815005097385
#define phi_2s 6.0620746573056303
#define nturns 10001
#define num_particles 10000
void calculate_next_values(double *zvec, double *delta_vec)
{
for (int i = 0; i < num_particles; i++) {
double z = zvec[i];
double delta = delta_vec[i];
delta += v1 / centerenergy * (sin(phi_s - (2 * h1 * PI * z) / circum) - sin(phi_s) + r * sin(phi_2s - (2 * h1 * h * PI * z) / circum) - r * sin(phi_2s));
zvec[i] = z - (alpha_c - 1 / (pow((centerenergy / restenergy), 2) * pow((1 + delta), 2))) * delta * circum;
delta_vec[i] = delta;
}
}
int main()
{
// 声明并输入初始zvec, delta_vec值
double zvec[num_particles];
double delta_vec[num_particles];
clock_t start_time, end_time;
double cpu_time_used;
// 读取CSV文件
FILE *file = fopen("data.csv", "r");
if (file == NULL) {
printf("无法打开文件！\n");
return 1;
}
char line[100];
for (int i = 0; i < num_particles; i++) {
if (fgets(line, sizeof(line), file) == NULL) {
printf("文件格式错误！\n");
return 1;
}
// 解析CSV文件中的数据
char *token = strtok(line, ",");
zvec[i] = atof(token);
token = strtok(NULL, ",");
delta_vec[i] = atof(token);
}
fclose(file);
// 记录开始时间
start_time = clock();
// 计算 z_{n+1} 和 delta_{n+1}
for(int i =0; i < nturns; i++){
calculate_next_values(zvec, delta_vec);
}
// 记录结束时间
end_time = clock();
// 计算代码片段的执行时间
cpu_time_used = ((double) (end_time - start_time)) / CLOCKS_PER_SEC;
// 输出执行时间
printf("代码片段执行时间：%lf 秒\n", cpu_time_used);
// 输出结果
printf("更新后的坐标矢量 zvec 和 delta_vec：\n");
for (int i = 0; i < 10; i++) {
printf("zvec[%d] = %lf\n", i, zvec[i]);
printf("delta_vec[%d] = %lf\n", i, delta_vec[i]);
}
return 0;
}
这段C代码，我只计算了calculate_next_values——本质上就是对矢量逐个更新——重复调用10001次的时间。结果时间花了：9.753000秒！
using BSON: @load
using BenchmarkTools
@load "zvec_δvec.bson" zvec δvec
const restenergy = 0.511e6 # [eV]
nturns = 10001
struct Par
v1::Float64
r::Float64
ϕs::Float64
ϕ2s::Float64
h1::Int64
h::Int64
circum::Float64
centerenergy::Float64
αc::Float64
end
Base.broadcastable(x::Par) = Ref(x)
function revolution(z, δ, p::Par)::Tuple{Float64,Float64}
δ = δ + p.v1 / p.centerenergy * (sin(p.ϕs - 2 * p.h1 * π / p.circum * z ) - sin(p.ϕs) + p.r * sin(p.ϕ2s - 2 * p.h1 * p.h * π / p.circum * z) - p.r * sin(p.ϕ2s))
η = p.αc - 1 / (p.centerenergy * (1 + δ) / restenergy)^2
z = z - η * δ * p.circum
z, δ
end
function main1!(zvec::Vector{Float64}, δvec::Vector{Float64}, nturns, p::Par)
npart = size(zvec, 1)
t1 = time()
for j = 1:nturns
for i = 1:npart
zvec[i], δvec[i] = revolution(zvec[i], δvec[i], p)
end
end
t2 = time()
println(t2-t1, "s")
end
p = Par(3.6395053705870048e+06, 0.1801760487622639, 2.1575815005097385, 6.0620746573056303, 756, 3, 1360.4, 6e9, 0.0000156)
main1!(zvec, δvec, nturns, p)
这里，main1和C里面的差不多，优先对矢量逐个更新，然后重复10001次，花费5.161999940872192秒。
Julia这里，我怎么还没开始优化，就比C语言快了呢？是C语言有优化空间吗？
如果有需要跑的，zvec, δvec用正态分布*1e-3就好了。