この場合, データがランダムだったとしても, その時点での最大値との比較になるので後半ほど最大値を更新しない確率が高い. つまり分岐予測があたりやすくなると思われる.
/* cmov版のコア */
cmp(eax, ptr [edx + ecx * 4]);
cmovl(eax, ptr [edx + ecx * 4]);
/* jmp版のコア */
cmp(eax, ptr [edx + ecx * 4]);
jge(".skip");
mov(eax, ptr [edx + ecx * 4]);
L(".skip");
| 乱数 | 最初が最大値 | 単調増加 | だいたい単調増加 | ||
|---|---|---|---|---|---|
| Core2Duo + Xp | jmp | 2.04 | 2.01 | 2.02 | 4.02 |
| cmov | 5.10 | 5.10 | 5.11 | 5.10 | |
| Core i7 + Win7 | jmp | 1.77 | 1.84 | 1.85 | 3.54 |
| cmov | 2.68 | 2.78 | 2.66 | 2.58 | |
| Opteron 2378 + Win2008 | jmp | 5.06 | 5.05 | 3.03 | 4.05 |
| cmov | 2.02 | 2.02 | 2.02 | 2.02 |
cmovを使うとデータによらず実行時間は一定. jmpは分岐予測のあたり外れにより時間が変動する. 最初に最大値があったり(=殆ど分岐する), 単調増加(=殆ど分岐しない)だったりするとIntel系では速い.
Opteronは分岐すると常にコストがかかる模様. ということは分岐の確率が偏っている場合, このCPUでは確率が高い方を分岐しないように配置したほうがいい(昔のIntel CPUもそうだった). cmovはかなり速い. 最近のAMD系CPUだとどうだろう. 知りたいけど手近なマシンが無かった.
データが"だいたい単調増加"のときは, 最大値の更新がときどき起こり, しかも最後まであるので分岐予測が難しい.
したがって分岐予測が外れる率も高くなり, jmp版のコストが高い.
case 2. 二つの配列を比較して左が大きいときの個数を数える
int countMax(const int *a, const int *b, size_t n)
{
int ret = 0;
for (size_t i = 0; i < n; i++) {
if (a[i] > b[i]) ret++;
}
return ret;
}
0. 条件分岐を使う, 1. setg + movzxを使う, 2. carryを使ったややトリッキーなコード, 3. setgで余計にレジスタを使う, で試してみた.
/* 0. 条件分岐 */
mov(ebx, ptr [esi + ecx * 4]); // a
cmp(ebx, ptr [edx + ecx * 4]); // b
jle(".skip");
add(eax, 1);
L(".skip");
/* 1. setg + movzx */
mov(ebx, ptr [esi + ecx * 4]);
cmp(ebx, ptr [edx + ecx * 4]);
setg(bl);
movzx(ebx, bl);
add(eax, ebx);
/* 2. carryを使う */
mov(ebx, ptr [esi + ecx * 4]);
cmp(ebx, ptr [edx + ecx * 4]);
sbb(eax, -1);
/* 3. setg + レジスタ */
mov(edi, ptr [esi + ecx * 4]);
cmp(edi, ptr [edx + ecx * 4]);
setg(bl);
add(eax, ebx);
| 割合[%] | 0 | 25 | 50 | 75 | 100 | |
|---|---|---|---|---|---|---|
| Core2Duo + Xp | jmp | 2.28 | 9.65 | 12.82 | 9.44 | 2.27 |
| setg + movzx | 3.02 | 3.02 | 3.02 | 3.02 | 3.02 | |
| carry | 2.40 | 2.41 | 2.41 | 2.40 | 2.40 | |
| setg | 2.35 | 2.35 | 2.35 | 2.35 | 2.35 | |
| Core i7 + Win7 | jmp | 2.62 | 9.33 | 12.71 | 9.15 | 1.85 |
| setg + movzx | 2.61 | 2.89 | 2.65 | 2.75 | 2.60 | |
| carry | 2.78 | 2.88 | 2.70 | 2.61 | 2.64 | |
| setg | 2.64 | 2.65 | 2.72 | 2.72 | 2.65 | |
| Opteron 2378 + Win2008 | jmp | 3.04 | 7.91 | 10.90 | 8.28 | 4.04 |
| setg + movzx | 2.88 | 2.94 | 2.95 | 3.15 | 2.90 | |
| carry | 3.03 | 3.03 | 3.02 | 3.03 | 3.03 | |
| setg | 2.25 | 2.04 | 2.04 | 2.04 | 2.04 |
carryを使ったコードはsetgよりも常によいだろうと思ったが, Opteronではそうとも言い切れなかった. Opteronではcmovと同様setgは速い.
setg + movzxとsetgを比べると, Core2DuoとOpteronで気持ち速くなった. レジスタが余ってるならcarryよりいいかも.
jmpでは"割合"が大きくなるほど個数をincrementする部分のコードはスキップされる.
割合が50の時のjmpの遅さは注意すべき. 本当に外れるときはこのぐらいのコストがかかっている.
#1/2のならし平均なので恐らく20clk以上.
これを見るとcase1の最大値のときの"だいたい単調増加"では実は結構, 分岐予測があたっていたと推測される.
ところで0以上n未満の二つの自然数a, bの和がnを超える確率はほぼ1/2なので, (a + b) mod nの計算で分岐予測はできない. よって上記の"割合=50"に相当する. というわけでateペアリングの実装ではcmovを使っていた. で, そもそもmodulo演算をできるだけしないようにロジックを変えたのであった.
const uint64_t n16 = 10000000000000000ULL;
const int N = 250250;
const int M = 250;
static int vec[N + 1];
static uint64_t num[M];
uint32_t ipowMod(uint32_t x, uint32_t a, uint32_t m)
{
if (a == 0) return 1;
if (a == 1) return x % m;
uint64_t t = x;
uint64_t out = 1;
while (a > 0) {
if (a & 1) out = (uint64_t(out) * t) % m;
t = (uint64_t(t) * t) % m;
a >>= 1;
}
return out;
}
template<class T>
T addm(T a, T b, T c) { T s = a + b; return (s >= c) ? s - c : s; }
int main()
{
for (int i = 1; i <= N; i++) {
vec[i] = ipowMod(i,i,M);
}
num[0] = 1;
for (int i = 1; i <= N; i++) {
int v = vec[i];
uint64_t ad[M];
std::copy(num, num + M, ad);
for (int j = 0; j < M; j++) {
int p = addm(v, j, M);
num[p] = addm(ad[p], ad[j], n16);
}
}
printf("%llu\n", num[0] - 1);
}
でHaskellでも書いてみた.
import Data.List
import Data.Array
import Data.Word
n=250250
m=250
n16::Word64
n16=10^16
pow mul x n = f x (n-1) x
where f _ 0 y = y
f x n y = g x n where
g x n | even n = g (mul x x) (n`quot`2)
| otherwise = f x (n-1) (mul x y)
{- x^n mod m -}
ipowMod x n m = pow mul x n
where mul x y = rem (x * y) m
vec::Array Int Int
vec = listArray (1,n).map f $ [1..n]
where m'::Word64
m' = fromIntegral m
f i = fromIntegral $ ipowMod i' i' m'
where i'::Word64
i' = fromIntegral i
{- a + b mod c -}
addm a b c = let s = a + b in if s >= c then s - c else s
up num i = array (0,m-1) $ map f [0..m-1]
where f j = (p, addm (num ! p) (num ! j) n16)
where p = addm (vec ! i) j m
num::Array Int Word64
num=listArray (0,m-1).take m $ 1:[0,0..]
calc = foldl' up num [1..n]
ans = (calc ! 0) - 1
main = print $ ans
n=25025と一桁落とした場合は2秒ぐらいだったのでこのままでもいけるかと思ったらn=250250ではメモリ不足で落ちた.
foldl'でstrict評価してるのになんでメモリが足りなくなるのだろう. GC回収してくれ.
悩むのも面倒なのでSTArrayを使ってup以下を書き直す.
import Control.Monad.ST
import Data.Array.ST
{- 上と同じ -}
up num i = do
let v = vec ! i
ad <- newArray_ (0,m-1) :: ST s (STArray s Int Word64)
mapM_ (cp ad) [0..m-1]
mapM_ (f v num ad) [0..m-1]
where cp a j = do
x<- readArray num j
writeArray a j $! x
f v ao ai j = do
let p = addm v j m
x <- readArray ai p
y <- readArray ai j
writeArray ao p $! (addm x y n16)
calc = do
num <- newArray (0,m-1) 0 :: ST s (STArray s Int Word64)
writeArray num 0 1
mapM_ (up num) [1..n]
getElems num
ans = runST calc
main = print $ (ans !! 0) - 1
これなら7.36秒でC++の10倍遅い程度ですんだ. でもどうみてもC++のループの方が簡単.
速くてわかりやすく書く方法が知りたいなあ. [2011/4/13 追記高速化遍]
import Data.List
import Data.Array
limit=2^50
{- primesは素数の無限列, isqrtは整数の平方根(略) -}
ps = map (^2).takeWhile(<isqrt limit) $ primes
pa = listArray (0,(length ps)-1) ps
f::Integer->Int->Integer
f 0 _ = 0
f n j | n == (fromIntegral j-1) = n
| otherwise = n - sum [f n' (i+1)|(n',i)<-sub]
where sub = takeWhile((/=0).fst) [(quot n (pa ! i), i)|i<-[j..(length ps)-1]]
main = print $ f limit 0
typedef std::vector<uint64_t> IntVec;
/* primesは2^25までの素数の列(略) */
IntVec ps(primes.size());
uint64_t f(uint64_t n, size_t j)
{
if (n == 0) return 0;
if (n == j-1) return n;
uint64_t sum = 0;
for (size_t i = j; i < ps.size(); i++) {
uint64_t x = n / ps[i];
if (x == 0) break;
sum += f(x, i + 1);
}
return n - sum;
}
int main()
{
for (size_t i = 0; i < primes.size(); i++) {
uint64_t p = primes.get(i);
ps[i] = uint64_t(p) * p;
}
printf("%llu\n", f(limit, 0));
}