Some horribly unoptimized quadrature rules, integrating the "black swan" function from Eric Lippert's "Fixing Random" series.

quadratures.cs

using System;
using System.Collections.Generic;
using System.Linq;

namespace Quadratures
{
    interface Integrator
    {
        double Integrate(Func<double, double> f, double a, double b);
    }

    public static class Functions
    {
        public static double WeirdFunction(double x)
        {
            return Math.Atan(1000 * (x - 0.45)) * 20 - 31.2;
        }

        static readonly double INVERSE_ROOT_PI = 1.0 / Math.Sqrt(2 * Math.PI);

        public static double NormalDistributionDensity(double x, double mu, double sigma)
        {
            return Math.Exp(-(x - mu) * (x - mu) / (2 * sigma * sigma)) * INVERSE_ROOT_PI / sigma;
        }
    }

    struct PointWeight
    {
        public double Point { get; }
        public double Weight { get; }

        public PointWeight(double point, double weight)
        {
            Point = point;
            Weight = weight;
        }
    }

    class GaussianQuadratureIntegrator : Integrator
    {
        static readonly PointWeight[] TwoPointsRuleData =
        {
            new PointWeight(-Math.Sqrt(1.0 / 3.0), 1.0),
            new PointWeight( Math.Sqrt(1.0 / 3.0), 1.0),
        };

        static readonly PointWeight[] ThreePointsRuleData =
        {
            new PointWeight(-Math.Sqrt(3.0 / 5.0), 5.0 / 9.0),
            new PointWeight(                  0.0, 8.0 / 9.0),
            new PointWeight( Math.Sqrt(3.0 / 5.0), 5.0 / 9.0),
        };

        static readonly PointWeight[] FivePointsRuleData =
        {
            new PointWeight(-0.9061798459386640, 0.2369268850561891),
            new PointWeight(-0.5384693101056831, 0.4786286704993665),
            new PointWeight(                0.0, 0.5688888888888889),
            new PointWeight( 0.5384693101056831, 0.4786286704993665),
            new PointWeight( 0.9061798459386640, 0.2369268850561891),
        };

        static readonly PointWeight[] SevenPointsRuleData =
        {
            new PointWeight(-0.9491079123427585, 0.1294849661688697),
            new PointWeight(-0.7415311855993944, 0.2797053914892767),
            new PointWeight(-0.4058451513773972, 0.3818300505051189),
            new PointWeight(                0.0, 0.4179591836734694),
            new PointWeight( 0.4058451513773972, 0.3818300505051189),
            new PointWeight( 0.7415311855993944, 0.2797053914892767),
            new PointWeight( 0.9491079123427595, 0.1294849661688697),
        };

        public static readonly IReadOnlyCollection<int> ValidPointCounts = new List<int> { 2, 3, 5, 7 }.AsReadOnly();

        PointWeight[] RuleData { get; }

        public GaussianQuadratureIntegrator(int pointsCount)
        {
            switch (pointsCount)
            {
                case 2: RuleData = TwoPointsRuleData; break;
                case 3: RuleData = ThreePointsRuleData; break;
                case 5: RuleData = FivePointsRuleData; break;
                case 7: RuleData = SevenPointsRuleData; break;
                default: throw new ArgumentOutOfRangeException(nameof(pointsCount));
            }
        }

        public double Integrate(Func<double, double> f, double a, double b)
        {
            double center = (a + b) / 2;
            double halfLength = (b - a) / 2;

            var result = 0.0;

            for (int i = 0; i < RuleData.Length; i++)
            {
                result += RuleData[i].Weight * f(center + halfLength * RuleData[i].Point);
            }
            result *= halfLength;

            return result;
        }
    }

    class DumbRiemannSummator : Integrator
    {
        int NumberOfSegments { get; }

        public DumbRiemannSummator(int numberOfSegments)
        {
            NumberOfSegments = numberOfSegments;
        }

        public double Integrate(Func<double, double> f, double a, double b)
        {
            var result = 0.0;

            for (int i = 0; i < NumberOfSegments; i++)
            {
                result += f(a + i * (b - a) / NumberOfSegments);
            }
            result /= NumberOfSegments;

            return result;
        }
    }

    class AdaptiveIntegrator : Integrator
    {
        Integrator PreciseIntegrator { get; }
        Integrator CrudeIntegrator { get; }
        double MaxError { get; }
        int MaxSubsegments { get; }

        public int SubsegmentCount { get; private set; }
        public double ErrorEstimate { get; private set; }

        public AdaptiveIntegrator(Integrator preciseIntegrator, Integrator crudeIntegrator, double maxError, int maxSubsegments)
        {
            PreciseIntegrator = preciseIntegrator;
            CrudeIntegrator = crudeIntegrator;
            MaxError = maxError;
            MaxSubsegments = maxSubsegments;

            SubsegmentCount = 0;
            ErrorEstimate = double.NaN;
        }

        struct EstimationResult {
            public double Result { get; }
            public double Error { get; }

