Binomial Theorem implementation, including complex numbers

Question

#include <iostream>
#include <complex>
#include <vector>
#include <cmath>
/*
 adapted from https://stackoverflow.com/questions/55421835/c-binomial-coefficient-is-too-slow
*/
template<typename T = int, typename U = unsigned long long int>
inline U computeBinomialCoefficient(const T& n, const T& k) {
 if (n == k || k == 0) {
 return 1;
 }
 
 U out = n - k + 1;
 for (T i = 1; i < k; ++i) {
 out = out * (n - k + 1 + i) / (i + 1);
 }
 return out;
}
/*
 compute the horizontal shift of a polynomial
 with coefficients coeffs (decreasing order)
 by shift units
*/
template<typename T = int>
const std::vector<T> computeHorizontalShiftBinomialThm(const std::vector<T>& coeffs, const T& shift) {
 std::size_t size = coeffs.size();
 std::vector<T> out (size, 0);
 
 for (std::size_t i = 0; i < size; i++) {
 T coefficient = coeffs[i];
 if (coefficient == 0) {
 continue;
 }
 
 std::size_t n = size - i - 1; // power of monomial
 for (std::size_t k = 0; k <= n; k++) {
 T kthPower = computeBinomialCoefficient(n, k) * powl(-1 * shift, n - k);
 out[size - k - 1] += coefficient * kthPower;
 }
 }
 
 return out;
}
/*
 Specialization for std::complex
*/
template<typename ComplexType>
const std::vector<std::complex<ComplexType>> computeHorizontalShiftBinomialThm(const std::vector<std::complex<ComplexType>>& coeffs, const std::complex<ComplexType>& shift) {
 using imag = std::complex<ComplexType>;
 
 std::size_t size = coeffs.size();
 std::vector<imag> out (size, 0);
 
 for (std::size_t i = 0; i < size; i++) {
 imag coefficient = coeffs[i];
 if (coefficient == std::complex<ComplexType>(0, 0)) {
 continue;
 }
 
 std::size_t n = size - i - 1; // power of monomial
 for (std::size_t k = 0; k <= n; k++) {
 ComplexType binomialCoeff = computeBinomialCoefficient(n, k);
 imag power = (imag) pow(shift * (ComplexType) -1, n - k);
 imag kthPower = (std::complex<decltype(binomialCoeff)>) power * binomialCoeff;
 out[size - k - 1] += coefficient * kthPower;
 }
 }
 
 return out;
}
template<typename T>
void print(const std::vector<T>& vec) {
 for (const T& n : vec) {
 std::cout << n << " ";
 }
 std::cout << "\n";
}
int main()
{
 std::vector<int> intCoeffs {2, 19, 0, 0, 1};
 int intShift {5};
 auto intExample = computeHorizontalShiftBinomialThm(intCoeffs, intShift);
 print(intExample);
 
 
 std::vector<std::complex<double>> imaginaryCoeffs {1, 0, 0, 0};
 std::complex<double> imaginaryShift (0, 2);
 
 auto complexExample = computeHorizontalShiftBinomialThm(imaginaryCoeffs, imaginaryShift);
 print(complexExample);
 
 return 0;
}

And it works! I get out

2 -21 15 425 -1124 
(1,0) (3.67394e-16,-6) (-12,-1.46958e-15) (-1.46958e-15,8)

How can I reduce/eliminate floating point error, increase speed, and make the code easier to read?

Answer 1

Avoid needing to specialize

It is possible to make just one template that can compute the desired result without having a specialization for std::complex. You just have to ensure that variables get their correct type. Some tricks you can use:

• Use T{} to get a default initialized value of type T; this is the way to get a correctly typed zero.
• Use auto where possible to let the compiler deduce the right type.
• Use operators or functions to get a result of the right type, for example:
  • Write -foo instead of -1 * foo.
  • Use decltype(std::abs(value)) to get a scalar type even if value is std::complex.

Here is your template rewritten to work with std::complex:

template<typename T = int>
const std::vector<T> computeHorizontalShiftBinomialThm(const std::vector<T>& coeffs, const T& shift) {
 std::size_t size = coeffs.size();
 std::vector<T> out (size, 0);
 
 for (std::size_t i = 0; i < size; i++) {
 T coefficient = coeffs[i];
 if (coefficient == T{}) {
 continue;
 }
 
 std::size_t n = size - i - 1;
 for (std::size_t k = 0; k <= n; k++) {
 T binomialCoeff = computeBinomialCoefficient(n, k);
 auto power = std::pow(-shift, n - k);
 auto kthPower = power * binomialCoeff;
 out[size - k - 1] += coefficient * kthPower;
 }
 }
 
 return out;
}

The only issue with the above code is that you want extended precision for some temporary variables like binomialCoeff and power. In particular you want to use long long if T is an integer type, long double if it's a floating point type, and std::complex<long double> if it's a complex number. To tackle this, I would create a separate template function to compute kthPower, and then this could be specialized for integers, floating point and complex types.

Use std:: versions of math functions

Use std::pow() instead of pow(). The C++ versions of the math functions have overloads for the various floating point types, so you don't have to worry about whether you are passing in a float, double, long double. Even better, std::pow() has overloads where the first argument is a floating point number but the second is an integer, so it can choose a more optimal implementation for raising values to an integer power.

The same goes for all other function from <cmath>.

Avoid unnecessary copies

You already pass the parameters by reference, which is great, especially if you pass containers or types of which you don't know the size. However, there are a few other places where copies could be made that might be expensive depending on the type.

The first is in the for-loop in computeBinomialCoefficient(). Use out *= ... instead of out = out * ....

The other is coefficient in computeHorizontalShiftBinomialThm(). Use this instead:

const T& coefficient = coeffs[i];

For simple types like float, double and even std::complex types, the compiler will likely produce exactly the same code either way, but consider for example someone using a BigInt library and wanting to process very large numbers.

Naming things

Your functions have very long names. Consider finding ways to make them more concise while still being clear and descriptive given the context they would be used in. I would drop the prefix compute; you won't see a compute_sin() function in the standard library for example. Use binomial() instead of computeBinomialCoefficient(). I think it is usually clear that a function with that name that takes two parameters produces the binomial coefficient as the output. As for the function that does the horizontal shift, it might be a bit hard to come up with a better name apart from dropping the compute part.

Be consistent when naming things. I see both coefficient and the abbreviation coeff being used. Choose one and use it everywhere.

Answer 2

One suggestion that I think is good is to cache the Ncr results.

struct Ncr {
 using T = long double;
 static std::vector<T> &fact() {
 static std::vector<T> _fact{1};
 return _fact;
 }
 static std::vector<T> &inv_fact() {
 static std::vector<T> _inv_fact{1};
 return _inv_fact;
 }
 static T C(int n, int k) {
 if (k < 0 || k > n) {
 return 0;
 }
 while ((int) fact().size() < n + 1) {
 fact().push_back(fact().back() * (int) fact().size());
 inv_fact().push_back(1 / fact().back());
 }
 return fact()[n] * inv_fact()[k] * inv_fact()[n - k];
 }
};

I have added an example code and to use it you can directly call it by Ncr::C(n, k) this will save all the previous results and you might not have to recompute every time. You have to make sure you don't actually run this with int or any other small-sized DS as factorials encreases polynomially and can overflow in case of int or long long. That's why using long double might be helpful here.
Also, notice the long double this is actually to improve the precision. When I run the code I got the following output.

2 -21 15 425 -1124 
(1.000000,0.000000) (0.000000,-6.000000) (-12.000000,-0.000000) (-0.000000,8.000000)