function result = fint(integrand, a, b, N)

% FFT integration of a periodic function
% a, b are the limits of integration.
% N is the number of abcissas.

% First 2 lines same as midpoint1.m
dx = (b - a)/N;
x = a + dx*( (1:N) - 0.5 );

f = integrand(x);

% The following is an inefficient way of summing the values
% of f at all the abcissas. But it has very low roundoff error.
F = fft(f);
result = F(1)*dx;