            public EstimationResult(double result, double error)
            {
                Result = result;
                Error = error;
            }
        }

        struct Subsegment
        {
            public double Left { get; }
            public double Right { get; }

            public Subsegment(double left, double right)
            {
                Left = left;
                Right = right;
            }
        }

        EstimationResult EstimateIntegral(Func<double, double> f, Subsegment s)
        {
            var result = PreciseIntegrator.Integrate(f, s.Left, s.Right);
            var crudeResult = CrudeIntegrator.Integrate(f, s.Left, s.Right);
            var error = Math.Abs(result - crudeResult);

            return new EstimationResult(result, error);
        }

        static int FindIndexOfMaxError(List<EstimationResult> estimations)
        {
            var maxError = double.NegativeInfinity;
            var result = -1;

            for (var index = 0; index < estimations.Count; index++)
            {
                if (maxError < estimations[index].Error)
                {
                    result = index;
                    maxError = estimations[index].Error;
                }
            }

            return result;
        }

        public double Integrate(Func<double, double> f, double a, double b)
        {
            var subsegment = new Subsegment(a, b);
            var estimation = EstimateIntegral(f, subsegment);

            SubsegmentCount = 1;
            ErrorEstimate = estimation.Error;

            var subsegments = new List<Subsegment> { subsegment };
            var estimations = new List<EstimationResult> { estimation };

            if (estimation.Error < MaxError)
            {
                return estimation.Result;
            }

            // Loop invariant: there are exactly SubsegmentCount subsegments and estimations.
            // Loop exit condition: there are exactly MaxSubsegments subsegments and estimations.
            while (SubsegmentCount < MaxSubsegments)
            {
                var maxErrorIndex = FindIndexOfMaxError(estimations);
                var dividedSubsegment = subsegments[maxErrorIndex];

                var subsegmentCenter = (dividedSubsegment.Left + dividedSubsegment.Right) / 2;

                var leftSubsegment = new Subsegment(dividedSubsegment.Left, subsegmentCenter);
                var rightSubsegment = new Subsegment(subsegmentCenter, dividedSubsegment.Right);

                var leftEstimation = EstimateIntegral(f, leftSubsegment);
                var rightEstimation = EstimateIntegral(f, rightSubsegment);

                subsegments[maxErrorIndex] = leftSubsegment;
                estimations[maxErrorIndex] = leftEstimation;

                subsegments.Add(rightSubsegment);
                estimations.Add(rightEstimation);

                var totalError = estimations.Sum(x => x.Error);

                SubsegmentCount++;
                ErrorEstimate = totalError;

                if (totalError < MaxError)
                {
                    var totalResult = estimations.Sum(x => x.Result);
                    return totalResult;
                }
            }

            throw new Exception("Exceeded subsegments limit");
        }
    }

    class InvocationCounter<TArg, TResult>
    {
        private Func<TArg, TResult> Func { get; }

        public int Counter { get; private set; }

        public InvocationCounter(Func<TArg, TResult> func)
        {
            Func = func;
            Counter = 0;
        }

        public TResult Invoke(TArg arg)
        {
            Counter++;
            return Func(arg);
        }
    }

    class Program
    {
        static double F(double x)
        {
            return Functions.NormalDistributionDensity(x, 0.75, 0.09) * Functions.WeirdFunction(x);
        }

        static void TryIntegrator(string name, Integrator integrator)
        {
            const double a = 0.0;
            const double b = 1.0;

            var invoker = new InvocationCounter<double, double>(F);
            var result = integrator.Integrate(invoker.Invoke, a, b);
            Console.Write("Method '{0}': result = {1}, after {2} evaluations of F", name, result, invoker.Counter);

            var adaptiveIntegrator = integrator as AdaptiveIntegrator;
            if (adaptiveIntegrator != null)
            {
                Console.Write(", with {0} subsegments, error estimation is {1}", 
                    adaptiveIntegrator.SubsegmentCount, adaptiveIntegrator.ErrorEstimate);
            }

            Console.WriteLine();
        }

        static void Main(string[] args)
        {
            const int RIEMANN_SEGMENTS_COUNT = 1000;
            const double MAX_ERROR = 1e-5;
            const int MAX_SUBSEGMENTS = 50;

            TryIntegrator($"Dumb Riemann summator ({RIEMANN_SEGMENTS_COUNT})", new DumbRiemannSummator(RIEMANN_SEGMENTS_COUNT));
            foreach (var pointCount in GaussianQuadratureIntegrator.ValidPointCounts)
            {
                TryIntegrator($"{pointCount}-point Gaussian quadrature", new GaussianQuadratureIntegrator(pointCount));
            }

            TryIntegrator($"Adaptive Gaussian quadrature with 5- and 7-points rules",
                new AdaptiveIntegrator(
                    new GaussianQuadratureIntegrator(7),
                    new GaussianQuadratureIntegrator(5),
                    MAX_ERROR, MAX_SUBSEGMENTS));

            Console.ReadLine();
        }
    }
}