@@ -17,40 +17,38 @@ immutable TrapezoidalEvenFast <: IntegrationMethod end
17
17
immutable SimpsonEven <: IntegrationMethod end # https://en.wikipedia.org/wiki/Simpson%27s_rule#Alternative_extended_Simpson.27s_rule
18
18
immutable SimpsonEvenFast <: IntegrationMethod end
19
19
20
- integrate {T<:AbstractFloat} (x:: Vector{T} , y:: Vector{T} , method:: IntegrationMethod = Trapezoidal ()) = integrate (x, y, method)
21
-
22
- function integrate {T<:AbstractFloat} (x:: Vector{T} , y:: Vector{T} , :: Trapezoidal )
23
- @assert length (x) == length (y) " x and y vectors must be of the same length!"
24
- retval = zero (eltype (x))
20
+ function integrate {X<:Real, Y<:Real} (x:: Vector{X} , y:: Vector{Y} , :: Trapezoidal )
21
+ @assert length (x) == length (y) " x and y vectors must be of the same length!"
22
+ retval = zero (promote (x[1 ], y[1 ])[1 ])
25
23
for i in 1 : length (y)- 1
26
- retval += (x[i+ 1 ] - x[i]) * (y[i] + y[i+ 1 ])
24
+ retval += (x[i+ 1 ] - x[i]) * (y[i] + y[i+ 1 ])
27
25
end
28
- return 0.5 * retval
26
+ return 0.5 * retval
29
27
end
30
28
31
- function integrate {T<:AbstractFloat } (x:: Vector{T } , y:: Vector{T } , :: TrapezoidalEven )
29
+ function integrate {X<:Real, Y<:Real } (x:: Vector{X } , y:: Vector{Y } , :: TrapezoidalEven )
32
30
@assert length (x) == length (y) " x and y vectors must be of the same length!"
33
31
return 0.5 * (x[end ] - x[1 ]) / (length (y) - 1 ) * (y[1 ] + y[end ] + sum (y[2 : end - 1 ]))
34
32
end
35
33
36
- function integrate {T<:AbstractFloat } (x:: Vector{T } , y:: Vector{T } , :: TrapezoidalFast )
37
- retval = zero (eltype (x) )
34
+ function integrate {X<:Real, Y<:Real } (x:: Vector{X } , y:: Vector{Y } , :: TrapezoidalFast )
35
+ retval = zero (promote (x[ 1 ], y[ 1 ])[ 1 ] )
38
36
@fastmath @simd for i in 1 : length (y)- 1
39
37
@inbounds retval += (x[i+ 1 ] - x[i]) * (y[i] + y[i+ 1 ])
40
38
end
41
39
return 0.5 * retval
42
40
end
43
41
44
- function integrate {T<:AbstractFloat } (x:: Vector{T } , y:: Vector{T } , :: TrapezoidalEvenFast )
45
- retval = zero (eltype (x) )
42
+ function integrate {X<:Real, Y<:Real } (x:: Vector{X } , y:: Vector{Y } , :: TrapezoidalEvenFast )
43
+ retval = zero (promote (x[ 1 ], y[ 1 ])[ 1 ] )
46
44
N = length (y) - 1
47
45
@fastmath @simd for i in 2 : N
48
46
@inbounds retval += y[i]
49
47
end
50
48
@inbounds return (x[end ] - x[1 ]) / N * (retval + 0.5 * y[1 ] + 0.5 * y[end ])
51
49
end
52
50
53
- function integrate {T<:AbstractFloat } (x:: Vector{T } , y:: Vector{T } , :: SimpsonEven )
51
+ function integrate {X<:Real, Y<:Real } (x:: Vector{X } , y:: Vector{Y } , :: SimpsonEven )
54
52
@assert length (x) == length (y) " x and y vectors must be of the same length!"
55
53
retval = (17 * y[1 ] + 59 * y[2 ] + 43 * y[3 ] + 49 * y[4 ] + 49 * y[end - 3 ] + 43 * y[end - 2 ] + 59 * y[end - 1 ] + 17 * y[end ]) / 48
56
54
for i in 5 : length (y) - 1
@@ -59,7 +57,7 @@ function integrate{T<:AbstractFloat}(x::Vector{T}, y::Vector{T}, ::SimpsonEven)
59
57
return (x[end ] - x[1 ]) / (length (y) - 1 ) * retval
60
58
end
61
59
62
- function integrate {T<:AbstractFloat } (x:: Vector{T } , y:: Vector{T } , :: SimpsonEvenFast )
60
+ function integrate {X<:Real, Y<:Real } (x:: Vector{X } , y:: Vector{Y } , :: SimpsonEvenFast )
63
61
@inbounds retval = 17 * y[1 ] + 59 * y[2 ] + 43 * y[3 ] + 49 * y[4 ]
64
62
@inbounds retval += 49 * y[end - 3 ] + 43 * y[end - 2 ] + 59 * y[end - 1 ] + 17 * y[end ]
65
63
retval /= 48
@@ -69,4 +67,6 @@ function integrate{T<:AbstractFloat}(x::Vector{T}, y::Vector{T}, ::SimpsonEvenFa
69
67
@inbounds return (x[end ] - x[1 ]) / (length (y) - 1 ) * retval
70
68
end
71
69
70
+ integrate {X<:Real, Y<:Real} (x:: Vector{X} , y:: Vector{Y} ) = integrate (x, y, TrapezoidalFast ())
71
+
72
72
end
0 commit comments