replace 0.5 by 1//2 (#10)

Author: jw3126

src/NumericalIntegration.jl

@@ -7,6 +7,7 @@ using Compat
 export integrate
 export Trapezoidal, TrapezoidalEven, TrapezoidalFast, TrapezoidalEvenFast
 export SimpsonEven, SimpsonEvenFast
+export IntegrationMethod
 
 @compat abstract type IntegrationMethod end
 
@@ -17,26 +18,28 @@ immutable TrapezoidalEvenFast <: IntegrationMethod end
 immutable SimpsonEven         <: IntegrationMethod end # https://en.wikipedia.org/wiki/Simpson%27s_rule#Alternative_extended_Simpson.27s_rule
 immutable SimpsonEvenFast     <: IntegrationMethod end
 
+const HALF = 1//2
+
 function integrate{X<:Number, Y<:Number}(x::AbstractVector{X}, y::AbstractVector{Y}, ::Trapezoidal)
   @assert length(x) == length(y) "x and y vectors must be of the same length!"
   retval = zero(promote(x[1], y[1])[1])
     for i in 1 : length(y)-1
       retval += (x[i+1] - x[i]) * (y[i] + y[i+1])
     end
-  return 0.5 * retval
+  return HALF * retval
 end
 
 function integrate{X<:Number, Y<:Number}(x::AbstractVector{X}, y::AbstractVector{Y}, ::TrapezoidalEven)
     @assert length(x) == length(y) "x and y vectors must be of the same length!"
-    return (x[2] - x[1]) * (0.5 * (y[1] + y[end]) + sum(y[2:end-1]))
+    return (x[2] - x[1]) * (HALF * (y[1] + y[end]) + sum(y[2:end-1]))
 end
 
 function integrate{X<:Number, Y<:Number}(x::AbstractVector{X}, y::AbstractVector{Y}, ::TrapezoidalFast)
     retval = zero(promote(x[1], y[1])[1])
     @fastmath @simd for i in 1 : length(y)-1
         @inbounds retval += (x[i+1] - x[i]) * (y[i] + y[i+1])
     end
-    return 0.5 * retval
+    return HALF * retval
 end
 
 function integrate{X<:Number, Y<:Number}(x::AbstractVector{X}, y::AbstractVector{Y}, ::TrapezoidalEvenFast)
@@ -45,7 +48,7 @@ function integrate{X<:Number, Y<:Number}(x::AbstractVector{X}, y::AbstractVector
     @fastmath @simd for i in 2 : N
         @inbounds retval += y[i]
     end
-    @inbounds return (x[2] - x[1]) * (retval + 0.5*y[1] + 0.5*y[end])
+    @inbounds return (x[2] - x[1]) * (retval + HALF*y[1] + HALF*y[end])
 end
 
 function integrate{X<:Number, Y<:Number}(x::AbstractVector{X}, y::AbstractVector{Y}, ::SimpsonEven)

test/runtests.jl

@@ -1,13 +1,20 @@
 using NumericalIntegration
 using Base.Test
 
-methods = [Trapezoidal(), TrapezoidalEven(), TrapezoidalFast(), TrapezoidalEvenFast(), SimpsonEven(), SimpsonEvenFast()]
 x = collect(-π : π/1000 : π)
 y = sin.(x)
 p = collect(0 : π/1000 : π)
 q = sin.(p)
-for method in methods
-    println(string("Testing method: ", typeof(method)))
-    @test abs(integrate(x, y, method)) < 1e-4
-    @test abs(integrate(p, q, method)-2.0) < 1e-4
+for M in subtypes(IntegrationMethod)
+    println("Testing method: ", M)
+    for T in [Float32, Float64, BigFloat]
+        for (xs,ys,val,atol) in [
+                                 (x,y,0,1e-4),
+                                 (p,q,2,1e-4),
+                              ]
+            result = @inferred integrate(T.(xs), T.(ys),M())
+            @test isapprox(result, val, atol=atol)
+            @test typeof(result) == T
+        end
+    end
 end