Showing posts with label Maths Utility Functions. Show all posts
Showing posts with label Maths Utility Functions. Show all posts

Friday, 1 November 2013

Fast Implementation of the Factorial function in C++

I’ve been investigating the fastest way to implement the factorial function in my MathsUtils library. I had originally implemented it as simply:

template <typename T>
emu64 MathUtils::RecursiveFactorial(const T& x)
{
  return (x <= 1 ? 1 : x * Factorial(x - 1));
}

although I suspected this wasn’t the most efficient way, and I wanted to measure it.
To measure the efficiency of this implementation I used the following benchmark test:

for (ulonglong i=1; i<10000000; i++) {
  for (ulonglong j=1; j<20; j++) {
    RecursiveFactorial(j);
  }
}

(I don't find the factorial for anything greater than 20 as it would be too large to represent with a 64-bit unsigned number)

First I took a closer look at what the compiler was doing with the above factorial implementation, and the assembly for a debug build was as expected, it implements recursion. With this implementation the above test takes about 14 seconds on my machine in debug mode.

I found that the assembly the compiler writes for this function doesn't change regardless of if a const value is passed in. This would change in the future if I mark the function with constexpr. VC++ 2012 doesn’t support this feature yet but it would allow the compiler to evaluate a function's return value so you get it for free at run time.

I then thought about inlining the function. If the compiler did this blindly then it might hang while writing an inline function, inside an inline function, etc… Compilers are smarter than this though and looking at the assembly showed that using inline didn't make any difference to the output produced by the compiler. This is understandable as the inline keyword is only a hint to the compiler, one it is free to ignore.

I found, however, that it is possible to make the compiler inline a recursive function using the inline_recursing pragma (for VC++, there are equivalents for other compilers). When this is used the recursive function will be unrolled and written inline up to a maximum of 16 times, and this is also the default inline recursion depth. The depth can be changed using inline_depth.This Stack Overflow answer gives a good explanation of how a recursive function can be inlined to a depth of three.

To compile using these pragmas you need to turn off some debugging options, so after I'd seen from the assembly that the function was indeed being inlined I switched to release mode. Here are the results for the original RecursiveFactorial() function on the benchmark test:
  • Without inlining - 7.9 seconds
  • With inlining and a depth of five - 6.7 seconds.
I found that an inline recursion depth of five gave the best results.

As an interesting aside, the above test, shown here again:

for (ulonglong i=1; i<10000000; i++) {
  for (ulonglong j=1; j<20; j++) {
    RecursiveFactorial(j);
  }
}

actually runs instantaneously in release mode. This is because, in release mode, the compiler detects that the call to RecursiveFactorial() is dead code, so it's eliminated from the compiled assembly. To force it to execute I had to add the result to a running total, and even then the calls were skipped unless I used the total afterwards, so I printed it to the console.

Having spent some time tweaking the recursive factorial function I still had a suspicion that a simple for loop would be faster, so I tested the following:

template <typename T>
ulonglong ForLoopFactorial(const T& value)
{
  ulonglong returnValue = 1;
  for (ulonglong i = 2; i <= value; i++) {
    returnValue *= i;
  }
  return returnValue;
}

In release mode this takes 5.1 seconds to run on my machine. So, despite digging into the first option, and learning a few things, the simpler for loop wins for speed, so that's what I've updated my MathUtils library to use. Unit tests still pass so I haven't broken anything.

I've put the code for this spike on my GitHub.

Monday, 21 October 2013

Some Math Utility Functions (on GitHub) to help out with Project Euler

I've been working through the problems on the Project Euler website recently, and as there were some maths functions that I needed over and over I decided to put them into a (C++) utility namespace, MathUtils, which I've put up on GitHub. This post is a short overview of what MathsUtils currently provides. The code should (hopefully) be self explanatory so I'll not provide much more detail here.

(I've also got a class up there called embiguint, which is an unsigned integer that can be arbitrarily large. It has proved useful in its current state but I intend to add a lot more functionality and make it behave more like a built in type at which point I'll do a blog post on it too.)

The functions currently included in MathsUtils are:

template <typename T> std::string
                     ToAlternateBaseRepresentation(
                     const T& value, const T& base);
template <typename T> std::vector<T>
                     GetDivisors(const T& number);
template <typename T> emu64 GetSumOfDivisors(const T& number);
template <typename T> emu64 Factorial(const T& x);
template <typename T> T
                     PowerOf(const T& base, const T& power);
template <typename T> std::vector<bool>
                     GetPrimesBoolArrayToN(const T& n);
template <typename T> std::vector<T> GetPrimesToN(const T& n);

I had originally developed the functions in a class with static members but I changed this to a namespace as it seemed more more natural. After all, classes are intended to define how to build objects whereas namespaces are used to group related functionality. In Effective C++, in an item entitled Perfer non-member non-friend functions to member functions, Scott Meyers argues that this approach helps to increase extensiblity and (counter-intuitively) encapsulation.

All the functions are templated as I want to make them as easy to use as possible. There's still lots more I could add to these functions, for example the PowerOf functions only works with non-negative exponents, however each function was written or extended with a particular problem in mind. I'll add more functionality as I need it and in the mean time the functions should throw a useful exception if you ask them to do something they can't do.

There's also lots more I can do with regard to optimisation, and this should give me an opportunity to try out (and learn) some new C++11 features. For example PowerOf and Factorial could benefit from constexpr, although I've just been looking at the release of Visual Studio 2013 and I don't think this feature has made it in.  I guess I can achieve the same result with template metaprogramming if I really want quicker run-time code.

Anyway, my MathUtils are in place and poised to help out on more Project Euler problems.