🚀 go-pugleaf

It may also be that if you are evaluating over a grid of points,
you can save a lot of time. For example, if you are investigating
the function at a given x over a range of values of a and b, then 
note that there is a component in a only (1/Gamma(a) * x**(a-1)) and 
a component in b only (exp(-x/b)), and only a very simple component
in both (b**a). It would be silly to recompute all three things
for the whole grid - you could compute a vector of the term in a
(T(a), say), a vector for the term in b (V(b), say) and then
you only need to compute b**a each time: f = T(a)*V(b)*b**a.
Even better, you can save a lot of time computing that b**a term 
by the fact that 'a' is over a grid, so b**a(n+1) = K * b**a(n); 
errors will slowly build up, but you could do it on widely spaced
values of a (say, every 16th or 8th value). You might even check
the centre estimate in the middle of each range both from above
and below as a check on loss of accuracy, though more formal
estimation and control of accuracy would also be possible.

[If you were instead doing it over a grid of x, then the constant 
term is ... constant, and you get that the term in exp(-x/b) is a
constant multiple of the same term for the previous x (in which
case similar comments apply as for the b**a term).  If you have
a grid for x, a and b then you're really flying.]

Glen