Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

I am using the following numerical function in my program and I would like to speed it up.

calculate[mm_, nn_, pp_, ff_, quiet_] := 
  Module[{n, y, m, z, f, k, p, l, \[Tau]},
   h[0, 0, 0] = 1;
   h[n_, 0, y_] := If[Mod[n, 2] == 0,
     n!*y^(n/2)/((n/2)!),
     0];

   If[quiet != 1, 
    If[Simplify[h[4, 0, -1] - HermiteH[4, 0]] === 0, , 
     Print["Bogus definition of h[n_,0,y_]"]]];
   If[quiet != 1, 
    If[Simplify[h[3, 0, -1] - HermiteH[3, 0]] === 0, , 
     Print["Bogus definition of h[n_,0,y_]"]]];

   hh[m_, n_, 0, y_, 0, z_, \[Tau]_] := 
    Sum[m!*n!/((m - k)!*(n - k)!*k!)*(\[Tau])^k*h[m - k, 0, y]*
      h[n - k, 0, z], {k, 0, Min[m, n]}];

   r[m_, n_, k_, f_] := (1/f)^k*
     Sum[Binomial[k, l]* m!*n!/((m - l)!*(n - k + l)!)*
       hh[m - l, n - k + l, 0, -1 + 1/f, 0, -1 + 1/f, 2/f], {l, 0, k}];

   integral[m_, n_, p_, f_] := 
    Sqrt[N[Pi]/f]*
     Sum[Binomial[p, k]*h[p - k, 0, 1/(4*f)]*r[m, n, k, f], {k, 0, p}];

   Return[integral[mm, nn, pp, ff]];
   ];

Performance isn't too bad at the moment,

Timing[Table[calculate[n, n, 0, f, 0], {n, 0, 100}]]
{0.096000, {....}}

but the function is called millions of times.

All the input to the function will be solemnly numeric. I don't need any symbolic evaluation and machine precision would be just fine. So I figured that there would be plenty of room for improvement, but I am fairly new to Mathematica and especially to its numerical capabilities.

I would greatly appreciate any suggestions on speeding up that code. Thanks so much.

Note: after taking into account the comments, this is my current code:

h[0., 0., 0.] = 1;
h[n_, 0., y_] := If[Mod[n, 2] == 0, n!*y^(n/2)/((n/2)!), 0];

hh[m_, n_, 0, y_, 0, z_, \[Tau]_] := 
  Sum[m!*n!/((m - k)!*(n - k)!*k!)*(\[Tau])^k*h[m - k, 0., y]*
    h[n - k, 0., z], {k, 0., Min[m, n]}];

r[m_, n_, k_, f_] := (1/f)^k*
   Sum[Binomial[k, l]*m!*n!/((m - l)!*(n - k + l)!)*
     hh[m - l, n - k + l, 0, -1 + 1/f, 0, -1 + 1/f, 2/f], {l, 0., 
     k}];

calculate[m_, n_, p_, f_] := 
  Sqrt[N[Pi]/f]*
   Sum[Binomial[p, k]*h[p - k, 0., 1/(4*f)]*r[m, n, k, f], {k, 0., 
     p}];

and the performance is

f = 1.;
Timing[Table[integrateGaussAndTwoHermite[n, n, 0, f], {n, 0, 100}]]
{0.080000, {...}}

Note: After taking into account even more comments, my code is

h[0., 0., 0.] = 1;
h[n_, 0., y_] := 
  If[Mod[n, 2] == 0, Gamma[n + 1.]*y^(n/2.)/((n/2.)!), 0];

hh[m_, n_, 0., y_, 0., z_, \[Tau]_] := 
  Sum[Binomial[m, k]*Gamma[n + 1.]/(Gamma[n - k + 1.])*(\[Tau])^k*
    h[m - k, 0., y]*h[n - k, 0., z], {k, 0., Min[m, n]}];

r[m_, n_, k_, f_] := (1./f)^k*
   Sum[Binomial[k, l]*Gamma[m + 1.]*
     Gamma[n + 1.]/(Gamma[m - l + 1.]*Gamma[n - k + l + 1.])*
     hh[m - l, n - k + l, 0., -1. + 1./f, 0., -1. + 1./f, 2./f], {l, 
     0., k}];

integrateGaussAndTwoHermite[m_, n_, p_, f_] := 
  Sqrt[N[Pi]/f]*
   Sum[Binomial[p, k]*h[p - k, 0., 1./(4.*f)]*r[m, n, k, f], {k, 0., 
     p}];

and the execution time is

{0.0720000, {...}}
share|improve this question
2  
You say you don't need symbolic evaluation, but calculate modifies global symbols every time it's called. Also, f doesn't seem to be defined anywhere. I think we could help you better if you describe what you want your code to do, rather than what the code currently does. –  phosgene May 28 at 0:20
2  
I don't think you need to define calculate at all. Simply move the functions you define inside calculate to top-level, and call integral wherever you need to in the rest of your code. You really don't want to define all those functions at every call to calculate. –  m_goldberg May 28 at 0:20
    
Thanks a lot for your comments. @phosgene The routine is the implementation of the analytic result of a special case of an integral involving some Gaussians and HermitePolynomials. And f=1 in my testcases. –  ftiaronsem May 28 at 15:26
    
@m_goldberg good call, after doing as you suggested, the execution time went down slightly to 0.082000. –  ftiaronsem May 28 at 15:26
    
get rid of the debugging code. You are running quiet!=1 but it will speed up a tad to not check that every time. –  george2079 May 28 at 15:29

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Browse other questions tagged or ask your own question.