using Symata;

 a = 3.0
NIntegrate(Exp(-t)*t^(a-1), [t,0,Infinity])

gamma(a_, x_:Infinity) := NIntegrate(Exp(-t)*t^(a-1), [t,0,x])[1]

Map(gamma, [3.0, 4.0, 5.0])

 gammaint = Compile([t] , Exp(-t)*t^(a-1));

gamma1(a1_, x_:Infinity) := (SetJ(a,a1), NIntegrate(gammaint, [0,x]))[1]

gamma1(3.0)

Timing(gamma(3.0))

Timing(gamma1(3.0))

VersionInfo()

InputForm(Now())