diff --git a/src/Interpolations.jl b/src/Interpolations.jl
index f1c50cad..c6f0d61e 100644
--- a/src/Interpolations.jl
+++ b/src/Interpolations.jl
@@ -31,21 +31,24 @@ export
 
     degree,
     boundarycondition,
-    gridrepresentation
+    gridrepresentation,
+
+    gradient,
+    gradient!
 
 abstract Degree{N}
 
 abstract GridRepresentation
-type OnGrid <: GridRepresentation end
-type OnCell <: GridRepresentation end
+immutable OnGrid <: GridRepresentation end
+immutable OnCell <: GridRepresentation end
 
 abstract BoundaryCondition
-type None <: BoundaryCondition end
-type Flat <: BoundaryCondition end
-type Line <: BoundaryCondition end
+immutable None <: BoundaryCondition end
+immutable Flat <: BoundaryCondition end
+immutable Line <: BoundaryCondition end
 typealias Natural Line
-type Free <: BoundaryCondition end
-type Periodic <: BoundaryCondition end
+immutable Free <: BoundaryCondition end
+immutable Periodic <: BoundaryCondition end
 immutable Reflect <: BoundaryCondition end
 
 abstract InterpolationType{D<:Degree,BC<:BoundaryCondition,GR<:GridRepresentation}
@@ -53,60 +56,69 @@ abstract InterpolationType{D<:Degree,BC<:BoundaryCondition,GR<:GridRepresentatio
 include("extrapolation.jl")
 
 abstract AbstractInterpolation{T,N,IT<:InterpolationType,EB<:ExtrapolationBehavior} <: AbstractArray{T,N}
-type Interpolation{T,N,IT<:InterpolationType,EB<:ExtrapolationBehavior} <: AbstractInterpolation{T,N,IT,EB}
-    coefs::Array{T,N}
+type Interpolation{T,N,TCoefs,IT<:InterpolationType,EB<:ExtrapolationBehavior} <: AbstractInterpolation{T,N,IT,EB}
+    coefs::Array{TCoefs,N}
 end
-function Interpolation{T,N,IT<:InterpolationType,EB<:ExtrapolationBehavior}(A::Array{T,N}, it::IT, ::EB)
+function Interpolation{N,TCoefs,TWeights<:Real,IT<:InterpolationType,EB<:ExtrapolationBehavior}(::Type{TWeights}, A::AbstractArray{TCoefs,N}, it::IT, ::EB)
     isleaftype(IT) || error("The interpolation type must be a leaf type (was $IT)")
-    
-    isleaftype(T) || warn("For performance reasons, consider using an array of a concrete type T (eltype(A) == $(eltype(A)))")
 
-    Interpolation{T,N,IT,EB}(prefilter(A,it))
+    isleaftype(TCoefs) || warn("For performance reasons, consider using an array of a concrete type T (eltype(A) == $(eltype(A)))")
+
+    c = one(TWeights)
+    for _ in 2:N
+        c *= c
+    end
+    T = typeof(c * one(TCoefs))
+
+    Interpolation{T,N,TCoefs,IT,EB}(prefilter(TWeights,A,it))
 end
+Interpolation(A::AbstractArray, it::InterpolationType, eb::ExtrapolationBehavior) = Interpolation(Float64, A, it, eb)
+Interpolation(A::AbstractArray{Float32}, it::InterpolationType, eb::ExtrapolationBehavior) = Interpolation(Float32, A, it, eb)
+Interpolation(A::AbstractArray{Rational{Int}}, it::InterpolationType, eb::ExtrapolationBehavior) = Interpolation(Rational{Int}, A, it, eb)
 
 # Unless otherwise specified, use coefficients as they are, i.e. without prefiltering
 # However, all prefilters copy the array, so do that here as well
-prefilter{T,N,IT<:InterpolationType}(A::AbstractArray{T,N}, ::IT) = copy(A)
+prefilter{TWeights,T,N,IT<:InterpolationType}(::Type{TWeights}, A::AbstractArray{T,N}, ::IT) = copy(A)
 
-size{T,N,IT<:InterpolationType}(itp::Interpolation{T,N,IT}, d::Integer) =
-    size(itp.coefs, d) - 2*padding(IT())
-size(itp::Interpolation) = tuple([size(itp,i) for i in 1:ndims(itp)]...)
+size{T,N}(itp::Interpolation{T,N}, d::Integer) = d > N ? 1 : size(itp.coefs, d) - 2*padding(interpolationtype(itp))
+size(itp::AbstractInterpolation) = tuple([size(itp,i) for i in 1:ndims(itp)]...)
 ndims(itp::Interpolation) = ndims(itp.coefs)
-eltype(itp::Interpolation) = eltype(itp.coefs)
 
 offsetsym(off, d) = off == -1 ? symbol(string("ixm_", d)) :
                     off ==  0 ? symbol(string("ix_", d)) :
                     off ==  1 ? symbol(string("ixp_", d)) :
                     off ==  2 ? symbol(string("ixpp_", d)) : error("offset $off not recognized")
 
+interpolationtype{T,N,TCoefs,IT<:InterpolationType}(itp::Interpolation{T,N,TCoefs,IT}) = IT()
+
 boundarycondition{D,BC<:BoundaryCondition}(::InterpolationType{D,BC}) = BC()
-boundarycondition{T,N,IT}(::Interpolation{T,N,IT}) = boundarycondition(IT())
 gridrepresentation{D,BC,GR<:GridRepresentation}(::InterpolationType{D,BC,GR}) = GR()
-gridrepresentation{T,N,IT}(::Interpolation{T,N,IT}) = gridrepresentation(IT())
 degree{D<:Degree,BC,GR}(::InterpolationType{D,BC,GR}) = D()
-degree{T,N,IT}(::Interpolation{T,N,IT}) = degree(IT())
+
+# If nothing else is specified, don't pad at all
+padding(::InterpolationType) = 0
+
+for op in (:boundarycondition, :gridrepresentation, :degree, :padding)
+    eval(:($(op)(itp::Interpolation) = $(op)(interpolationtype(itp))))
+end
 
 include("constant.jl")
 include("linear.jl")
 include("quadratic.jl")
 
-# If nothing else is specified, don't pad at all
-padding(::InterpolationType) = 0
-padding{T,N,IT<:InterpolationType}(::Interpolation{T,N,IT}) = padding(IT())
-
-function pad_size_and_index(sz::Tuple, pad)
+function padded_index(sz::Tuple, pad)
     sz = Int[s+2pad for s in sz]
     N = length(sz)
     ind = cell(N)
     for i in 1:N
         ind[i] = 1+pad:sz[i]-pad
     end
-    sz, ind
+    ind,sz
 end
 function copy_with_padding(A, it::InterpolationType)
     pad = padding(it)
-    sz,ind = pad_size_and_index(size(A), pad)
-    coefs = fill(convert(eltype(A), 0), sz...)
+    ind,sz = padded_index(size(A), pad)
+    coefs = zeros(eltype(A), sz...)
     coefs[ind...] = A
     coefs, pad
 end
@@ -139,11 +151,9 @@ for IT in (
         it = IT()
         eb = EB()
         gr = gridrepresentation(it)
-        eval(ngenerate(
-            :N,
-            :(promote_type(T, x...)),
-            :(getindex{T,N}(itp::Interpolation{T,N,$IT,$EB}, x::NTuple{N,Real}...)), 
-            N->quote
+
+        function getindex_impl(N)
+                quote
                 # Handle extrapolation, by either throwing an error,
                 # returning a value, or shifting x to somewhere inside
                 # the domain
@@ -160,13 +170,38 @@ for IT in (
                 @inbounds ret = $(index_gen(degree(it), N))
                 ret
             end
+        end
+
+        # Resolve ambiguity with linear indexing,
+        #   getindex{T,N}(A::AbstractArray{T,N}, i::Real)
+        eval(ngenerate(:N, :T, :(getindex{T,N,TCoefs}(itp::Interpolation{T,N,TCoefs,$IT,$EB}, x::Real)),
+            N->:(error("Linear indexing is not supported for interpolation objects"))
         ))
 
-        eval(ngenerate(
-            :N,
-            :(Array{promote_type(T,typeof(x)...),1}),
-            :(gradient!{T,N}(g::Array{T,1}, itp::Interpolation{T,N,$IT,$EB}, x::NTuple{N,Real}...)),
+        # Resolve ambiguity with real multilinear indexing
+        #   getindex{T,N}(A::AbstractArray{T,N}, NTuple{N,Real}...)
+        eval(ngenerate(:N, :T, :(getindex{T,N,TCoefs}(itp::Interpolation{T,N,TCoefs,$IT,$EB}, x::NTuple{N,Real}...)), getindex_impl))
+
+        # Resolve ambiguity with general array indexing,
+        #   getindex{T,N}(A::AbstractArray{T,N}, I::AbstractArray{T,N})
+        eval(ngenerate(:N, :T, :(getindex{T,N,TCoefs}(itp::Interpolation{T,N,TCoefs,$IT,$EB}, I::AbstractArray{T})),
+            N->:(error("Array indexing is not defined for interpolation objects."))
+        ))
+
+        # Resolve ambiguity with colon indexing for 1D interpolations
+        #   getindex{T}(A::AbstractArray{T,1}, C::Colon)
+        eval(ngenerate(:N, :T, :(getindex{T,TCoefs}(itp::Interpolation{T,1,TCoefs,$IT,$EB}, c::Colon)),
+            N->:(error("Colon indexing is not supported for interpolation objects"))
+        ))
+
+        # Allow multilinear indexing with any types
+        eval(ngenerate(:N, :(promote_type(T,TIndex)), :(getindex{T,N,TCoefs,TIndex}(itp::Interpolation{T,N,TCoefs,$IT,$EB}, x::NTuple{N,TIndex}...)), getindex_impl))
+
+        # Define in-place gradient calculation
+        #   gradient!(g::Vector, itp::Interpolation, NTuple{N}...)
+        eval(ngenerate(:N, :(Vector{TOut}), :(gradient!{T,N,TCoefs,TOut}(g::Vector{TOut}, itp::Interpolation{T,N,TCoefs,$IT,$EB}, x::NTuple{N,Any}...)),
             N->quote
+                length(g) == $N || error("g must be an array with exactly N elements (length(g) == "*string(length(g))*", N == "*string(N)*")")
                 $(extrap_transform_x(gr,eb,N))
                 $(define_indices(it,N))
                 @nexprs $N dim->begin
@@ -184,7 +219,7 @@ for IT in (
     end
 end
 
-gradient{T}(itp::Interpolation{T}, x...) = gradient!(Array(T,ndims(itp)), itp, x...)
+gradient{T,N}(itp::Interpolation{T,N}, x...) = gradient!(Array(T,N), itp, x...)
 
 # This creates prefilter specializations for all interpolation types that need them
 for IT in (
@@ -199,8 +234,8 @@ for IT in (
         Quadratic{Periodic,OnGrid},
         Quadratic{Periodic,OnCell},
     )
-    @ngenerate N promote_type_grid(T, x...) function prefilter{T,N}(A::Array{T,N},it::IT)
-        ret, pad = copy_with_padding(A,it)
+    @ngenerate N Array{TWeights,N} function prefilter{TWeights,TCoefs,N}(::Type{TWeights}, A::Array{TCoefs,N}, it::IT)
+        ret, pad = copy_with_padding(A, it)
 
         szs = collect(size(ret))
         strds = collect(strides(ret))
@@ -209,7 +244,7 @@ for IT in (
             n = szs[dim]
             szs[dim] = 1
 
-            M, b = prefiltering_system(eltype(A), n, it)
+            M, b = prefiltering_system(TWeights, TCoefs, n, it)
 
             @nloops N i d->1:szs[d] begin
                 cc = @ntuple N i
diff --git a/src/constant.jl b/src/constant.jl
index 6e029bea..15f8b22a 100644
--- a/src/constant.jl
+++ b/src/constant.jl
@@ -1,9 +1,9 @@
-type ConstantDegree <: Degree{0} end
-type Constant{GR<:GridRepresentation} <: InterpolationType{ConstantDegree,None,GR} end
+immutable ConstantDegree <: Degree{0} end
+immutable Constant{GR<:GridRepresentation} <: InterpolationType{ConstantDegree,None,GR} end
 Constant{GR<:GridRepresentation}(::GR) = Constant{GR}()
 
 function define_indices(::Constant, N)
-    :(@nexprs $N d->(ix_d = clamp(round(Int,x_d), 1, size(itp,d))))
+    :(@nexprs $N d->(ix_d = clamp(round(real(x_d)), 1, size(itp,d))))
 end
 
 function coefficients(c::Constant, N)
@@ -12,12 +12,12 @@ end
 
 function coefficients(::Constant, N, d)
     sym, symx = symbol(string("c_",d)), symbol(string("x_",d))
-    :($sym = one(typeof($symx)))
+    :($sym = 1)
 end
 
 function gradient_coefficients(::Constant, N, d)
     sym, symx = symbol(string("c_",d)), symbol(string("x_",d))
-    :($sym = zero(typeof($symx)))
+    :($sym = 0)
 end
 
 function index_gen(degree::ConstantDegree, N::Integer, offsets...)
diff --git a/src/extrapolation.jl b/src/extrapolation.jl
index 9e5395a1..02461739 100644
--- a/src/extrapolation.jl
+++ b/src/extrapolation.jl
@@ -1,30 +1,30 @@
 abstract ExtrapolationBehavior
 
-type ExtrapError <: ExtrapolationBehavior end
+immutable ExtrapError <: ExtrapolationBehavior end
 function extrap_transform_x(::OnGrid, ::ExtrapError, N)
     quote
-        @nexprs $N d->(1 <= x_d <= size(itp,d) || throw(BoundsError()))
+        @nexprs $N d->(1 <= real(x_d) <= size(itp,d) || throw(BoundsError()))
     end
 end
 function extrap_transform_x(::OnCell, ::ExtrapError, N)
     quote
-        @nexprs $N d->(.5 <= x_d <= size(itp,d) + .5 || throw(BoundsError()))
+        @nexprs $N d->(1//2 <= real(x_d) <= size(itp,d) + 1//2 || throw(BoundsError()))
     end
 end
 
-type ExtrapNaN <: ExtrapolationBehavior end
+immutable ExtrapNaN <: ExtrapolationBehavior end
 function extrap_transform_x(::OnGrid, ::ExtrapNaN, N)
     quote
-        @nexprs $N d->(1 <= x_d <= size(itp,d) || return convert(T, NaN))
+        @nexprs $N d->(1 <= real(x_d) <= size(itp,d) || return convert(T, NaN))
     end
 end
 function extrap_transform_x(::OnCell, ::ExtrapNaN, N)
     quote
-        @nexprs $N d->(.5 <= x_d <= size(itp,d) + .5 || return convert(T, NaN))
+        @nexprs $N d->(1//2 <= real(x_d) <= size(itp,d) + 1//2 || return convert(T, NaN))
     end
 end
 
-type ExtrapConstant <: ExtrapolationBehavior end
+immutable ExtrapConstant <: ExtrapolationBehavior end
 function extrap_transform_x(::OnGrid, ::ExtrapConstant, N)
     quote
         @nexprs $N d->(x_d = clamp(x_d, 1, size(itp,d)))
@@ -32,11 +32,11 @@ function extrap_transform_x(::OnGrid, ::ExtrapConstant, N)
 end
 function extrap_transform_x(::OnCell, ::ExtrapConstant, N)
     quote
-        @nexprs $N d->(x_d = clamp(x_d, .5, size(itp,d)+.5))
+        @nexprs $N d->(x_d = clamp(x_d, 1//2, size(itp,d)+1//2))
     end
 end
 
-type ExtrapReflect <: ExtrapolationBehavior end
+immutable ExtrapReflect <: ExtrapolationBehavior end
 function extrap_transform_x(::OnGrid, ::ExtrapReflect, N)
     quote
         @nexprs $N d->begin
@@ -80,23 +80,12 @@ function extrap_transform_x(::OnCell, ::ExtrapReflect, N)
     end
 end
 
-type ExtrapPeriodic <: ExtrapolationBehavior end
+immutable ExtrapPeriodic <: ExtrapolationBehavior end
 function extrap_transform_x(::GridRepresentation, ::ExtrapPeriodic, N)
-    quote
-        @nexprs $N d->begin
-            # translate x_d to inside the domain
-            n = convert(typeof(x_d), size(itp,d))
-            while x_d < .5
-                x_d += n
-            end
-            while x_d >= n + .5
-                x_d -= n
-            end
-        end
-    end
+    :(@nexprs $N d->(x_d = mod1(x_d, size(itp,d))))
 end
 
-type ExtrapLinear <: ExtrapolationBehavior end
+immutable ExtrapLinear <: ExtrapolationBehavior end
 function extrap_transform_x(::OnGrid, ::ExtrapLinear, N)
     quote
         @nexprs $N d->begin
diff --git a/src/linear.jl b/src/linear.jl
index 0543db2e..661dc554 100644
--- a/src/linear.jl
+++ b/src/linear.jl
@@ -1,13 +1,13 @@
-type LinearDegree <: Degree{1} end
-type Linear{GR<:GridRepresentation} <: InterpolationType{LinearDegree,None,GR} end
+immutable LinearDegree <: Degree{1} end
+immutable Linear{GR<:GridRepresentation} <: InterpolationType{LinearDegree,None,GR} end
 Linear{GR<:GridRepresentation}(::GR) = Linear{GR}()
 
 function define_indices(::Linear, N)
     quote
         @nexprs $N d->begin
-            ix_d = clamp(floor(Int,x_d), 1, size(itp,d)-1)
+            ix_d = clamp(floor(real(x_d)), 1, size(itp,d)-1)
             ixp_d = ix_d + 1
-            fx_d = x_d - convert(typeof(x_d), ix_d)
+            fx_d = x_d - ix_d
         end
     end
 end
@@ -19,7 +19,7 @@ end
 function coefficients(::Linear, N, d)
     sym, symp, symfx = symbol(string("c_",d)), symbol(string("cp_",d)), symbol(string("fx_",d))
     quote
-        $sym = one(typeof($symfx)) - $symfx
+        $sym = 1 - $symfx
         $symp = $symfx
     end
 end
@@ -27,8 +27,8 @@ end
 function gradient_coefficients(::Linear,N,d)
     sym, symp, symfx = symbol(string("c_",d)), symbol(string("cp_",d)), symbol(string("fx_",d))
     quote
-        $sym = -one(typeof($symfx))
-        $symp = one(typeof($symfx))
+        $sym = -1
+        $symp = 1
     end
 end
 
diff --git a/src/quadratic.jl b/src/quadratic.jl
index 1fe07237..03e72a89 100644
--- a/src/quadratic.jl
+++ b/src/quadratic.jl
@@ -1,5 +1,5 @@
-type QuadraticDegree <: Degree{2} end
-type Quadratic{BC<:BoundaryCondition,GR<:GridRepresentation} <: InterpolationType{QuadraticDegree,BC,GR} end
+immutable QuadraticDegree <: Degree{2} end
+immutable Quadratic{BC<:BoundaryCondition,GR<:GridRepresentation} <: InterpolationType{QuadraticDegree,BC,GR} end
 
 Quadratic{BC<:BoundaryCondition,GR<:GridRepresentation}(::BC, ::GR) = Quadratic{BC,GR}()
 
@@ -7,11 +7,11 @@ function define_indices(q::Quadratic, N)
     quote
         pad = padding($q)
         @nexprs $N d->begin
-            ix_d = clamp(round(Integer, x_d), 1, size(itp,d)) + pad
+            ix_d = clamp(round(real(x_d)), 1, size(itp,d)) + pad
             ixp_d = ix_d + 1
             ixm_d = ix_d - 1
 
-            fx_d = x_d - convert(typeof(x_d), ix_d - pad)
+            fx_d = x_d - (ix_d - pad)
         end
     end
 end
@@ -19,11 +19,11 @@ function define_indices(q::Quadratic{Periodic}, N)
     quote
         pad = padding($q)
         @nexprs $N d->begin
-            ix_d = clamp(round(Integer, x_d), 1, size(itp,d)) + pad
+            ix_d = clamp(round(real(x_d)), 1, size(itp,d)) + pad
             ixp_d = mod1(ix_d + 1, size(itp,d))
             ixm_d = mod1(ix_d - 1, size(itp,d))
 
-            fx_d = x_d - convert(typeof(x_d), ix_d - pad)
+            fx_d = x_d - (ix_d - pad)
         end
     end
 end
@@ -36,9 +36,9 @@ function coefficients(q::Quadratic, N, d)
     symm, sym, symp =  symbol(string("cm_",d)), symbol(string("c_",d)), symbol(string("cp_",d))
     symfx = symbol(string("fx_",d))
     quote
-        $symm = .5 * ($symfx - .5)^2
-        $sym  = .75 - $symfx^2
-        $symp = .5 * ($symfx + .5)^2
+        $symm = 1//2 * ($symfx - 1//2)^2
+        $sym  = 3//4 - $symfx^2
+        $symp = 1//2 * ($symfx + 1//2)^2
     end
 end
 
@@ -46,9 +46,9 @@ function gradient_coefficients(q::Quadratic, N, d)
     symm, sym, symp =  symbol(string("cm_",d)), symbol(string("c_",d)), symbol(string("cp_",d))
     symfx = symbol(string("fx_",d))
     quote
-        $symm = $symfx-.5
-        $sym = -2*$symfx
-        $symp = $symfx+.5
+        $symm = $symfx - 1//2
+        $sym = -2 * $symfx
+        $symp = $symfx + 1//2
     end
 end
 
@@ -75,81 +75,87 @@ padding(::Quadratic) = 1
 padding(::Quadratic{Periodic}) = 0
 
 function inner_system_diags{T}(::Type{T}, n::Int, ::Quadratic)
-    du = fill(convert(T,1/8), n-1)
-    d = fill(convert(T,3/4),n)
+    du = fill(convert(T, 1//8), n-1)
+    d = fill(convert(T, 3//4), n)
     dl = copy(du)
     (dl,d,du)
 end
 
-function prefiltering_system{T,BC<:Union(Flat,Reflect)}(::Type{T}, n::Int, q::Quadratic{BC,OnCell})
+function prefiltering_system{T,TCoefs,BC<:Union(Flat,Reflect)}(::Type{T}, ::Type{TCoefs}, n::Int, q::Quadratic{BC,OnCell})
     dl,d,du = inner_system_diags(T,n,q)
     d[1] = d[end] = -1
     du[1] = dl[end] = 1
-    lufact!(Tridiagonal(dl, d, du)), zeros(T, n)
+    lufact!(Tridiagonal(dl, d, du)), zeros(TCoefs, n)
 end
 
-function prefiltering_system{T,BC<:Union(Flat,Reflect)}(::Type{T}, n::Int, q::Quadratic{BC,OnGrid})
+function prefiltering_system{T,TCoefs,BC<:Union(Flat,Reflect)}(::Type{T}, ::Type{TCoefs}, n::Int, q::Quadratic{BC,OnGrid})
     dl,d,du = inner_system_diags(T,n,q)
-    d[1] = d[end] = convert(T, -1)
-    du[1] = dl[end] = convert(T, 0)
+    d[1] = d[end] = -1
+    du[1] = dl[end] = 0
 
     rowspec = zeros(T,n,2)
     # first row     last row
-    rowspec[1,1] = rowspec[n,2] = convert(T, 1)
+    rowspec[1,1] = rowspec[n,2] = 1
     colspec = zeros(T,2,n)
     # third col     third-to-last col
-    colspec[1,3] = colspec[2,n-2] = convert(T, 1)
+    colspec[1,3] = colspec[2,n-2] = 1
     valspec = zeros(T,2,2)
-    # val for [1,3], val for [n,n-2]
-    valspec[1,1] = valspec[2,2] = convert(T, 1)
+    # [1,3]         [n,n-2]
+    valspec[1,1] = valspec[2,2] = 1
 
-    Woodbury(lufact!(Tridiagonal(dl, d, du)), rowspec, valspec, colspec), zeros(T, n)
+    Woodbury(lufact!(Tridiagonal(dl, d, du)), rowspec, valspec, colspec), zeros(TCoefs, n)
 end
 
-function prefiltering_system{T}(::Type{T}, n::Int, q::Quadratic{Line})
+function prefiltering_system{T,TCoefs}(::Type{T}, ::Type{TCoefs}, n::Int, q::Quadratic{Line})
     dl,d,du = inner_system_diags(T,n,q)
-    d[1] = d[end] = convert(T, 1)
-    du[1] = dl[end] = convert(T, -2)
+    d[1] = d[end] = 1
+    du[1] = dl[end] = -2
 
     rowspec = zeros(T,n,2)
     # first row     last row
-    rowspec[1,1] = rowspec[n,2] = convert(T, 1)
+    rowspec[1,1] = rowspec[n,2] = 1
     colspec = zeros(T,2,n)
     # third col     third-to-last col
-    colspec[1,3] = colspec[2,n-2] = convert(T, 1)
+    colspec[1,3] = colspec[2,n-2] = 1
     valspec = zeros(T,2,2)
-    # val for [1,3], val for [n,n-2]
-    valspec[1,1] = valspec[2,2] = convert(T, 1)
+    # [1,3]         [n,n-2]
+    valspec[1,1] = valspec[2,2] = 1
 
-    Woodbury(lufact!(Tridiagonal(dl, d, du)), rowspec, valspec, colspec), zeros(T, n)
+    Woodbury(lufact!(Tridiagonal(dl, d, du)), rowspec, valspec, colspec), zeros(TCoefs, n)
 end
 
-function prefiltering_system{T}(::Type{T}, n::Int, q::Quadratic{Free})
+function prefiltering_system{T,TCoefs}(::Type{T}, ::Type{TCoefs}, n::Int, q::Quadratic{Free})
     dl,d,du = inner_system_diags(T,n,q)
-    d[1] = d[end] = convert(T, 1)
-    du[1] = dl[end] = convert(T, -3)
+    d[1] = d[end] = 1
+    du[1] = dl[end] = -3
 
     rowspec = zeros(T,n,4)
-    rowspec[1,1] = rowspec[1,2] = rowspec[n,3] = rowspec[n,4] = convert(T,1)
+    # first row     first row       last row       last row
+    rowspec[1,1] = rowspec[1,2] = rowspec[n,3] = rowspec[n,4] = 1
     colspec = zeros(T,4,n)
-    colspec[1,3] = colspec[2,4] = colspec[3,n-2] = colspec[4,n-3] = convert(T,1)
+    # third col     fourth col     third-to-last col  fourth-to-last col
+    colspec[1,3] = colspec[2,4] = colspec[3,n-2] = colspec[4,n-3] = 1
     valspec = zeros(T,4,4)
-    valspec[1,1] = valspec[3,3] = convert(T, 3)
-    valspec[2,2] = valspec[4,4] = convert(T, -1)
+    # [1,3]          [n,n-2]
+    valspec[1,1] = valspec[3,3] = 3
+    # [1,4]          [n,n-3]
+    valspec[2,2] = valspec[4,4] = -1
 
-    Woodbury(lufact!(Tridiagonal(dl, d, du)), rowspec, valspec, colspec), zeros(T, n)
+    Woodbury(lufact!(Tridiagonal(dl, d, du)), rowspec, valspec, colspec), zeros(TCoefs, n)
 end
 
-function prefiltering_system{T}(::Type{T}, n::Int, q::Quadratic{Periodic})
+function prefiltering_system{T,TCoefs}(::Type{T}, ::Type{TCoefs}, n::Int, q::Quadratic{Periodic})
     dl,d,du = inner_system_diags(T,n,q)
 
     rowspec = zeros(T,n,2)
-    rowspec[1,1] = rowspec[n,2] = convert(T,1)
+    # first row       last row
+    rowspec[1,1] = rowspec[n,2] = 1
     colspec = zeros(T,2,n)
-    colspec[1,n] = colspec[2,1] = convert(T,1)
+    # last col         first col
+    colspec[1,n] = colspec[2,1] = 1
     valspec = zeros(T,2,2)
-    valspec[1,1] = convert(T, 1/8)
-    valspec[2,2] = convert(T, 1/8)
+    # [1,n]            [n,1]
+    valspec[1,1] = valspec[2,2] = 1//8
 
-    Woodbury(lufact!(Tridiagonal(dl, d, du)), rowspec, valspec, colspec), zeros(T, n)
+    Woodbury(lufact!(Tridiagonal(dl, d, du)), rowspec, valspec, colspec), zeros(TCoefs, n)
 end
diff --git a/test/gradient.jl b/test/gradient.jl
index 77d809db..e7f6d311 100644
--- a/test/gradient.jl
+++ b/test/gradient.jl
@@ -9,8 +9,13 @@ g1(x) = 2pi/(nx-1) * cos((x-3)*2pi/(nx-1) - 1)
 # Gradient of Constant should always be 0
 itp1 = Interpolation(Float64[f1(x) for x in 1:nx-1],
             Constant(OnGrid()), ExtrapPeriodic())
+
+g = Array(Float64, 1)
+
 for x in 1:nx
     @test gradient(itp1, x)[1] == 0
+    @test gradient!(g, itp1, x)[1] == 0
+    @test g[1] == 0
 end
 
 # Since Linear is OnGrid in the domain, check the gradients between grid points
@@ -18,6 +23,8 @@ itp1 = Interpolation(Float64[f1(x) for x in 1:nx-1],
             Linear(OnGrid()), ExtrapPeriodic())
 for x in 2.5:nx-1.5
     @test_approx_eq_eps g1(x) gradient(itp1, x)[1] abs(.1*g1(x))
+    @test_approx_eq_eps g1(x) gradient!(g, itp1, x)[1] abs(.1*g1(x))
+    @test_approx_eq_eps g1(x) g[1] abs(.1*g1(x))
 end
 
 # Since Quadratic is OnCell in the domain, check gradients at grid points
@@ -25,6 +32,8 @@ itp1 = Interpolation(Float64[f1(x) for x in 1:nx-1],
             Quadratic(Periodic(),OnGrid()), ExtrapPeriodic())
 for x in 2:nx-1
     @test_approx_eq_eps g1(x) gradient(itp1, x)[1] abs(.05*g1(x))
+    @test_approx_eq_eps g1(x) gradient!(g, itp1, x)[1] abs(.05*g1(x))
+    @test_approx_eq_eps g1(x) g[1] abs(.1*g1(x))
 end
 
 end
diff --git a/test/linear.jl b/test/linear.jl
index 503e20b2..d116bb09 100644
--- a/test/linear.jl
+++ b/test/linear.jl
@@ -42,6 +42,11 @@ xlo, xhi = itp3[.9], itp3[xmax+.2]
 @test xlo == A[1]
 @test xhi == A[end]
 
+# Rational element types
+A = Rational{Int}[x//10+2 for x in 1:10]
+itp = Interpolation(A, Linear(OnGrid()), ExtrapNaN())
+@test itp[11//10] == (11//10)//10+2
+
 end
 
 module Linear2DTests
diff --git a/test/quadratic.jl b/test/quadratic.jl
index 5ca7d094..e89452c5 100644
--- a/test/quadratic.jl
+++ b/test/quadratic.jl
@@ -12,7 +12,7 @@ A = Float64[f(x) for x in 1:xmax]
 
 itp1 = Interpolation(A, Quadratic(Flat(),OnCell()), ExtrapError())
 
-for x in [3.1:.2:4.3]  
+for x in [3.1:.2:4.3]
     @test_approx_eq_eps f(x) itp1[x] abs(.1*f(x))
 end
 
@@ -44,7 +44,7 @@ for x in [3.1:.2:4.3]
     @test_approx_eq_eps f(x) itp1[x] abs(.1*f(x))
 end
 
-xlo, xhi = itp3[.9], itp3[xmax+.2]
+xlo, xhi = itp3[.2], itp3[xmax+.7]
 @test_approx_eq xlo A[1]
 @test_approx_eq xhi A[end]
 
diff --git a/test/runtests.jl b/test/runtests.jl
index 27fd6b0e..333cb5b4 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -15,6 +15,9 @@ include("on-grid.jl")
 # test gradient evaluation
 include("gradient.jl")
 
+# test interpolation with specific types
+include("typing.jl")
+
 # Tests copied from Grid.jl's old test suite
 #include("grid.jl")
 
diff --git a/test/typing.jl b/test/typing.jl
new file mode 100644
index 00000000..f4ab3b8f
--- /dev/null
+++ b/test/typing.jl
@@ -0,0 +1,37 @@
+module TypingTests
+
+using Interpolations, Base.Test
+
+nx = 10
+f(x) = convert(Float32, x^3/(nx-1))
+g(x) = convert(Float32, 3x^2/(nx-1))
+
+A = Float32[f(x) for x in 1:nx]
+
+itp = Interpolation(A, Quadratic(Flat(), OnCell()), ExtrapConstant())
+
+# display(plot(
+#     layer(x=1:nx,y=[f(x) for x in 1:1//1:nx],Geom.point),
+#     layer(x=1:.1:nx,y=[itp[x] for x in 1:1//10:nx],Geom.path),
+# ))
+
+for x in [3.1:.2:4.3]
+    @test_approx_eq_eps float(f(x)) float(itp[x]) abs(.1*f(x))
+end
+
+@test typeof(itp[3.5f0]) == Float32
+
+for x in [3.1:.2:4.3]
+    @test_approx_eq_eps g(x) gradient(itp, x) abs(.1*g(x))
+end
+
+@test typeof(gradient(itp, 3.5)[1]) == Float32
+
+if Base.VERSION >= v"0.4-" # This doesn't work on 0.3 due to rational errors in linalg
+    # Rational element types
+    R = Rational{Int}[x^2//10 for x in 1:10]
+    itp = Interpolation(R, Quadratic(Free(),OnCell()), ExtrapNaN())
+    @test itp[11//10] == (11//10)^2//10
+end
+
+end