C言語で学ぶシンプソンの公式による積分の近似計算

#シンプソンの公式とは
台形公式(https://qiita.com/stmtk/items/0d43372950fe3ee674be) では積分したい関数を何分割かして区間ごとに台形の面積を求めてそれらを足すと積分した値と近くなるという方法でした。
しかし、台形公式のページの絵を見ていただくとわかると思いますが、区間ごとに1次関数に近似して積分してるのと同じですね。1次関数で近似するとカーブを表現できない分誤差が出てしまいます。なので区間ごとに2次関数で近似するためのものがシンプソンの公式です。
$ f(x) $が3次以下の1変数多項式の場合

\int_a^b f(x) dx = \frac{b - a}{6}\left(f(a) + 4 f\left(\frac{a + b}{2}\right) + f(b)\right)

導き方

まず、ラグランジュの補間多項式(https://qiita.com/stmtk/items/e622bfb1e65b76d33fa1) を使って$ \left(a, f(a)\right), \left(\frac{a + b}{2}, f(\frac{a + b}{2})\right) , \left(b, f(b)\right) $の三つの点を通る2次関数を求めます。そして積分すると、シンプソンの公式と一致するはずです。証明は書くのがめんどくさいので省きます。

#ソースコード

simpson.c

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

#define EPS (1.0e-6)
#define f(x) ((x)*(x))

int main(){
    double a = 0, b = 2, sum = 0;
    for(double x = a; x < b; x += EPS * 2){
        sum += (f(x) + 4 * f(x + EPS) + f(x + EPS * 2));
    }
    sum *= (b - a) * (EPS / 6);

    printf("%f\n", sum);
}