あと、>>505や>>540で言われていた局所性ということもやっと分かりました
最初はb~z[i],b~z[j]を使っていたので、キャッシュが追い出されて十分に使えなかったんだと思います

しかし、試しに

for( int j = i+1; j < max; j++)
　　XY[c[j].r] += c[j].b;

を消して計算したところ15秒ほどかかったので、
局所性を使ってもそれほど速度は変わらないかと思いました
>>482で95%がX[r]へのランダムアクセスと書かれていましたが
もしかしたら、どこか他にボトルネックがあるのかもしれません

また、試しに

int B = 0;
for( int i = 0; i < max; i++)
　for( int j = i+1; j < max; j++)
B++;

と、書いてみたところ、計測不可能なくらい時間がかかりました
b~z[N]を定義していないプログラムだと3秒程度で計算が終わったので
ヒープ領域が関係あるのかと思いましたが、もし原因がわかる方がいたら教えて下さい。
連投申し訳ありません。
520さんのコードでももう少し調べてみようと思います本当にありがとうございます。