Group Abstract Group Abstract

Message Boards Message Boards

Make code faster and in C-compatible manner?

Posted 9 years ago
POSTED BY: Ulrich Utiger
7 Replies

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} *)
POSTED BY: Daniel Lichtblau

If you change Sum to: Total@Table** code is 20%** faster.

POSTED BY: Mariusz Iwaniuk
Posted 1 year ago

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?

POSTED BY: Ulrich Utiger

You might experiment with taking logs, PowerExpanding using LogGamma instead of Log[Factorial[...]], and exponentiating at the end. If there is intermediate coefficient growth followed by cancelation, the LogGamma usage should ameliorate the issue numerically.

POSTED BY: Daniel Lichtblau
Posted 1 year ago

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?

POSTED BY: Ulrich Utiger

Maybe post your question here. They will definitely answer your question and you will be satisfied

I'm not a good programmer who knows tricks to speed up the code and I won't help you in your situation.

Regards M.I.

POSTED BY: Mariusz Iwaniuk
Posted 9 years ago
POSTED BY: Ulrich Utiger
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard