随便写写
一如既往的，看到 HPCGame 又开赛了。想去玩玩并且顺便再整点奖金，于是我又参加了。只可惜今年的比赛周实在是太忙了（周内在赶大作业，周末在打另外一个 CTF），于是就只能随便做了做。不过好在 o1 足够强大，帮我拿了不少分，最后也是一如既往地拿了二等奖（不过其实今年稍微多用一点时间的话就能拿到一等了）。
本次比赛的评测环境上手难度有点高。不过成功上手后用起来还算是挺方便的。主要差评点是 pod 有时候会重启（似乎是资源耗尽之类的情况），以及服务器上硬件 profiler 跑不起来，所以只能大概从软的 profiling 记录看热点去猜测可以优化的地方。
参加完比赛之后一段时间人都挺忙的（或者虽然不忙但是不方便开电脑），好不容易闲一点的情况下也只想摸鱼。所以过了一个月才开始写这篇文章。
读完说明后不难找到答案。不过其实直接塞给 LLM 让 LLM 去找也没问题。另外经验丰富的人看到前两位应该就能猜出答案是 1898 了吧。
问一问 GPT/o1，不难得到答案：
4 8
否
128 2048
3 0.5
BCE
5.0.6 3.0.7
DABC
53124
2.0
0.5 1
当然个别题 GPT 是做不出来的。比如第四题，有不小概率算错。第六题还是得自己去翻版本列表以及 Changelog。非常好玩的是，第七题不用给 o1 图片，它也能做出来，可能因为图里面这个顺序确实太有规律了。
有了上一年的经验，我们不难看出题目的示例代码中，即便是让编译器自动向量化了，浮点流水线也没有打满。因此我们手动循环展开一下，让编译器能编译出打满流水线的向量化操作即可。
当然，具体编译选项我们还是问一下 GPT，对每个编译器都让 GPT 输出一个最好的选项。之后我们手动测一测比较一下选一选就好了。
void filter_run(double x[113][2700], double wgt[113][1799], int ngrid[113], int is, int ie, int js, int je)
{
double tmp[ie - is + 10];
for(int j = js; j <= je; j++){
int n = ngrid[j];
int hn = (n - 1) >> 1;
for(int i = is; i <= ie; i++){
tmp[i - is] = 0;
}
for(int i = is; i <= ie; i+=8)
{
for(int p = 0; p < n; p++)
{
double w= wgt[j][p];
tmp[i + 0 - is] += w * x[j][i + 0 - hn + p];
tmp[i + 1 - is] += w * x[j][i + 1 - hn + p];
tmp[i + 2 - is] += w * x[j][i + 2 - hn + p];
tmp[i + 3 - is] += w * x[j][i + 3 - hn + p];
tmp[i + 4 - is] += w * x[j][i + 4 - hn + p];
tmp[i + 5 - is] += w * x[j][i + 5 - hn + p];
tmp[i + 6 - is] += w * x[j][i + 6 - hn + p];
tmp[i + 7 - is] += w * x[j][i + 7 - hn + p];
}
}
for(int i = is; i <= ie; i++)
{
x[j][i] = tmp[i - is];
}
}
}
cmake_minimum_required(VERSION 3.20)
project(FilterProject LANGUAGES C CXX)
SET(CMAKE_C_FLAGS "-Ofast -O3 -ffast-math -march=native -funroll-loops -unroll=16 -fstrict-aliasing")
SET(CMAKE_EXE_LINKER_FLAGS "")
add_executable(program main.cpp filter.F90)
set_source_files_properties(filter.F90 PROPERTIES LANGUAGE C)
一眼：这不是我们 OI 题吗。众所周知，如果是排列求 LCS，可以对其中一个排列做值-位置的转换，之后在这个位置数组上做 的 LIS 就好了。那么本题不是排列，能不能这么做呢。当然也是可以的，但是代价就是如果出现了同一个值在另一个数组中有多个可能的位置，每个位置都需要被添加到数组中。这就意味着如果某个数字 出现了 次，那么在新的位置数组中，其对应的位置总共会出现 次。那么在本题中，由于输入是在 65536 内均匀随机的，所以有新数组的长度 ，总复杂度。我们算一下不难发现，单线程就能过题了。
另外代码也几乎不用自己写，有人写过这样的 LCS，网上能查到，改一改接口就好了。
#include <utility>
#include <cstdlib>
#include <algorithm>
#include <iostream>
typedef int element_t;
size_t lcs(element_t* arr_1, element_t* arr_2, size_t len_1, size_t len_2) {
std::vector <int> vec[0xffff+1];
size_t mlen = (len_1>len_2?len_1:len_2)+1;
long long tt = mlen;
int* num = (int*)calloc(tt*tt/50000+1000, sizeof(int));
int* tmp = (int*)calloc(mlen + 10, sizeof(int));
for(int i=0;i<len_1;i++){
vec[arr_1[i]].push_back(i);
}
int cnt=0;
for(int i=0;i<len_2;i++){
if(vec[arr_2[i]].size()!=0){
for(int j=vec[arr_2[i]].size()-1;j>=0;j--){
num[cnt++]=vec[arr_2[i]][j];
}
}
}
int counter=0;
for(int i=0;i<cnt;i++){
if(counter==0){
tmp[0]=num[0];
counter++;
}
else{
if(tmp[counter-1]<num[i]){
tmp[counter]=num[i];
counter++;
}
else{
int s=std::lower_bound(tmp,tmp+counter,num[i])-tmp;
tmp[s]=num[i];
}
}
}
free(num);
free(tmp);
return counter;
}
本题是一个 MPI 题。做法也是非常的显然：将整个棋盘按行均分成多个块分别给 n 个进程，每个进程只维护自己的当前块就好了。当然为了处理蔓延的情况，每个进程还要额外维护当前块的上下相邻行。每个 tick 后不同进程之间还需要同步一下。使用循环数组的方法压一下内存。这个题 o1 肯定是会写的，我们把题面和大概思路告诉它然后复制代码就好了（实际上我做的时候是给了一个简化的情况的，自己后面又改了改）。
#include <mpi.h>
#include <iostream>
#include <vector>
#include <cstdlib>
#include <fstream>
using namespace std;
const int TREE = 1;
const int FIRE = 2;
const int ASH = 3;
const int EMPTY = 0;
struct Event {
int ts; // 事件触发的时间步
int type; // 事件类型：1（天降惊雷），2（妙手回春）
int x1, y1; // 事件的坐标或区域范围
int x2, y2; // 仅用于“妙手回春”事件
};
const int MAXN = 16384;
const int MAXT = 2000;
int n, m, t;
char board[2][MAXN][MAXN];
char outboard[MAXN][MAXN];
int getnext(istream &fin, int cur) {
static int i = 0;
static int las = -1;
if(i>=m) return t+114514;
if(las>=cur) {
return las;
}
fin >> las;
i++;
return las;
}
int main(int argc, char **argv) {
if (argc < 3) {
std::cerr << "Usage: " << argv[0] << " <input_file> <output_file>" << std::endl;
return 1;
}
const char *input_file = argv[1];
const char *output_file = argv[2];
cout<<"A"<<endl;
MPI_Init(&argc, &argv);
int rank, size;
MPI_Comm_rank(MPI_COMM_WORLD, &rank);
MPI_Comm_size(MPI_COMM_WORLD, &size);
cout<<"IN"<<endl;
// 读取输入文件
std::ifstream fin(input_file);
if (!fin) {
std::cerr << "Error: cannot open file '" << input_file << "'." << std::endl;
MPI_Abort(MPI_COMM_WORLD, 1);
}
cout<<"C"<<endl;
fin >> n >> m >> t;
cout<<"?"<<endl;
for (int i = 0; i < n; i++) {
for (int j = 0; j < n; j++) {
int x;
fin >> x;
board[0][i][j] = x;
}
}
cout<<"B"<<endl;
int mysx = rank * n / size;
int myex = (rank + 1) * n / size;
int nesx = std::max(0, mysx - 1);
int neex = std::min(n, myex + 1);
for(int ts = 0; ts < t; ts++) {
// 1. 本地更新
// 2. 蔓延操作
// 3. 同相邻进程交换边界行，更新幽灵行
cout<<"T "<<ts<<endl;
int curb = ts&1;
if(getnext(fin,ts+1)==ts+1)
{
int op;
fin >> op;
cout<<"OP "<<ts<<' '<<op<<endl;
if(op==1)
{
int x, y;
fin >> x >> y;
if(x >= nesx && x < neex) {
if(board[curb][x][y] == TREE)
{
board[curb][x][y] = FIRE;
}
}
}
else if(op==2)
{
int x1, y1, x2, y2;
fin >> x1 >> y1 >> x2 >> y2;
for(int i = std::max(x1, nesx); i <= std::min(x2, neex-1); i++) {
for(int j = y1; j <= y2; j++) {
if(board[curb][i][j] == ASH) {
board[curb][i][j] = TREE;
}
}
}
}
}
for(int i = mysx; i < myex; i++) {
for(int j = 0; j < n; j++) {
board[!curb][i][j] = board[curb][i][j];
if(board[curb][i][j] == FIRE) {
board[!curb][i][j] = ASH;
}
else if(board[curb][i][j] == TREE) {
if(i > 0 && board[curb][i-1][j] == FIRE) {
board[!curb][i][j] = FIRE;
}
if(i < n-1 && board[curb][i+1][j] == FIRE) {
board[!curb][i][j] = FIRE;
}
if(j > 0 && board[curb][i][j-1] == FIRE) {
board[!curb][i][j] = FIRE;
}
if(j < n-1 && board[curb][i][j+1] == FIRE) {
board[!curb][i][j] = FIRE;
}
}
}
}
MPI_Request reqs[4];
MPI_Status stats[4];
int cnt = 0;
if(rank > 0) {
MPI_Isend(board[!curb][mysx], n, MPI_CHAR, rank-1, 0, MPI_COMM_WORLD, &reqs[cnt++]);
MPI_Irecv(board[!curb][nesx], n, MPI_CHAR, rank-1, 0, MPI_COMM_WORLD, &reqs[cnt++]);
}
if(rank < size-1) {
MPI_Isend(board[!curb][myex-1], n, MPI_CHAR, rank+1, 0, MPI_COMM_WORLD, &reqs[cnt++]);
MPI_Irecv(board[!curb][neex-1], n, MPI_CHAR, rank+1, 0, MPI_COMM_WORLD, &reqs[cnt++]);
}
MPI_Waitall(cnt, reqs, stats);
// for(int i=mysx; i<myex; i++) {
// for(int j=0; j<n; j++) {
// cout<<ts<<" "<<i<<" "<<j<<" "<<(int)board[!curb][i][j]<<endl;
// }
// }
}
// 收集所有进程计算结果（去掉幽灵行），放到 rank=0 进程进行输出
MPI_Gather(
&board[t&1][mysx][0], // 发送缓冲区首地址(跳过顶部幽灵行)
(myex-mysx) * 16384, // 发送元素条数
MPI_CHAR,
&outboard[0][0], // 接收缓冲区首地址(仅在 root=0 进程上使用)
(n/size) * 16384,
MPI_CHAR,
0, // root 进程
MPI_COMM_WORLD);
if(rank == 0) {
std::ofstream fout(output_file);
if (!fout) {
std::cerr << "Error: cannot open file '" << output_file << "'." << std::endl;
MPI_Abort(MPI_COMM_WORLD, 1);
}
for (int i = 0; i < n; i++) {
for (int j = 0; j < n; j++) {
fout << (int)outboard[i][j] << " ";
}
fout << std::endl;
}
fout.close();
}
MPI_Finalize();
return 0;
}
神神算法题。一眼：这不是今年 Hackergame2024 的题吗。如果没有空方块的话，本题做法就完全一样了，只需要设首行值为未知数，之后逐行递推到最后一行，再解约束就好了。整体复杂度是 的。对于空间，该方法也可以使用滚动数组（因为下一行只会依赖上面两行），因此空间可以被压到 。当我们解约束解出第一行的值后，再递推一遍就可以得到全图了。
但可惜有空方块。空方块的个数 的规模大约是 ，蛮大的。注意到通常情况下，上述做法当遇到一个空方块时，空方块的上面一格会被额外添加一条约束，而空方块的下面一格则会失去递推的约束，进而多了一个自由元。如果我们正常保留这个约束和自由元，那么最终解约束的复杂度就变成了 ，同时空间占用也非常大。
因此我们考虑使用动态消元的方法。即当添加进一个额外约束时，此时该约束中涉及的变量均是已被定义过的变量，使用这个约束把某一个变量解出来，即用其他变量表示出。同时再把这个变量回代到棋盘中。
在具体实现的时候，有一些需要注意的设计细节。注意到当一个变量被消元后，如果我们希望节省空间，则应当重新把这个变量的位置分配给其他新的变量，而这就需要我们做一些变量释放分配的方案。同时由于自由元可能出现在任意一行，如果使用滚动数组的话，则当某个变量被消元后，维护棋盘信息，以及在解出结果回代棋盘会比较复杂。为了维护简单，我没有使用滚动数组的方法，当然这会导致时间上略慢。另外一种可能的方案是，对每个自由元维护两套编号，原始编号以及压缩后的存储用编号，这样使用滚动数组时可能会比较方便。
在使用了重分配策略后，简单测一下可以观测到，同时使用的最大变量也就只有几百个，但每次代入消元时需要扫描之前的整个棋盘非常慢。我们也可以想到一个非常有趣的优化方案：如果一个自由元对应位置在第 x 行，那么这意味着它是在第 x 行扫描时才被引入的，也就意味着 x 行的棋盘不会涉及到这个变量，除非它在其他变量换元时被代入了。当我们严格使用额外约束中最新被添加的变量时，这个例外情况就不会被触发（当然我们需要维护一下最新时间了）。
之后再观测一下热点，发现在递推执行到最后一行时，每一个格子都会触发自动消元，而且消元的过程非常慢。那么显然我们参照原始的方法，不要在最后一行做代入了，直接把最后一行的约束列出来并统一做高斯消元，最后回代棋盘即可。此时程序已经跑的非常快了，快能满分了。简单加一下 omp，发现没什么用。这是因为能加的循环都是小循环，并且只有两个线程也不会有什么大的提升。
只能再观测热点了，发现热点主要在回代。在写代码时，为了安全性我写了很多 if-branch。在后续优化时发现有很多的是不需要的，被其他条件已经覆盖过的。我们把这些优化了后能快不少。同时代码中有大量的地方都对 3 取模了。如果我们把取模的次数尽可能减少，也可以优化速度。做了优化后飞快，就过了。
#include<iostream>
#include<cstdio>
#include<algorithm>
#include<queue>
#include<cstdlib>
#include<cstring>
#include<omp.h>
using namespace std;
#define DEBUG 0
#define STDINOUT 0
//free var manager
const int MAXVAR = 1024;
priority_queue<int, vector<int>, greater<int> > free_vars;
int curr_idx = 0;
int max_idx = 0;
bool var_used[MAXVAR+1];
vector<int> vars_time_order;
int get_free_var()
{
int ret;
if(free_vars.empty())
{
if(curr_idx == MAXVAR)
{
cout << "Error: No more free variables" << endl;
exit(1);
}
ret=++curr_idx;
max_idx = max(max_idx, curr_idx);
}
else
{
ret = free_vars.top();
free_vars.pop();
}
var_used[ret] = true;
vars_time_order.push_back(ret);
#if DEBUG
cout << "New Var " << ret << endl;
#endif
return ret;
}
void free_var(int idx)
{
free_vars.push(idx);
var_used[idx] = false;
int pos = -1;
for(int i=vars_time_order.size()-1; i>=0; i--)
{
if(vars_time_order[i] == idx)
{
pos = i;
break;
}
}
#if DEBUG
cout<<"Free Var "<<idx<<endl;
#endif
if(pos == -1)
{
cout << "Error: Freeing non-allocated variable" << endl;
cout << idx << endl;
exit(1);
}
vars_time_order.erase(vars_time_order.begin()+pos);
if(idx == curr_idx)
{
while(!free_vars.empty() && free_vars.top() == curr_idx)
{
free_vars.pop();
--curr_idx;
}
}
}
int var_x[MAXVAR+1], var_y[MAXVAR+1];
int n,m;
int smat[1024][512];
int amat[1024][512];//answer matrix
// char cureq[4][512][MAXVAR+1];
// bool haseq[4][512];
char cureq[1024][512][MAXVAR+1];
bool haseq[1024][512];
int fxx[]={0,0,0,-1,1};
int fyy[]={0,-1,1,0,0};
char tmpeq[MAXVAR+1];
int eqs[MAXVAR+1][MAXVAR+1];
int eqcnt=0;
int ans[MAXVAR+1];
int main()
{
omp_set_num_threads(2);
#if STDINOUT
cin>>m>>n;
for(int i=0; i<n; i++)
{
for(int j=0; j<m; j++)
{
cin>>smat[i][j];
}
}
#else
//read binary file
FILE *fp = fopen("in.data", "rb");
if(fp == NULL)
{
cout << "Error: File not found" << endl;
return 1;
}
//read the number of rows and columns
fread(&m, sizeof(int), 1, fp);
fread(&n, sizeof(int), 1, fp);
//read the matrix
for(int i=0; i<n; i++)
{
fread(smat[i], sizeof(int), m, fp);
}
#endif
var_used[0] = true;
for(int i=0; i<n; i++)
{
// int curl=i&3;
// memset(haseq[nexl],0,sizeof(haseq[nexl]));
// memset(cureq[nexl],0,sizeof(cureq[nexl]));
for(int j=0; j<m; j++)
{
if(smat[i][j]==0)
continue;
// cout<<"SS "<<i <<' '<<j<<' '<<smat[i][j]<<endl;
int freecnt=0;
int freek=-1;
int lasvar=-1;
for(int k=0; k<5; k++)
{
int nx=i+fxx[k];
int ny=j+fyy[k];
if(nx<0 || nx>=n || ny<0 || ny>=m)
continue;
if(smat[nx][ny]==0)
continue;
if(!haseq[nx][ny])
{
freecnt++;
freek=k;
}
lasvar=k;
}
if(freecnt==0)//do gauss
{
if(lasvar==-1)
{
cout << "Error: No variable" << endl;
return 1;
}
memset(tmpeq,0,sizeof(tmpeq));
tmpeq[0]=(smat[i][j])%3;
for(int k=0; k<5; k++)
{
int nx=i+fxx[k];
int ny=j+fyy[k];
if(nx<0 || nx>=n || ny<0 || ny>=m)
continue;
if(smat[nx][ny]==0)
continue;
// #pragma omp parallel for
for(int l=0; l<=curr_idx; l++)
{
if(cureq[nx][ny][l]==0)
continue;
tmpeq[l]=(tmpeq[l]+cureq[nx][ny][l]);
}
}
for(int l=0;l<=curr_idx;l++)
{
tmpeq[l]=(tmpeq[l]+3*16384)%3;
}
//when not last line
bool addeq=false;
if(i<n-1)
{
int elmvar=-1;
for(int ll=vars_time_order.size()-1; ll>=0; ll--)
{
int l=vars_time_order[ll];
if(tmpeq[l]!=0)
{
elmvar=l;
break;
}
}
// for(int l=curr_idx; l>=1; l--)
// {
// if(tmpeq[l]!=0)
// {
// elmvar=l;
// break;
// }
// }
if(elmvar<=0)
{
cout << "Error: No solution" << endl;
return 1;
}
if(tmpeq[elmvar]!=1)
{
int inv=tmpeq[elmvar];
for(int l=0; l<=curr_idx; l++)
{
tmpeq[l]=(tmpeq[l]*inv)%3;
}
}
int elmx=var_x[elmvar];
int stx=elmx-2;
if(stx<0)
stx=0;
// cout<<elmvar<<' '<<var_x[elmvar]<<' '<<var_y[elmvar]<<endl;
// int ddcnt=0;
for(int di=stx;di<n;di++)
{
if(di>i+1)
break;
// #pragma omp parallel for
for(int dj=0;dj<m;dj++)
{
// if(smat[di][dj]==0)
// continue;
if(!haseq[di][dj])
continue;
int inv=cureq[di][dj][elmvar];
if(inv==0)
continue;
// cout<<i<<' '<<j<<' '<<di<<' '<<dj<<' '<<inv<<endl;
for(int l=0; l<=curr_idx; l++)
{
cureq[di][dj][l]=(cureq[di][dj][l]-tmpeq[l]*inv+9)%3;
}
// ddcnt++;
}
}
// if(i>=511)
// cout<<"DDCNT "<<i<<' '<<j<<' '<<ddcnt<<endl;
if(!addeq)
{
free_var(elmvar);
}
#if DEBUG
cout << "Equation " << eqcnt <<" By " << i << " " << j << endl;
for(int l=0; l<=curr_idx; l++)
{
cout << int(tmpeq[l]) << " ";
}
cout << endl;
cout << "Free Var " << elmvar << endl;
#endif
}
if(i==n-1 || addeq)
{
for(int l=0; l<=curr_idx; l++)
{
eqs[eqcnt][l]=tmpeq[l];
}
eqcnt++;
#if DEBUG
cout << "Equation " << eqcnt <<" By " << i << " " << j << endl;
for(int l=0; l<=curr_idx; l++)
{
cout << int(tmpeq[l]) << " ";
}
cout << endl;
#endif
}
}
else//create equation
{
for(int k=0; k<5; k++)
{
int nx=i+fxx[k];
int ny=j+fyy[k];
if(nx<0 || nx>=n || ny<0 || ny>=m)
continue;
if(smat[nx][ny]==0)
continue;
if(!haseq[nx][ny])
{
if(freek!=k)
{
int var=get_free_var();
cureq[nx][ny][var]=1;
haseq[nx][ny]=1;
var_x[var]=nx;
var_y[var]=ny;
#if DEBUG
cout<<"Var "<<var<<" = "<<nx<<" "<<ny<<endl;
#endif
}
}
}
int fkx=i+fxx[freek];
int fky=j+fyy[freek];
haseq[fkx][fky]=1;
cureq[fkx][fky][0]=(3-smat[i][j])%3;
for(int k=0; k<5; k++)
{
int nx=i+fxx[k];
int ny=j+fyy[k];
if(nx<0 || nx>=n || ny<0 || ny>=m)
continue;
if(smat[nx][ny]==0)
continue;
if(k==freek)
continue;
// #pragma omp parallel for
for(int l=0; l<=curr_idx; l++)
{
if(cureq[nx][ny][l]==0)
continue;
cureq[fkx][fky][l]=(cureq[fkx][fky][l]-cureq[nx][ny][l]+3)%3;
}
}
#if DEBUG
cout << "Expression " << fkx << " " << fky << " By " << i << " " << j << endl;
for(int l=0; l<=curr_idx; l++)
{
cout << int(cureq[fkx][fky][l]) << " ";
}
cout << endl;
#endif
}
}
}
//gauss
cout<<"Gauss "<<eqcnt<<' '<<vars_time_order.size()<<endl;
for(int i=0; i<eqcnt; i++)
{
int k=i;
int cv=vars_time_order[i];
eqs[k][cv]=(eqs[k][cv]+3)%3;
if(eqs[k][cv]==0)
{
for(int j=i+1; j<eqcnt; j++)
{
eqs[j][cv]=(eqs[j][cv]+3)%3;
if(eqs[j][cv]!=0)
{
k=j;
break;
}
}
}
if(k!=i)
{
for(int j=0; j<=curr_idx; j++)
{
swap(eqs[i][j],eqs[k][j]);
}
}
for(int j=0;j<=curr_idx;j++)
{
eqs[i][j]=(eqs[i][j]%3+3)%3;
}
if(eqs[i][cv]==0)
{
cout << "Error: No solution" << endl;
return 1;
}
if(eqs[i][cv]!=1)
{
int inv=eqs[i][cv];
for(int j=0; j<=curr_idx; j++)
{
if(!var_used[j])
continue;
eqs[i][j]=(eqs[i][j]*inv)%3;
}
}
for(int j=0; j<eqcnt; j++)
{
if(j==i)
continue;
eqs[j][cv]=(eqs[j][cv]+3)%3;
if(eqs[j][cv]==0)
continue;
int inv=eqs[j][cv];
for(int l=0; l<=curr_idx; l++)
{
if(!var_used[l])
continue;
eqs[j][l]=(eqs[j][l]-eqs[i][l]*inv);
}
}
}
//get answer
for(int i=0; i<eqcnt; i++)
{
int cv=vars_time_order[i];
ans[cv]=(3-eqs[i][0]%3)%3;
#if DEBUG
cout << "Ans x" << i+1 << " = " << ans[i+1] << " At " << var_x[i+1] << " " << var_y[i+1] << endl;
#endif
}
//write answer matrix
for(int i=0; i<n; i++)
{
for(int j=0; j<m; j++)
{
int a=0;
// #pragma omp parallel for reduction(+:a)
for(int k=1; k<=curr_idx; k++)
{
a=(a+ans[k]*cureq[i][j][k]);
}
a=(a+cureq[i][j][0])%3;
amat[i][j]=(a+3)%3;
}
}
#if STDINOUT
for(int i=0; i<n; i++)
{
for(int j=0; j<m; j++)
{
cout << amat[i][j] << " ";
}
cout << endl;
}
#else
//write binary file
FILE *fp2 = fopen("out.data", "wb");
//write the matrix
for(int i=0; i<n; i++)
{
fwrite(amat[i], sizeof(int), m, fp2);
}
#endif
cout<<max_idx<<' '<<eqcnt<<endl;
return 0;
}
怎么是个 julia 题，但是我不会 julia 啊。不过我知道 GPT 肯定会。首先本题的做法非常显然：把数组按线程均分，之后每个线程内求 topK ，之后不同线程的归并起来再一下序就好了。
我们直接让 o1 写代码。o1 告诉我用堆排非常快，但是我测了下发现它写的非常慢，怀疑是堆排根本不行，于是我让他只用官方库的 kth，就调理好了。事后发现原来是 o1 写的烂，不是堆排的问题。
之后我还尝试了一些奇技淫巧。比如说我们知道数据是均匀分布的情况下，可以直接筛掉很多小的数据，只在比较大的数据内求 topK。可惜我照猫画虎是写不出高性能 julia 的，写出来的都更慢了，于是只能作罢。不过本题还是拿了 98 分。
using Base.Threads
mutable struct Atomic{T} # 题目自带的原子类型定义，可无视
@atomic x::T
end
# 并行分块地对 data 做部分排序，保留最靠后的 K 个元素(即全局 top-k)
# 然后把所有线程的 top-k 拼成一个数组，再做一次 partialsort!，
# 最后返回按题目要求顺序(值大的优先；值相等时索引小的优先)排好的下标。
@inbounds function topk(data::AbstractVector{T}, k::Int) where T
n = length(data)
nt = nthreads()
chunk_size = cld(n, nt)
partial_idx_list = Vector{Vector{Int}}(undef, nt)
@threads for t in 1:nt
start_i = (t-1)*chunk_size + 1
end_i = min(t*chunk_size, n)
if start_i <= end_i
# 用局部的下标数组
local_idxs = collect(start_i:end_i)
# 定义 lt，使得 partialsort! 会将“值大”的元素放到数组 local_idxs 的末尾
# partialsort!(A, range, lt=...) 会把“最小”(由 lt 判定) 排到 A 的前面，
# 因此要让值大/索引小 的元素视作“更大”，就需要 lt(i,j)=true 表示 i“更小”。
# 因此 lt(i,j) = data[i]<data[j] 或 相等时 i>j
# 这样就能确保 top-k 在 local_idxs[end-k+1:end]。
partialsort!(local_idxs, (end_i - start_i + 2 - k):(end_i - start_i + 1),
lt = (i, j) -> (data[i] < data[j]) || (data[i] == data[j] && i > j)
)
# 取出末尾 k 个(它们是本分块里“最大的 k 个”)
partial_idx_list[t] = last(local_idxs, k)
else
partial_idx_list[t] = Int[]
end
end
# 收集所有线程各自的 top-k
merged_idxs = reduce(vcat, partial_idx_list)
# 在 merged_idxs 里只需再次求全局 top-k。
# 同理，用 partialsort! 并把 top-k 放到 merged_idxs[-k+1:end]。
partialsort!(merged_idxs, (length(merged_idxs)-k+1):length(merged_idxs),
lt = (i, j) -> (data[i] < data[j]) || (data[i] == data[j] && i > j)
)
# 现在 merged_idxs[-k+1:end] 就是全局前 k 大元素的下标，按照“从小到大”放在后面。
# 接下来按题目要求(值大的优先，若值相等则索引小优先)对这 k 个下标排好序输出。
topk_idxs = last(merged_idxs, k)
# 逆向排序：值大排前，若值相等则索引小排前
sort!(topk_idxs; lt = (i,j) -> begin
if data[i] == data[j]
return i < j
else
return data[i] > data[j]
end
end)
return topk_idxs
end
没做。看着就不是能短时间内拿到分的题。不过当时忘了去试试 o1 能写出什么样的代码了，毕竟我连题目都懒得看。
我直接把参考代码喂给 o1，让 o1 给我写了个 cuda。o1 写的初版最开始不太满意，我做了一番简单指导，让它自己想想并修改，比如做些预处理什么的，它就写了一个更快的代码，拿了 27 分。
之后尝试让 o1 做一些更厉害的优化，发现实际上跑的都更慢了。或者用混合精度计算，但最终结果的精度很差，过不了题。我也尝试让 o1 做误差分析（正解），但 o1 做的不太行。不过答案的误差分析确实还是有一定难度的，我估计没一定时间也做不出来。
#include <cuda_runtime.h>
#include <math.h>
#include <stdio.h>
#include <stdint.h>
#include <stdlib.h>
// 与原先相同，用于文件读写(double)
typedef double d_t;
struct d3_t {
d_t x, y, z;
};
// 内联距离函数(双精度)
__host__ __device__ inline double dist_d3(const d3_t& a, const d3_t& b)
{
double dx = a.x - b.x;
double dy = a.y - b.y;
double dz = a.z - b.z;
return sqrt(dx*dx + dy*dy + dz*dz);
}
// CUDA 的内置函数 sincos，一次性计算 sin 和 cos
__device__ inline void dsincos(double x, double &s, double &c)
{
sincos(x, &s, &c);
}
//------------------------------------------------------------------------------
// 预计算 kernel：mirror->src 的距离
// 不做越界检查，假定 blockSize×gridSize = mirn
//------------------------------------------------------------------------------
__global__ void kernelPreComputeSrcDist(const d3_t* __restrict__ d_mir,
d3_t src,
d_t* __restrict__ d_Lsrc)
{
// 全局线程 ID
int64_t tid = blockIdx.x * (int64_t)blockDim.x + threadIdx.x;
// 假设 tid < mirn
d_Lsrc[tid] = dist_d3(d_mir[tid], src);
}
//------------------------------------------------------------------------------
// 主核函数：分块(tiled)；不用越界检查，假设 mirn%tileSize=0
// 每个线程负责一个 sensor，下标 i = blockIdx.x*blockDim.x+threadIdx.x
//------------------------------------------------------------------------------
__global__ void kernelComputeDataTiled(
const d3_t* __restrict__ d_mir, // [mirn]
const d_t* __restrict__ d_Lsrc, // [mirn]
const d3_t* __restrict__ d_sen, // [senn]
d_t* __restrict__ d_out, // [senn]
int64_t mirn, // 例如 1048576
int64_t tileSize // 例如 1024
) // 2π*2000
{
// 每个线程对应一个 sensor i
const double C = 6.283185307179586 * 2000.0;
int64_t i = (int64_t)blockIdx.x * blockDim.x + threadIdx.x;
// 读出该传感器坐标
d3_t s = d_sen[i];
// 寄存器累加器
double accA = 0.0;
double accB = 0.0;
// 共享内存: 缓存本 tile 的 mirror 坐标+mirror->src 距离
// mirn, tileSize 都是 2^n，可假设 tileSize=1024
__shared__ d3_t shMir[1024];
__shared__ d_t shL1 [1024];
// tileCount = mirn/tileSize
// for(...) 没越界检查
for (int64_t base = 0; base < mirn; base += tileSize) {
// 将 tileSize 个 mirror 坐标和 mirror->src 距离搬到共享内存
// 每个 block 有 blockDim.x=256 线程，保守循环：t += blockDim.x
for (int t = threadIdx.x; t < tileSize; t += blockDim.x) {
shMir[t] = d_mir [base + t];
shL1[t] = d_Lsrc[base + t];
}
__syncthreads();
// 对 tileSize=1024 个 mirror 做累加
// 加上 #pragma unroll，可让编译器展开
#pragma unroll 8
for (int k = 0; k < 1024; k++) {
double dist = shL1[k] + dist_d3(shMir[k], s);
double val;
val = C * dist; // 频率相位
double ss, cc;
dsincos(val, ss, cc);
accA += cc;
accB += ss;
}
__syncthreads();
}
// 写回结果
d_out[i] = sqrt(accA*accA + accB*accB);
}
//------------------------------------------------------------------------------
// 主函数：假设 mirn, senn 均为 2^n，且 blockSize×gridSize= senn, tileSize 整除 mirn
//------------------------------------------------------------------------------
int main()
{
//================= 1. 读入 in.data (double) =================
FILE* fi = fopen("in.data","rb");
if (!fi) {
fprintf(stderr, "Error: cannot open in.data\n");
return 1;
}
d3_t src;
fread(&src, 1, sizeof(d3_t), fi);
int64_t mirn;
fread(&mirn, 1, sizeof(int64_t), fi);
d3_t* mir = (d3_t*)malloc(mirn * sizeof(d3_t));
fread(mir, 1, mirn * sizeof(d3_t), fi);
int64_t senn;
fread(&senn, 1, sizeof(int64_t), fi);
d3_t* sen = (d3_t*)malloc(senn * sizeof(d3_t));
fread(sen, 1, senn * sizeof(d3_t), fi);
fclose(fi);
//================= 2. 分配输出数组 (double) =================
d_t* data = (d_t*)malloc(senn * sizeof(d_t));
//================= 3. 在 GPU 上分配内存 =================
d3_t *d_mir = nullptr, *d_sen = nullptr;
d_t *d_Lsrc = nullptr, *d_out = nullptr;
cudaMalloc(&d_mir, mirn * sizeof(d3_t));
cudaMalloc(&d_sen, senn * sizeof(d3_t));
cudaMalloc(&d_Lsrc, mirn * sizeof(d_t));
cudaMalloc(&d_out, senn * sizeof(d_t));
//================= 4. 把数据拷贝到 GPU =================
cudaMemcpy(d_mir, mir, mirn * sizeof(d3_t), cudaMemcpyHostToDevice);
cudaMemcpy(d_sen, sen, senn * sizeof(d3_t), cudaMemcpyHostToDevice);
free(mir);
free(sen);
//================= 5. 配置参数并调用预计算 kernel =================
// mirn=1048576；用 blockSize=256 => gridSize=4096，无须越界判断
int blockSize = 256;
int gridSize_mir = (int)(mirn / blockSize);
kernelPreComputeSrcDist<<<gridSize_mir, blockSize>>>(d_mir, src, d_Lsrc);
cudaDeviceSynchronize();
//================= 6. 调用主核函数 (分块/tiled) =================
int gridSize_sen = (int)(senn / blockSize);
// tileSize 也可设置 512、2048，但必须整除 mirn
int tileSize = 1024;
// 2π × 2000
kernelComputeDataTiled<<<gridSize_sen, blockSize>>>(d_mir, d_Lsrc, d_sen, d_out, mirn, tileSize);
cudaDeviceSynchronize();
//================= 7. 拷回结果并写出 out.data =================
cudaMemcpy(data, d_out, senn * sizeof(d_t), cudaMemcpyDeviceToHost);
FILE* fo = fopen("out.data", "wb");
if (fo) {
fwrite(data, sizeof(d_t), senn, fo);
fclose(fo);
} else {
fprintf(stderr, "Warning: cannot open out.data for writing.\n");
}
//================= 8. 释放资源 =================
cudaFree(d_mir);
cudaFree(d_sen);
cudaFree(d_Lsrc);
cudaFree(d_out);
free(data);
return 0;
}
题目看起来是让优化一个 gemm 算法。但看了题面后发现题面疯狂暗示不用自己写。我们查一下题目给的文档发现 KML 里面就有能做 gemm 的函数。这个函数我们可以相信他是足够优化的。我们先自己跑一下 profiler，确认了热点在 sgemm 和 strsm 两个函数。我们问一下 o1，发现 o1 居然懂 KML 的用法。于是我们直接把函数体发给它，让它改成调 KML 库的即可。
#include "hpl-ai.h"
#include "kblas.h"
void sgemm(char transa, char transb, int m, int n, int k,
float alpha, float *A, int lda, float *B, int ldb,
float beta, float *C, int ldc) {
int i, j, l;
// Only supprt transa=='N', trabsb=='N'
if (transa != 'N' || transb != 'N') {
printf("Not supported in SGEMM.\n");
return;
}
if (m == 0 || n == 0) {
return;
}
if ((alpha == 0.0 || k == 0) && beta == 1.0) {
return;
}
cblas_sgemm(
CblasColMajor, // HPC Kit 中指定列优先
CblasNoTrans, // 不转置 A
CblasNoTrans, // 不转置 B
m, // M
n, // N
k, // K
alpha, // alpha
A, // A
lda, // lda
B, // B
ldb, // ldb
beta, // beta
C, // C
ldc // ldc
);
return;
}
void strsm(char side, char uplo, char transa, char diag, int m, int n,
float alpha, float *A, int lda, float *B, int ldb) {
int i, j, k;
// Only support side=='L', transa=='N', alpha==1.0.
if (side != 'L' || transa != 'N' || alpha != 1.0) {
printf("Not supported in STRSM.\n");
return;
}
if (m == 0 || n == 0) {
return;
}
CBLAS_SIDE cblas_side = CblasLeft; // 因为只支持 side=='L'
CBLAS_UPLO cblas_uplo = (uplo == 'U') ? CblasUpper : CblasLower;
CBLAS_TRANSPOSE cblas_transA = CblasNoTrans; // 因为只支持 transa=='N'
CBLAS_DIAG cblas_diag = (diag == 'N') ? CblasNonUnit : CblasUnit;
// 调用 KML 提供的 cblas_strsm，使用列优先(CblasColMajor)
cblas_strsm(
CblasColMajor, // 列优先布局
cblas_side,
cblas_uplo,
cblas_transA,
cblas_diag,
m,
n,
alpha,
A, lda,
B, ldb
);
return;
}
当然我们还需要改一下 Makefile。可以比较一下用哪种模式的并行更快，发现果然是 pthread 好。我不太会用 module，于是写了个大力的 Makefile。
CC = gcc
CFLAGS = -O3 -march=native -I/opt/huawei/HPCKit/24.12.30/kml/bisheng/include -L/opt/huawei/HPCKit/24.12.30/kml/bisheng/lib/sve/kblas/pthread -lkblas
LDFLAGS = -I/opt/huawei/HPCKit/24.12.30/kml/bisheng/include -L/opt/huawei/HPCKit/24.12.30/kml/bisheng/lib/sve/kblas/pthread -lkblas
.PHONY: all clean
all: hpl-ai
# hpl-ai.o 是预编译的
hpl-ai: blas.o gmres.o sgetrf_nopiv.o
$(CC) $(CFLAGS) hpl-ai.o $^ -o $@ -lm $(LDFLAGS)
# Pattern rule
%.o: %.c
$(CC) $(CFLAGS) -c $< -o $@
clean:
rm -f hpl-ai blas.o gmres.o sgetrf_nopiv.o
本题要优化的是一个暴力枚举私钥的代码。我们看示例代码，发现出题人真的是费尽心思劣化性能了。我们把它先改成一个正常人类会写出来的代码，比如先枚举私钥，之后对公钥匹配完整的列表，看有没有匹配上的，当列表中所有的都匹配上时退出枚举；为了避免随机数浪费时间，枚举私钥时才用递增而不是重新随机的方式。
之后我们加上 openmp 并行就好了。注意到不同线程匹配到私钥时需要对全局的匹配状态做维护，因此我们最好使用一个 critical 来保证互锁。注意到只有公钥在被匹配上时才会进入边界，而这个是很少出现的，因此几乎不会带来效率损失。
#include <iostream>
#include <iomanip>
#include <sstream>
#include <cstring>
#include <fstream>
#include <openssl/sha.h>
#include <openssl/evp.h>
#include <secp256k1.h>
#include <omp.h>
std::string toHex(const uint8_t* data, size_t size) {
std::ostringstream oss;
for (size_t i = 0; i < size; ++i) {
oss << std::hex << std::setfill('0') << std::setw(2) << (int)data[i];
}
return oss.str();
}
std::string sha3256(const uint8_t* data, size_t size) {
EVP_MD_CTX* context = EVP_MD_CTX_new();
const EVP_MD* md = EVP_get_digestbyname("sha3-256");
EVP_DigestInit_ex(context, md, nullptr);
EVP_DigestUpdate(context, data, size);
uint8_t hash[EVP_MAX_MD_SIZE];
unsigned int hashLen;
EVP_DigestFinal_ex(context, hash, &hashLen);
EVP_MD_CTX_free(context);
return toHex(hash, hashLen);
}
void generateRandomPrivateKey(uint8_t privateKey[32]) {
FILE* urandom = fopen("/dev/urandom", "rb");
int res = fread(privateKey, 1, 32, urandom);
if (res != 32) {
std::cerr << "Failed to read random data" << std::endl;
exit(1);
}
fclose(urandom);
}
void inc(uint8_t privateKey[32]) {
for (int i = 31; i >= 0; --i) {
if (privateKey[i] == 0xff) {
privateKey[i] = 0;
} else {
++privateKey[i];
break;
}
}
}
std::string computeEthereumAddress(const secp256k1_context* ctx, const uint8_t privateKey[32]) {
secp256k1_pubkey pubkey;
secp256k1_ec_pubkey_create(ctx, &pubkey, privateKey);
uint8_t pubkeySerialized[65];
size_t pubkeySerializedLen = 65;
secp256k1_ec_pubkey_serialize(ctx, pubkeySerialized, &pubkeySerializedLen, &pubkey, SECP256K1_EC_UNCOMPRESSED);
std::string hash = sha3256(pubkeySerialized + 1, pubkeySerializedLen - 1);
return hash.substr(24);
}
std::string prefixs[10];
bool found[10];
uint8_t privateKeys[10][32];
std::string addresses[10];
int main(int argc, char* argv[]) {
std::ifstream infile("vanity.in");
std::ofstream outfile("vanity.out");
for(int i = 0; i < 10; ++i){
infile >> prefixs[i];
}
int foundcnt = 0;
omp_set_num_threads(8);
#pragma omp parallel
{
uint8_t privateKey[32];
generateRandomPrivateKey(privateKey);
secp256k1_context* ctx = secp256k1_context_create(SECP256K1_CONTEXT_SIGN);
std::string address;
while (foundcnt < 10) {
inc(privateKey);
address = computeEthereumAddress(ctx, privateKey);
for(int i = 0; i < 10; ++i){
if(address.substr(0, prefixs[i].size()) == prefixs[i]){
#pragma omp critical
{
if(!found[i]){
found[i] = true;
memcpy(privateKeys[i], privateKey, 32);
addresses[i] = "0x"+address;
foundcnt++;
std::cout << "Found " << addresses[i] << " with private key " << toHex(privateKeys[i], 32) << std::endl;
std::cout << "Got " << foundcnt << " private keys" << std::endl;
}
}
}
}
}
secp256k1_context_destroy(ctx);
}
for(int i = 0; i < 10; ++i){
outfile << addresses[i] << std::endl;
outfile << toHex(privateKeys[i], 32) << std::endl;
}
return 0;
}
本来根本不想做的，后来发现有一车人都拿了十来分。于是尝试交了一下样例，发现也拿了 10 分。简单改了一下编译选项，看看能不能多拿点分，果然多拿了 1 分，变成了 11 分。
日期： 2025-02-23