I am trying to solve https://projecteuler.net/problem=499 by constructing a transition matrix and finding its steady state. Here is my code:
pmatrix[s_, m_, size_, iteration_] := Module[{matrix, matrix2},
matrixtransformer[n_] := 0;
Do[matrixtransformer[2^(i - 1) - m] = N[1/2^i,20], {i,
Floor[Log[2, size]]}];
matrixtransformer2[n_] := 0; matrixtransformer2[0] = 1;
matrix = Join[
Map[matrixtransformer2,
Table[Table[i - j, {i, size}], {j, m - 1}], {2}],
Map[matrixtransformer,
Table[Table[i - j, {i, size}], {j, m, size}], {2}]];
(*construct the initial matrix*)
fastmul[n_, t_] :=
Module[{matrix3}, matrix3 = n;
Do[matrix3 = Dot[matrix3, matrix3], {t}]; matrix3];
matrix = fastmul[matrix, iteration];
matrix2 =
Map[matrixtransformer2, Table[s - i, {i, 1, size}]].matrix;
N[1 - Plus @@ matrix2[[1 ;; m - 1]], 20]]
It works just fine with small parameters like $pmatrix[2,2,100,500]=0.25220849$ However, it is impossible to do the problem with matrix size of $10^9$. What can be done to improve the code? (PS:I am new with Wolfram Language, any suggestions on optimising the code will be appreciated also)
Update: I wrote a new algorithm seems to have a much better time complexity. But the ten thousand case would take a long time also, I don't know where the problem is. Here is my code:
pmatrix2[s_, m_, acc1_, acc2_, iteration_] :=
Module[{}, ptemp[n_] := 0;
ft[n_] := Sum[1/2^i ptemp[n + 2^(i - 1) - m], {i, acc2}];
Do[ptemp[i] = 1, {i, m - 1}];
Do[Do[ptemp[i] = ft[i], {i, m, acc1}], {iteration}];
1 - ptemp[s] // N]