# Integrator 2.0: A Simple Integrator in C++17

This is a follow up of A simple definite integrator class of a single variable in C++

I took most of the advice from the user: emma-x and some from user: sudo-rm-rf-slash

Here is my fully revised code:

Integrator.h

#pragma once

#include <type_traits>

template <typename Field>
struct Limits {
Field lower;
Field upper;

constexpr Limits(Field a = 0, Field b = 0) :
lower{ a },
upper{ b }
{}
};

template <typename LimitType, typename Func>
class Integrator {
private:
Limits<LimitType> limits_;
size_t step_size_;
Func integrand_;

public:
Integrator(Limits<LimitType> limits, size_t stepSize, Func integrand) :
limits_{ limits },
step_size_{ stepSize },
integrand_{ integrand }
{}

Limits<LimitType> limits() const { return limits_; }
Func* integrand() { return &integrand_; }

// This is always a 1st order integration!
constexpr auto evaluate() {
if (limits_.lower == limits_.upper) {
return 0.0;
} else if (limits_.lower > limits_.upper) {
limits_ = Limits{ limits_.upper, limits_.lower };
auto distance = limits_.upper - limits_.lower;
auto dx = distance / step_size_;
return -calculate(dx);
} else {
auto distance = limits_.upper - limits_.lower;
auto dx = distance / step_size_;
return calculate(dx);
}
}

// This will perform a second order of integration where the inner limits are defined
// by [lower, y] where "upper" is not used directly. This may be expanded in the future...
template<typename ValueType>
constexpr ValueType doubleIntegral(Limits<LimitType> innerLimits) {
// Currently the inner upper bound is not being used, if it is <= lower bound,
// it will default to lower bound + 1, just to ensure that an integration is
// is calculated. Eventually this may be expanded to return true values
// such as 0 when the lower and upper bounds are equal or the (-)integration if a > b.
if (innerLimits.upper <= innerLimits.lower) innerLimits.upper = innerLimits.lower + static_cast<LimitType>(1);

ValueType sum = 0;
ValueType dy = (ValueType)(limits_.upper - limits_.lower) / step_size_;
for (size_t i = 0; i < step_size_; i++) {
ValueType yi = (ValueType)(limits_.lower + i * dy);
ValueType dx = (yi - innerLimits.lower) / step_size_;
Integrator inner{ innerLimits, step_size_, integrand_ };
sum += inner.calculate(dx) * dy;
}
return sum;
}

private:
template<typename ValueType>
constexpr auto calculate(ValueType dx) ->
std::enable_if_t<std::is_invocable_v<Func&, ValueType&>, ValueType> {
ValueType result = 0.0;
for (size_t i = 0; i < step_size_; ++i) {
auto dy = integrand_(limits_.lower + i * dx);
auto area = dy * dx;
result += area;
}
return result;
}
};


main.cpp

#include <iostream>
#include <exception>
#include <cmath>

#include "Integrator.h"

constexpr double PI = 3.14159265358979;

double funcA0(double x) {
return x;
}

template <typename T>
constexpr T funcA(T x) {
return x;
}

template <typename T>
constexpr T funcB(T x) {
return x * x;
}

template <typename T>
T funcC(T x) {
return sin(x);
}

template <typename T>
constexpr T funcD(T x) {
return (0.5*(x*x) + (3 * x) - (1 / x));
}

double funcE(double x) {
return cos(x);
}

