After 8 years, I am again struggling with my matrix B to get it still faster... The problem is that the code involves big numbers generated by Binomial[t, i]. For instance, Binomial[t, t/2]=2.04815162699310^600 for t=2000. So this is well above double floating of 2^1023=8.9884710^307. In the code, there are also factorials t!, Anyone has an idea how such numbers can be handled to prevent the compiled code to generate errors and reverse back to the slow uncompiled code?

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} *)

Thank you Daniel and Mariusz for your suggestions. Meanwhile I found that the following code is the fastest version:

  Bmatrix[a_, t_, m_] := Block[{n, p, q, B},
   n = 1 - m; p = m/a; q = m (1 - 1/a); B = ConstantArray[0., {t, t}];
   Do[B[[i, j]] = 
     j! /i! p^(-i + j) (n + p)^(-j + t) (n + q)^
      i Hypergeometric2F1Regularized[-i, j - t, 1 - i + j, (
       p q)/((n + p) (n + q))], {i, 1, Floor[t/2]}, {j, i, t - i}];
   Do[B[[i, j]] = 
     j! /i! p^(-i + j) (n + p)^(-j + t) (n + q)^
      i Hypergeometric2F1Regularized[-i, j - t, 1 - i + j, (
       p q)/((n + p) (n + q))], {j, Ceiling[t/2], t - 1}, {i, 
     t - j + 1, j}];
   Do[B[[i, t]] = p^(t - i) (n + q)^i Binomial[t, i], {i, 1, t}];
   Do[B[[i, j]] = (q/p)^(i - j) Binomial[t, i]/
      Binomial[t, j] B[[j, i]], {i, 2, t}, {j, 1, i - 1}];
   Return[B]];

It can be C-compiled but j! /i! and similar expressions produce too big numbers. To prevent this, I could use

((p^j) j! )/(p^i i!) = FactorialPower[p j, j - i, p]

in the first Do loop because there j>=i and 0<p<<1. But this doesn't help because strangely the thus compiled code is slower than the uncompiled one...

Edit: I just found out that Product[p (j - k), {k, 0, j - i - 1}] is much faster than FactorialPower[p j, j - i, p] if compiled. Shouldn't built-in functions be more rapid?