Improve polynomials

src/calculus/derivation.h

@@ -20,10 +20,16 @@ namespace theoretica {
 	template<typename T = real>
 	inline polynomial<T> deriv_polynomial(const polynomial<T>& p) {
 
+		if(p.coeff.size() == 0) {
+			TH_MATH_ERROR("deriv_polynomial", p.coeff.size(), INVALID_ARGUMENT);
+			return polynomial<T>({T(nan())});
+		}
+
 		polynomial<T> Dp;
+		Dp.coeff.resize(p.coeff.size() - 1);
 
 		for (unsigned int i = 1; i < p.size(); ++i)
-			Dp.coeff.push_back(p.get(i) * i);
+			Dp.coeff[i - 1] = p[i] * i;
 
 		return Dp;
 	}

src/calculus/integration.h

@@ -24,10 +24,11 @@ namespace theoretica {
 	inline polynomial<T> integrate_polynomial(const polynomial<T>& p) {
 
 		polynomial<T> P;
-		P.coeff.push_back(0);
+		P.coeff.resize(p.size() + 1);
+		P[0] = 0;
 
 		for (unsigned int i = 0; i < p.size(); ++i)
-			P.coeff.push_back(p.get(i) / T(i + 1));
+			P[i + 1] = p[i] / T(i + 1);
 
 		return P;
 	}

src/polynomial/polynomial.h

@@ -54,6 +54,12 @@ namespace theoretica {
 			}
 
 
+			/// Return the nth order coefficient
+			inline const T& operator[](unsigned int i) const {
+				return coeff[i];
+			}
+
+
 			/// Return the nth order coefficient
 			inline T& operator[](unsigned int i) {
 				return coeff[i];

test/test_autodiff.cpp

@@ -137,7 +137,7 @@ int main(int argc, char const *argv[]) {
 			[=](interval k, Real tol, unsigned int n) {
 				return test_operator(
 					[=](vec<> v) {
-						return gradient(f, v) * J * gradient(H, v);
+						return gradient(f, v) * (J * gradient(H, v));
 					},
 					PRNG::xoshiro(time(nullptr)),
 					k, 1E-8, n, 4