Add point-wise derivative and integral of polynomials

src/calculus/deriv.h

@@ -38,6 +38,27 @@ namespace theoretica {
 	}
 
 
+	/// Compute the exact derivative of a polynomial function at the given point.
+	/// In-place calculation with Horner's evaluation scheme is used,
+	/// with linear complexity in the coefficients \f$O(n)\f$.
+	///
+	/// @param p The polynomial to differentiate
+	/// @param x The point to compute the derivative at
+	/// @return The derivative of the polynomial at the given point
+	inline real deriv(const polynomial<real>& p, real x) {
+
+		real dp = 0.0;
+
+		for (unsigned int i = 0; i + 1 < p.size(); ++i) {
+
+			const unsigned int pos = p.size() - i - 1;
+			dp = pos * p[pos] + x * dp;
+		}
+
+		return dp;
+	}
+
+
 	/// Approximate the first derivative of a real function
 	/// using the central method.
 	///

src/calculus/integral.h

@@ -16,26 +16,50 @@
 namespace theoretica {
 
 
-	/// Compute the indefinite integral of a polynomial
+	/// Compute the indefinite integral of a polynomial.
 	///
 	/// @param p The polynomial to integrate
-	/// @return The indefinite polynomial integral
-	template<typename T = real>
-	inline polynomial<T> integral(const polynomial<T>& p) {
+	/// @return The indefinite polynomial integral, with zero constant term.
+	template<typename Field = real>
+	inline polynomial<Field> integral(const polynomial<Field>& p) {
 
-		polynomial<T> P;
+		polynomial<Field> P;
 		P.coeff.resize(p.size() + 1);
-		P[0] = 0;
+		P[0] = Field(0.0);
 
 		for (unsigned int i = 0; i < p.size(); ++i)
-			P[i + 1] = p[i] / T(i + 1);
+			P[i + 1] = p[i] / Field(i + 1);
 
 		return P;
 	}
 
 
+	/// Compute the definite integral of a polynomial over an interval.
+	/// In-place calculation with Horner's evaluation scheme is used,
+	/// with linear complexity in the coefficients \f$O(n)\f$.
+	///
+	/// @param p The polynomial to integrate
+	/// @param a The lower extreme of integration
+	/// @param b The upper extreme of integration
+	/// @return The definite polynomial integral
+	inline real integral(const polynomial<real>& p, real a, real b) {
+
+		real P_a = 0.0;
+		real P_b = 0.0;
+
+		for (unsigned int i = 0; i < p.size(); ++i) {
+
+			const unsigned int pos = p.size() - i - 1;
+			P_a = a * (p[pos] / (pos + 1) + P_a);
+			P_b = b * (p[pos] / (pos + 1) + P_b);
+		}
+
+		return P_b - P_a;
+	}
+
+
 	/// Approximate the definite integral of an arbitrary function
-	/// using the midpoint method
+	/// using the midpoint method.
 	///
 	/// @param f The function to integrate
 	/// @param a The lower extreme of integration
@@ -62,7 +86,7 @@ namespace theoretica {
 
 
 	/// Approximate the definite integral of an arbitrary function
-	/// using the trapezoid method
+	/// using the trapezoid method.
 	///
 	/// @param f The function to integrate
 	/// @param a The lower extreme of integration
@@ -127,7 +151,7 @@ namespace theoretica {
 
 
 	/// Approximate the definite integral of an arbitrary function
-	/// using Romberg's method accurate to the given order
+	/// using Romberg's method accurate to the given order.
 	///
 	/// @param f The function to integrate
 	/// @param a The lower extreme of integration

src/optimization/roots.h

@@ -185,11 +185,12 @@ namespace theoretica {
 	/// the algorithm (defaults to OPTIMIZATION_NEWTON_ITER).
 	/// @return The coordinate of the root of the function,
 	/// or a complex NaN if the algorithm did not converge.
-	template<typename Type = real>
+	template <
+		typename Type = real,
+		typename ComplexFunction = std::function<complex<Type>(complex<Type>)>
+	>
 	inline complex<Type> root_newton(
-		complex<Type>(*f)(complex<Type>),
-		complex<Type>(*df)(complex<Type>),
-		complex<Type> guess = complex<Type>(0.0, 0.0),
+		ComplexFunction f, ComplexFunction Df, complex<Type> guess,
 		real tol = OPTIMIZATION_TOL,
 		unsigned int max_iter = OPTIMIZATION_NEWTON_ITER) {
 
@@ -261,7 +262,8 @@ namespace theoretica {
 	/// @param f The real function to search the root of,
 	/// with dual2 argument and return value.
 	/// @param guess The initial guess (defaults to 0).
-	/// @param tol The \f$\epsilon\f$ tolerance value to stop the algorithm when \f$|f(x_n)| < \epsilon\f$
+	/// @param tol The \f$\epsilon\f$ tolerance value to stop the algorithm
+	/// when \f$|f(x_n)| < \epsilon\f$.
 	/// @param max_iter The maximum number of iterations to perform, after which
 	/// the algorithm is assumed to not have converged.
 	/// @return The coordinate of the root of the function,

test/prec/test_calculus.cpp

@@ -154,6 +154,22 @@ int main(int argc, char const *argv[]) {
 		);
 	}
 
+
+	{
+		polynomial<real> P = {};
+		prec::equals("deriv(P{}, x)", deriv(P, 1.0), 0.0);
+
+		P = {1.0};
+		prec::equals("deriv(P{const.}, x)", deriv(P, 1.0), 0.0);
+
+		P = {1.0, 1.0, 1.0, 1.0};
+		prec::equals("deriv(P, x)", deriv(P, 1.0), 6.0);
+
+		P = {+1.4, -3.2, 12.7, 7.5, 11,7};
+		prec::equals("deriv(P, x)", deriv(P, 1.0), deriv(P)(1.0));
+	}
+
+
 	{
 		auto deriv_opt = prec::estimate_options<real, real>(
 			prec::interval(0.001, 0.5),
@@ -192,6 +208,52 @@ int main(int argc, char const *argv[]) {
 		// Compare integral quadrature to primitives
 
 
+	{
+		auto opt = prec::equation_options<polynomial<real>>(
+			1E-08, distance_polyn
+		);
+
+		polynomial<real> p = {1, 1};
+		polynomial<real> expected = {0, 1, 0.5};
+
+		prec::equals(
+			"integral(P)",
+			integral(p),
+			expected,
+			opt
+		);
+
+		prec::equals(
+			"integral(0)",
+			integral(polynomial<real>({0.0})),
+			polynomial<real>({0.0}),
+			opt
+		);
+
+		prec::equals(
+			"integral(P{})",
+			integral(polynomial<real>({})) == polynomial<real>({}),
+			true
+		);
+	}
+
+
+	{
+		polynomial<real> P = {};
+		prec::equals("integral(P{}, x)", integral(P, 0.0, 1.0), 0.0);
+
+		P = {1.0};
+		prec::equals("integral(P{const.}, x)", integral(P, 0.0, 1.0), 1.0);
+
+		P = {1.0, 2.0, 3.0, 4.0};
+		prec::equals("integral(P, x)", integral(P, 0.0, 1.0), 4.0);
+
+		P = {1.2, 4.5, -8.8, 11.3, -9.6};
+		auto I_P = integral(P);
+		prec::equals("integral(P, x)", integral(P, -7.0, 5.0), I_P(5.0) - I_P(-7.0));
+	}
+
+
 	{
 		auto integ_opt = prec::estimate_options<real, real>(
 			prec::interval(0.1, 3.0),
@@ -260,35 +322,6 @@ int main(int argc, char const *argv[]) {
 		);
 	}
 
-	{
-		auto opt = prec::equation_options<polynomial<real>>(
-			1E-08, distance_polyn
-		);
-
-		polynomial<real> p = {1, 1};
-		polynomial<real> expected = {0, 1, 0.5};
-
-		prec::equals(
-			"integral(P)",
-			integral(p),
-			expected,
-			opt
-		);
-
-		prec::equals(
-			"integral(0)",
-			integral(polynomial<real>({0.0})),
-			polynomial<real>({0.0}),
-			opt
-		);
-
-		prec::equals(
-			"integral(P{})",
-			integral(polynomial<real>({})) == polynomial<real>({}),
-			true
-		);
-	}
-
 
 		// ode.h
 		// Integrate the simple harmonic oscillator