int main() {
try {

std::cout << "Integration of f(x) = x from a=3.0 to b=3.0 with an expected output of 0\n";
std::cout << "Testing non template function object\n";
Integrator integratorA0A{ Limits{3.0, 3.0}, 10000, funcA0 };
std::cout << integratorA0A.evaluate() << '\n';

std::cout << "Testing non template function object by reference\n";
Integrator integratorA0B{ Limits{3.0, 3.0}, 10000, &funcA0 };
std::cout << integratorA0B.evaluate() << '\n';

std::cout << "Testing function template object\n";
Integrator integratorA0C{ Limits{3.0, 3.0}, 10000, funcA<double> };
std::cout << integratorA0C.evaluate() << '\n';

std::cout << "Testing function template object by reference\n";
Integrator integratorA0D{ Limits{3.0, 3.0}, 10000, &funcA<double> };
std::cout << integratorA0D.evaluate() << '\n';

std::cout << "\nIntegration of f(x) = x from a=3.0 to b=5.0 with an expected output of 8\n";
std::cout << "Testing non template function object\n";
Integrator integratorA1A{ Limits{3.0, 5.0}, 10000, funcA0 };
std::cout << integratorA1A.evaluate() << '\n';

std::cout << "Testing non template function object by reference\n";
Integrator integratorA1B{ Limits{3.0, 5.0}, 10000, &funcA0 };
std::cout << integratorA1B.evaluate() << '\n';

std::cout << "Testing function template object\n";
Integrator integratorA1C{ Limits{3.0, 5.0}, 10000, funcA<double> };
std::cout << integratorA1C.evaluate() << '\n';

std::cout << "Testing function template object by reference\n";
Integrator integratorA1D{ Limits{3.0, 5.0}, 10000, &funcA<double> };
std::cout << integratorA1D.evaluate() << '\n';

std::cout << "\nIntegration of f(x) = x from a=5.0 to b=3.0 with an expected output of -8\n";
std::cout << "Testing non template function object\n";
Integrator integratorA2A{ Limits{5.0, 3.0}, 10000, funcA0 };
std::cout << integratorA2A.evaluate() << '\n';

std::cout << "Testing non template function object by reference\n";
Integrator integratorA2B{ Limits{5.0, 3.0}, 10000, &funcA0 };
std::cout << integratorA2B.evaluate() << '\n';

std::cout << "Testing function template object\n";
Integrator integratorA2C{ Limits{5.0, 3.0}, 10000, funcA<double> };
std::cout << integratorA2C.evaluate() << '\n';

std::cout << "Testing function template object by reference\n";
Integrator integratorA2D{ Limits{5.0, 3.0}, 10000, &funcA<double> };
std::cout << integratorA2D.evaluate() << '\n';

std::cout << "\n\nIntegration of f(x) = x^2 from a=2.0 to b=20.0 with an expected output of 2664\n";
Integrator integratorB{ Limits{2.0, 20.0}, 10000, &funcB<double> };
std::cout << integratorB.evaluate() << '\n';

std::cout << "\n\nIntegration of f(x) = sin(x) from a=0.0 to b=" << PI << " with an expected output of 2\n";
Integrator integratorC{ Limits{0.0, PI}, 10000, &funcC<double> };
std::cout << integratorC.evaluate() << '\n';

std::cout << "\n\nIntegration of f(x) = (1\\2)x^2 + 3x - (1\\x) from a=1.0 to b=10.0 \nwith an expected output of 312.6974\n";
Integrator integratorD{ Limits{1.0, 10.0 }, 10000, &funcD<double> };
std::cout << integratorD.evaluate() << '\n';

std::cout << "\n\nTesting Double Integration of f(x) = (1\\2)x^2 + 3x - (1\\x) from [3,5] and [1,y]\nwith an expected output of 65.582\n";
Integrator integrator2D{ Limits{3, 5}, 10000, &funcD<double> };
std::cout << integrator2D.doubleIntegral<double>(1.0) << '\n';

std::cout << "\n\nTesting Double Integration of f(x) = cos(x) from [0," << PI / 2 << "] and [0.25,y]\nwith an expected output of 0.61137\n";
Integrator integratorE{ Limits{0.0, PI / 2}, 10000, &funcE };
std::cout << integratorE.doubleIntegral<double>(0.25) << '\n';

}
catch (const std::exception& e) {
std::cerr << e.what() << std::endl;
return EXIT_FAILURE;
}

return EXIT_SUCCESS;
}


Output

Integration of f(x) = x from a=3.0 to b=3.0 with an expected output of 0
Testing non template function object
0
Testing non template function object by reference
0
Testing function template object
0
Testing function template object by reference
0

Integration of f(x) = x from a=3.0 to b=5.0 with an expected output of 8
Testing non template function object
7.9998
Testing non template function object by reference
7.9998
Testing function template object
7.9998
Testing function template object by reference
7.9998

Integration of f(x) = x from a=5.0 to b=3.0 with an expected output of -8
Testing non template function object
-7.9998
Testing non template function object by reference
-7.9998
Testing function template object
-7.9998
Testing function template object by reference
-7.9998

Integration of f(x) = x^2 from a=2.0 to b=20.0 with an expected output of 2664
2663.64

Integration of f(x) = sin(x) from a=0.0 to b=3.14159 with an expected output of 2
2

