感謝wys和小火車普及這些技巧qwq 這篇文章大概沒什么營養
我們來看一道十分簡單的題目:
設n=131072,輸入兩個長度為n的數列$a_0,a_1...a_{n-1}$和$b_0,b_1...b_{n-1}$,要求輸出一個長度為n的數列$c_0,c_1...c_{n-1}$ 。
其中$c_i=max(a_j+b_{i-j})~(0 \leq j \leq i)$ ,$1 \leq a_i, b_i \leq 10^9$ 。
首先我們來講講這題怎么做。
如果數據是隨機的,那么有一種神奇的做法:在a和b中分別挑出最大的p個元素,對於每個i暴力枚舉每個p進行更新,這樣的復雜度是O(np)的,正確性我不會分析= =
那么數據不是隨機的...那么估計沒有什么快速的算法,不如暴力!
以下的運行時間均為在我的渣渣筆記本中測試得到,僅供參考。測試環境Ubuntu,編譯選項只有-O2。
0.
#define SZ 666666 int a[SZ],b[SZ],c[SZ]; const int n=131072; int main() { for(int i=0;i<n;i++) scanf("%d",a+i); for(int i=0;i<n;i++) scanf("%d",b+i); for(int i=0;i<=n;i++) { for(int j=0;j<=i;j++) c[i]=max(c[i],a[j]+b[i-j]); } for(int i=0;i<n;i++) printf("%d ",c[i]);puts(""); cerr<<clock()<<"ms\n"; }
simple and stupid。我們來測試一下...跑了8s。雖然不是太糟,但是還是很慢...我們來進行一波有理有據的常數優化吧。
1. 加上輸入輸出優化
const int n=131072; char ch,B[1<<20],*S=B,*T=B; #define getc() (S==T&&(T=(S=B)+fread(B,1,1<<20,stdin),S==T)?0:*S++) #define isd(c) (c>='0'&&c<='9') int aa,bb;int F(){ while(ch=getc(),!isd(ch)&&ch!='-');ch=='-'?aa=bb=0:(aa=ch-'0',bb=1); while(ch=getc(),isd(ch))aa=aa*10+ch-'0';return bb?aa:-aa; } #define gi F() #define BUFSIZE 5000000 namespace fob {char b[BUFSIZE]={},*f=b,*g=b+BUFSIZE-2;} #define pob (fwrite(fob::b,sizeof(char),fob::f-fob::b,stdout),fob::f=fob::b,0) #define pc(x) (*(fob::f++)=(x),(fob::f==fob::g)?pob:0) struct foce {~foce() {pob; fflush(stdout);}} _foce; namespace ib {char b[100];} inline void pint(int x) { if(x==0) {pc(48); return;} //if(x<0) {pc('-'); x=-x;} //如果有負數就加上 char *s=ib::b; while(x) *(++s)=x%10, x/=10; while(s!=ib::b) pc((*(s--))+48); } int main() { for(int i=0;i<n;i++) a[i]=gi; for(int i=0;i<n;i++) b[i]=gi; for(int i=0;i<=n;i++) for(int j=0;j<=i;j++) c[i]=max(c[i],a[j]+b[i-j]); for(int i=0;i<n;i++) pint(c[i]),pc(' ');pc(10); cerr<<clock()<<"ms\n"; }
雖然看起來只有10w多個數,我們還是加一波輸入輸出優化試試...
居然跑了10s。比原來還慢...這和預期不太相符啊...在windows上加了輸入輸出確實會變快,但是ubuntu下變慢了...大概輸入輸出少的時候最好還是不要加優化?
以下的測試全部基於輸入輸出優化,就假裝加了優化跑的更快好了。
2. 手寫stl
雖然這段代碼非常短,但是我們還是使用了一個stl:max。我們來常數優化一波!
for(int i=0;i<=n;i++) for(int j=0;j<=i;j++) if(a[j]+b[i-j]>c[i]) c[i]=a[j]+b[i-j];
測了測,跑了5.3s,比原來快了快一半!可喜可賀。
3. 把if改成三目?
這時候我想起了wys的教導:少用if,多用三目。
for(int i=0;i<=n;i++) for(int j=0;j<=i;j++) (a[j]+b[i-j]>c[i]) ?(c[i]=a[j]+b[i-j]):0;
這樣寫跑了6.1s,居然比if還慢?
有理有據的分析:正常情況下,if改成三目會變快的原因是因為消除了分支預測,分支預測錯誤跳轉的代價很大,而上面那段代碼預測錯誤幾率很小,所以if就比較快了。
4. 循環展開
為了寫起來方便,首先我們將b數組反序,這樣可以減少運算量,接下來把內層j循環展開。
int main() { for(register int i=0;i<n;i++) a[i]=gi; for(register int i=0;i<n;i++) b[n-i]=gi; for(register int i=0;i<n;i++) { int*r=b+n-i; for(register int j=0;j<=i;j+=8) { #define chk(a,b,c) if(a+b>c) c=a+b; #define par(p) chk(a[p],b[p],c[i]) par(j) par(j+1) par(j+2) par(j+3) par(j+4) par(j+5) par(j+6) par(j+7) } } for(register int i=0;i<n;i++) pint(c[i]),pc(' ');pc(10); cerr<<clock()<<"ms\n"; }
這樣理論上cpu可以對中間的代碼亂序執行,就是一次執行很多條,從而提高運行速度。
實測優化效果非常好,只跑了2.9s,比原來快了1倍多。
此外我還了解到openmp和cache blocking這兩種優化方法,但是對程序提速不明顯,這里就不提了,有興趣的自行度娘。
5. Intrinsic
這是真正的黑科技了= =orz小火車
#include "immintrin.h" #include "emmintrin.h" static __m256i a_m[SZ],b_m[8][SZ]; static int a[SZ],b[SZ],c[SZ]; __attribute__((target("avx2"))) inline int gmax(__m256i qwq) { int*g=(int*)&qwq,ans=0; (g[0]>ans)?(ans=g[0]):0; (g[1]>ans)?(ans=g[1]):0; (g[2]>ans)?(ans=g[2]):0; (g[3]>ans)?(ans=g[3]):0; (g[4]>ans)?(ans=g[4]):0; (g[5]>ans)?(ans=g[5]):0; (g[6]>ans)?(ans=g[6]):0; (g[7]>ans)?(ans=g[7]):0; return ans; } __attribute__((target("avx2"))) int main() { const int n=131072; memset(a,-127/3,sizeof(a)); memset(b,-127/3,sizeof(b)); for(register int i=0;i<n;i++) a[i]=gi; for(register int i=0;i<n;i++) b[n-i]=gi; for(register int i=0;i<=n+5;i+=8) a_m[i>>3]=_mm256_set_epi32 (a[i],a[i+1],a[i+2],a[i+3], a[i+4],a[i+5],a[i+6],a[i+7]); for(register int r=0;r<8;++r) for(register int i=0;i<=n+67;i+=8) b_m[r][i>>3]=_mm256_set_epi32 (b[i+r],b[i+1+r],b[i+2+r],b[i+3+r], b[i+4+r],b[i+5+r],b[i+6+r],b[i+7+r]); __m256i zero=_mm256_set_epi32(0,0,0,0,0,0,0,0); for(register int i=0,lj;i<n;i++) { __m256i*r=b_m[(n-i)&7]+((n-i)>>3), qwq=zero; lj=(i>>3); for(register int j=0;j<=lj;j+=8) { #define par(p) qwq=_mm256_max_epi32(\ qwq,_mm256_add_epi32(a_m[p],r[p])); par(j) par(j+1) par(j+2) par(j+3) par(j+4) par(j+5) par(j+6) par(j+7) } c[i]=gmax(qwq); } for(register int i=0;i<n;i++) pint(c[i]),pc(' ');pc(10); cerr<<clock()<<"ms\n"; }
這段代碼有點長,這里解釋一下原理。
大家都知道stl中有一個很好用的庫叫bitset,它的原理是將32/64個bit(取決於字長)壓成一個數(long),從而使常數/=32或64。
在intel部分指令集中,有類似的數據類型,可以將多個int/float/double等等壓在一個128/256/512位的數據類型里,從而一起進行計算。
大致有三種數據類型:
__m128i __m256i __m512i
分別對應壓在128/256/512位內。
我們可以對這三種數據類型中壓的數進行並行計算!例如一個__m256i里可以包8個int。這些數據類型的方法有點多,intel提供了一個可以查找這些方法的頁面:https://software.intel.com/sites/landingpage/IntrinsicsGuide/
實現起來相當於手寫bitset,細節詳見代碼吧。
這段代碼只跑了1.4s!又比循環展開快了一倍。
我實測了一下,在uoj上__m256i無法使用,__m128i只能使用部分指令,例如_mm_max_epi32這個指令就不支持......洛谷上可以正常運行。
upd:一些卡常方面的tips:
register基本沒用。
全局數組開static(有可能)可以放進L1~L3,明顯加快速度。
多維數組后面幾維別開2的次冪,可能導致cache miss。
inline的速度比手動展開或__attribute__((always_inline))慢。
一般什么卡常都比不上輸入輸出優化效果好= =