# How to tweak the following purely numerical code

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, {...}}

-
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 '14 at 0:20
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 '14 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 '14 at 15:26
@m_goldberg good call, after doing as you suggested, the execution time went down slightly to 0.082000. – ftiaronsem May 28 '14 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 '14 at 15:29