f(x):=x回目でONになっている確率とすると、漸化式は
f(x)=k∑x−1k回目でOFFで、あとはずっとONである確率=k∑x−1(1−f(k))p(1−p)x−k−1
このままでは間に合わないので差分を計算すると、f(x)−f(x−1)=p(1−2f(x−1))になっていて、f(x)=p+(1−2p)f(x−1)という線形漸化式が従う。f(x)−1/2=(1−2p)(f(x−1)−1/2)より、f(x)=1/2+(1−2p)x(f(0)−1/2)=1/2+(1−2p)x⋅(−1/2)と表せて、これなら計算できる形になっている。
……と迂遠な導出をしてしまったが、よく見ると自明にf(x)=p(1−f(x−1))+(1−p)f(x−1)だった。直前から求まるかどうかは常に最初に考えるべき。
p=1の場合に計算誤差が問題になるようだった。(expt -1d0 奇数)
は-1になってほしいが、指数が大きいと1になる。手元だと境界はこの辺:
> (expt -1d0 9007199254740991)
-1.0d0
> (expt -1d0 9007199254740993)
1.0d0
前者(=253−1)がdouble-float
で正確に表現できる最大の整数なので、そこを超えると精度が落ちるのは明らかとして、そこまではOKなのかは実装次第だが、たぶん普通はOK。