double integrate_error(
    double f( double x ),
    double g( double x ),
    double a,
    double b,
    unsigned int n
) {
    assert( b > a );
    h = (b - a)/n;
    double integral{ 0.0 };

    for ( unsigned int k{0}; k < n; ++k ) {
        double x0{a + k*h};
        double x1{a + (k + 0.5)*h};
        double x2{a + (k + 1.0)*h};

        double d0{ f( x0 ) - g( x0 ) };
        double d1{ f( x1 ) - g( x1 ) };
        double d2{ f( x2 ) - g( x2 ) };

        double disc{ d0*d0 - 8.0*d0*d1 - 2.0*d0*d2
                          + 16.0*d1*d1 - 8.0*d1*d2
                                           + d2*d2 };

        if ( disc < 0.0 ) {
            if ( d0 > 0.0 ) {
                integral += d0 + 4.0*d1 + d2;
            } else {
                integral -= d0 + 4.0*d1 + d2;
            }
        } else {
            double a1{ (4.0*(-d0 + 2.0*d1 - d2)*x0 - h*(3.0*d0 - 4.0*d1 + d2)) };
            double a2{ (2*d0 - 4*d1 + 2*d2) };

            disc = std::sqrt( disc );

            double r_min{ (-a1 - disc)/2.0*a2 };
            double r_max{ (-a1 + disc)/2.0*a2 };

            if ( r_min > r_max ) {
                double t{ r_min };
                r_min = r_max;
                r_max = t;
            }

            if ( (r_max <= x0) || (r_min >= x2) ) {
                // -o---o--|---------------|-------
                // --------|---------------|-o--o--
                //         x0              x2
                if ( d0 > 0.0 ) {
                    integral += d0 + 4.0*d1 + d2;
                } else {
                    integral -= d0 + 4.0*d1 + d2;
                }
            } else if ( (r_min <= x0) || (r_max >= x2) ) {
                // ---o----|---------------|-o-----
                //         x0              x2
                if ( d0 > 0.0 ) {
                    integral -= d0 + 4.0*d1 + d2;
                } else {
                    integral += d0 + 4.0*d1 + d2;
                }
            } else if ( (r_min < x0) && (r_max > x0) && (r_max < x2) ) {
                // -o------|---------o-----|-------
                //         x0              x2
            } else if ( (r_min > x0) && (r_max > x2) ) {
                // --------|-----o---o-----|-------
                //         x0              x2
            } else if ( (r_min > x0) && (r_min < x2) && (r_max > x2) ) {
                // --------|-----o---------|---o---
                //         x0              x2
            }

    }

    return integral*h/3.0;
}