Simpson's Rule Implementation trick?

[Namir]

I found C code in an article on the Internet that implements Simpson's rule using a single for loop. Here is the pseudo-code for the C code:

Code:


Simpson's Rule for Numerical Integration Ver 1
======================================

Given f(x) and limit [a, b} with n points

if n modulo 2 = 0 then n = n + 1
h = (b - a) / (n + 1)
sum = f(a) + f(b) + 4*f(a+h)

for i=2 to n step 2
  sum = sum + 2*f(a+i*h) + 4*f(a + (i+1)*h)
next
area = h/3*sum

Here is a version that I modified to avoid using a+i*h and a+(i+1)*h:

Code:


Simpson's Rule for Numerical Integration Ver 2
======================================

Given f(x) and limit [a, b} with n points

if n modulo 2 = 0 then n = n + 1
h = (b - a) / (n + 1)
x = a+h
sum = f(a) + f(b) + 4*f(x)

for i=2 to n step 2
  x = x + h
  sum = sum + 2*f(x)
  x = x + h
  sum = sum + 4*f(x)
next
area = h/3*sum

Consolidating the for loop to have a single update to sum:

Code:

Simpson's Rule for Numerical Integration Ver 3
==============================================

Given f(x) and limit [a, b} with n points

if n modulo 2 = 0 then n = n + 1
h = (b - a) / (n + 1)
x = a+h
sum = f(a) + f(b) + 4*f(x)
c = 2
for i=2 to n
  x = x + h
  sum = sum + c*f(x)
  c = 6 - c
next
area = h/3*sum

Enjoy!

Namir