>>684
512✕512の遷移行列を作って計算
慣れたRで書いたので結果は小数表示。
> cbind(n,'P[n]'=solve(n))
n P[n]
[1,] 5 0.002032211
[2,] 7 0.002985593
[3,] 9 0.003421399
[4,] 11 0.003635893
[5,] 13 0.003749860
[6,] 15 0.003813862
[7,] 17 0.003851047
[8,] 19 0.003873067
[9,] 21 0.003886242

コードのサラダ
library(expm)
n2b <- function(num, N=2, digit = 9){
r=num%%N
q=num%/%N
while(q > 0 | digit > 1){
r=append(q%%N,r)
q=q%/%N
digit=digit-1
}
return(rev(r))
}
b2n=\(n) (n %*% 2^(0:8))[1,1]
M=matrix(0,ncol=512,nrow=512)
colnames(M)=0:511
rownames(M)=0:511
a1=c(1,1,0,1,0,0,0,0,0)
a2=c(1,1,1,0,1,0,0,0,0)
a3=c(0,1,1,0,0,1,0,0,0)
a4=c(1,0,0,1,1,0,1,0,0)
a5=c(0,1,0,1,1,1,0,1,0)
a6=c(0,0,1,0,1,1,0,0,1)
a7=c(0,0,0,1,0,0,1,1,0)
a8=c(0,0,0,0,1,0,1,1,1)
a9=c(0,0,0,0,0,1,0,1,1)
A=rbind(a1,a2,a3,a4,a5,a6,a7,a8,a9)
f1=\(i) apply(A,1,\(v) b2n(bitwXor(n2b(i),v)))
for(i in 0:511) M[as.character(i),as.character(f1(i))]=1/9
solve=\(n) (M %^%n)["0","511"]
solve=Vectorize(solve)
n=c(5,7,9,11,13,15,17,19,21)
cbind(n,'P[n]'=solve(n))