単純にキャッシュの使い分けでいいんじゃね？
struct CACHE {
int r;
int b;
};
CACHE c[N];
を追加して

for(int i = 0; i<40000; i++) {
　for(int j=i+1; j<40000; j++) {
　　// 読み込み:キャッシュをフルに使う
　　// 書き込み:追い出されてもＯＫ
　　int r = static_cast<int>(sqrt(
　　　(x[i]-x[j])*(x[i]-x[j]) +
　　　(y[i]-y[j])*(y[i]-y[j]) +
　　　(z[i]-z[j])*(z[i]-z[j]))+0.5);
　　c[j].r = r;
　　c[j].b = b[i]*b[j];
　}

　for(int j=i+1; j<40000; j++) {
　　// 読み込み:追い出されてもＯＫ
　　// 書き込み:キャッシュをフルに使う
　　XY[c[j].r] += c[j].b;
　}
}

少しは速くなるはず