Simpson 3/8 pseudo-code

Simpson's 3/8 rule needs the number of sub-intervals n to be a multiple of 3. With h = (b − a)/n:

I ≈ (3h/8)[y₀ + yₙ + 3(y₁ + y₂ + y₄ + y₅ + …) + 2(y₃ + y₆ + y₉ + …)]

So every interior point whose index is a multiple of 3 gets weight 2; all other interior points get weight 3.

Pseudo-code

START
  DEFINE function f(x)
  READ a, b, n
  IF n mod 3 ≠ 0 THEN
      PRINT "n must be a multiple of 3"
      STOP
  ENDIF
  h   ← (b − a) / n
  sum ← f(a) + f(b)
  FOR i ← 1 TO n − 1
      x ← a + i·h
      IF i mod 3 = 0 THEN
          sum ← sum + 2·f(x)
      ELSE
          sum ← sum + 3·f(x)
      ENDIF
  ENDFOR
  I ← 3·h·sum / 8
  PRINT I
STOP

C implementation

#include <stdio.h>
#include <math.h>

double f(double x)
{
    return 1.0 / (1.0 + x * x);          /* function to integrate */
}

int main(void)
{
    double a, b, h, sum, x;
    int n, i;

    printf("Enter a, b and n (n multiple of 3): ");
    scanf("%lf %lf %d", &a, &b, &n);
    if (n % 3 != 0)
    {
        printf("n must be a multiple of 3\n");
        return 1;
    }
    h = (b - a) / n;
    sum = f(a) + f(b);
    for (i = 1; i < n; i++)
    {
        x = a + i * h;
        if (i % 3 == 0) sum += 2 * f(x);
        else            sum += 3 * f(x);
    }
    printf("Integral = %.4f\n", 3 * h / 8 * sum);
    return 0;
}

Sample run: for f = 1/(1+x²), a = 0, b = 1, n = 6 the program prints Integral = 0.7854 (exact π/4 = 0.7854).