Integration of f(x) = (1\2)x^2 + 3x - (1\x) from a=1.0 to b=10.0
with an expected output of 312.6974
312.663

Testing Double Integration of f(x) = (1\2)x^2 + 3x - (1\x) from [3,5] and [1,y]
with an expected output of 65.582
65.5725

Testing Double Integration of f(x) = cos(x) from [0,1.5708] and [0.25,y]
with an expected output of 0.61137
0.611325


List of Fixed Issues:

• I removed the default destructor
• I removed the unnecessary member variables and their equivalent functions.
• I'm now caching them on the stack.
• I removed the member function from the Limits struct.
• I also removed the swap from the limits.
• if a==b then the integral is 0 and
• if a>b then the integral is (-) the integral from b to a.
• These checks for the limits of integration I moved into the Integrator class.
• I tried to make sure everything follows RAII.
• The class now works in a constexpr context provided the function objects are constexpr.
• I moved the loop calculations for the approximation from the evaluate function and moved it in its own private member function.
• I'm using {} initialization for all object instantiations.
• I templated both classes.
• I can also pass in a function object that is also templated.

All of the changes to my classes are all based on the advice of the user Emma X. I did take some of the advice from user sudo-rm-rf-slash in that I read up on the links that they had provided. As for changing the structure of the code from a class to be just a set of functions I did not. I have my own reasons why I want this to be a class object. I plan on expanding this later for doing different types of integrations and I will need for this class to keep track of its integrand and its limits of integration.

Yes, I do understand that the current implementation is based on Riemann Sums. This is the intentional base structure of my class since it is the simplest form in approximate integrals. Eventually, I plan on expanding this to incorporate different methods and algorithms. These will be in the form of member functions naming them according to their algorithm.

Other changes that were made:

• I renamed integrate to doubleIntegral
• I also changed its signature to take a Limits struct instead of 2 individual parameters.
• I also removed the local innerLimits variable since it is now being passed in, and just updated it accordingly if needed.
• In both the evaluate and doubleIntegral function, I have early returns if the limits are equal.
• I only have a check to see if a > b for the 1st order integration to return the (-)integral. I did not apply this behavior to the doubleIntegral function yet.
• I removed a lot of unnecessary local variables in my functions and inlined many of the calculations.
• Finally I added a sanity check to make sure that the template argument that is passed into the construct is an invokable object.
• The last of the changes, I modified my driver program slightly to show different cases of the same function being used in different contexts and to show the results when the limits are equal, less than and greater than... and I added a second example of a double integration using a trig function.

Where Do I go from here?

• I'd like to know how to incorporate the innerLimit's upper into the calculations when calculating doubleIntegral.
• I'd like to know if there is any room for improvements, efficiency, readability, portability, etc.
• Are there still any existing code smells?
• -Note- I mentioned this in the previous post, I'm not concerned with namespace, that can be done at any time!
• User: Emma X had mentioned the use of [[nodiscard]] however, I'm not familiar with its use... I've never really used that feature of C++ before. It is new to me.
• I have not yet tried to use lambdas within this code yet as that is on my todo list.
• I would also like to be able to add in the ability or the concept to incorporate infinity as the limits of integration and since this is a numerical type integrator, only the upper or lower bound can go to either + or - infinity, but not both. I'm not at the stage of returning back a family of functions with a constant of integration... That's a little beyond my skills at the moment since that would involved a function lookup table and an extensive parser.

Let me know what you think, I look forward to your advice, constructive criticism, opinions, positions, and all feedback is welcomed!

• How do you argue that initializing the loop variable before the for-loop helps with optimization? It compiles to the same assembly, but your way is worse in terms of readability and scoping. Jun 3 '20 at 21:26
• @EmmaX Well, it may not make a difference these days, compilers have improved quite a bit, but I only have to declare and initialize it once outside the loop, then just iterate through the loop. This is an old trick going back to when I first started learning the C/C++ languages... These days, the compilers will probably optimize it anyway... Then again this trick is commonly used for extreme efficiency when working with "nested loops" and it usually pertains to the inner loop's counter. Jun 3 '20 at 23:17
• @EmmaX After thinking about it for a little while, yeah it's not necessary to have the counters outside of the for loops... If there were to be any kind of gain it would be minimal at best and the tradeoff for readability simply outweighs it. So I reverted my code to put them back in their respective for loops and I removed the highlighted bullet that mentions them. Jun 4 '20 at 1:20