I rewrote a bit. This still uses seq
but not as a Sequence
object. It has a final clause that is never actually needed.
B2[a_, m_, t_] := Module[{n, p, q, seq, s, f, b},
n = 1 - m;
p = m/a;
q = m - p;
s[i_, j_, K_] :=
Sum[Binomial[i, i - k] Binomial[t - i, j - i + k] (p q)^
k ((n + p) (n + q))^(K - k), {k, 0, K}];
f[i_, j_] :=
If[(t - i - j) >= 0, (n + p)^(t - i - j), (n + q)^-(t - i - j)];
seq[i_, j_] := Which[
i <= Round[t/2] && t - j <= Round[t/2],
s[i, j, Min[i, t - j]]
,
i <= Round[t/2],
s[i, j, i]
,
t - j <= Round[t/2],
s[i, j, t - j]
,
True, 0];
b[i_, j_] := Which[
j == t,
p^(j - i) f[i, j] Binomial[t, i]
,
i == t,
q^(i - j) f[i, j] Binomial[t, j]
,
i <= j,
p^(j - i) f[i, j] Binomial[t, i] seq[i, j]
,
i > j,
q^(i - j) f[i, j] Binomial[t, j] seq[j, i]];
Table[b[i, j], {i, 1, t}, {j, 1, t}]]
As you note, this can get into trouble with machine numbers.
Timing[B2[3, .22, 400];]
(* During evaluation of In[127]:= General::munfl: 3.75211*10^-299 5.02373*10^-10 is too small to represent as a normalized machine number; precision may be lost.
During evaluation of In[127]:= General::munfl: 2.75154*10^-300 5.88718*10^-10 is too small to represent as a normalized machine number; precision may be lost.
During evaluation of In[127]:= General::munfl: 2.0178*10^-301 6.89904*10^-10 is too small to represent as a normalized machine number; precision may be lost.
During evaluation of In[127]:= General::stop: Further output of General::munfl will be suppressed during this calculation.
Out[127]= {30.3564, Null} *)
It handles higher precision without signs of distress though, and this is not hugely slower.
Timing[b2big = B2[3, .22`20, 400];]
(* Out[129]= {51.2323, Null} *)