Numerical-Symbolical Integration (Calculus)

Question

I created a simple numeric-symbolic integration. Here you can use symbolical and numerical techniques at the same time. You can also interpolate numerical integrals.

The problem with my function is the use of TimeConstrained for testing if the integral can be done fully symbolically. This wastes time; without it the result is much much faster.

How can it be optimized? Maybe there is some other implementation like this I haven't found? Note: It can easily be extended to Sum-s, just exchange Integrate for Sum :)

Options@NSIntegrate = {
    InterpolationVariable -> {}
};

(* Linearity *)
NSIntegrate[f_ + g__, s_, options___] := NSIntegrate[f, s, options] + NSIntegrate[+g, s, options]
NSIntegrate[α_*f_, {x_, limits__}, options___] /; FreeQ[α, x] := α*NSIntegrate[f, {x, limits}, options]

(* Constant Function *)
NSIntegrate[α_, {x_, limits__}, options___] /; FreeQ[α, x] \[And] α =!= 1 := α*NSIntegrate[1, {x, limits}, options]

(* Exponentials *)
NSIntegrate[α_^(f_+g__), {x_, limits__}, options___] /; FreeQ[f, x] := α^f*NSIntegrate[α^g, {x, limits}, options]
NSIntegrate[α_^(f_+g__)*h_, {x_, limits__}, options___] /; FreeQ[f, x] := α^f*NSIntegrate[α^g*h, {x, limits}, options]

(* Some Basic Integrals - Polynomials *)
NSIntegrate[1, {x_, a_, b_}, ___] := b - a
NSIntegrate[x_, {x_, a_, b_}, ___] := (b^2 - a^2)/2
NSIntegrate[x_^-1, {x_, a_, b_}, ___] /; If[NumericQ@a \[And] NumericQ@b, Sign@a === Sign@b, True] \[And] a =!= b := Log[b/a]
NSIntegrate[x_^n_, {x_, a_, b_}, ___] /; n =!= -1 := (b^(n+1) - a^(n+1))/(n + 1)

(* Some Basic Integrals - Exponentials *)
NSIntegrate[α_^(A_*x_), {x_, a_, b_}, ___] /; α =!= 0 \[And] α =!= 1 := (α^(A*b) - α^(A*a))/(A*Log@α)
NSIntegrate[α_^x_, {x_, a_, b_}, ___] /; α =!= 0 \[And] α =!= 1 := (α^b - α^a)/Log@α
NSIntegrate[0^x_, s__, ___] := 0

(* Scaling and offset *)
NSIntegrate[f_[α_*x_ + β_, args___], {x_, a_, b_}, options___] /; (And @@ (FreeQ[#, x] & /@ Flatten@{args, α, β})) := 
    NSIntegrate[f[x, args]/α, {x, a*α + β, b*α + β}, options]
NSIntegrate[f_[α_*x_, args___], {x_, a_, b_}, options___] /; (And @@ (FreeQ[#, x] & /@ Flatten@{args, α})) := 
    NSIntegrate[f[x, args]/α, {x, a*α, b*α}, options]
NSIntegrate[f_[x_ + β_, args___], {x_, a_, b_}, options___] /; (And @@ (FreeQ[#, x] & /@ Flatten@{args, β})) := 
NSIntegrate[f[x, args], {x, a + β, b + β}, options]

(* Interpolation *)
NSIntegrate[f_, {x_, a_, b_}, options1___, InterpolationVariable -> varRange_List, options2___] /; Length@varRange != 0 := Module[{
    vars, range = varRange, table
    },
    If[Depth@varRange == 2, vars = varRange[[{1}]]; range = List@range, vars = varRange[[All, 1]]];

    table = Table[{vars, NSIntegrate[f, {x, a, b}, options1, options2]}, ##1] & @@ range;
    InterpolatingPolynomial[table, vars] /. s_Real :> Chop[s, 10^-6] // Expand
]

(* If not in database, use NIntegrate *)
NSIntegrate[f_, s:{x_, a_, b_}, opts___] /; ExpandAll@f =!= f := NSIntegrate[ExpandAll@f, s, opts]
NSIntegrate[f_, s:{x_, a_?NumericQ, b_?NumericQ}, opts___] := NIntegrate[f, s, opts]

*Usage:

NSIntegrate[Cos[a*x], {x, 0, 1}, InterpolationVariable -> {a, 1, 5, 1/2}]

InterpolationVariable work just like in a Table. In case or more than 1 interpolation use in the form {{a, 1, 5}, {b, 1, 10, 0.3}}

** Added (* Scaling and offset *) section.