Wednesday, 27 November 2013

Some STL functions are Simpler Than You Might Think

The STL functions generate and generate_n are quite useful, but when I first used them I thought I must be missing something. It turns out though that these functions, like many in the standard library, are really simple.

The first time I used generate_n was to initialise a vector. I looked at the documentation and quickly found how to achieve this, but sub-consciously I assumed I was only scratching the surface of what it could do. I probably thought that way about all STL functions, and in many cases I would have been right. For an example try stepping into the line:

std::cout << "Hello, world\n";

If you understand everything that happens in the above line of code you're probably someone who gets paid to implement the standard library. The fact that the operator<< is from the std namespace should have you scratching your head until you've read about Name Lookup, more specifically Koenig Lookup.

However, with generate_n, things are much more straight forward. In fact, generate_n is most easily understood, and in my case demystified, by looking at the implementation. Here's the crux of the function from VC++ 2012:

template<class _OutIt,
  class _Diff,
  class _Fn0> inline
  _OutIt _Generate_n(_OutIt _Dest, _Diff _Count, _Fn0 _Func)
  {
  for (; 0 < _Count; --_Count, ++_Dest)
    *_Dest = _Func();
  return (_Dest);
  }


That's as difficult as it gets

The above code is probably as easy to understand as any documentation on generate_n. But, more importantly, when you get this, you understand the function completely. There's nothing else to it, there's no smoke nor mirrors, no __asm__, and it doesn't link to any binaries to do the real work. It's just a small bit of C++ code.

The generate function is similar except it expects two input iterators and the for loop looks like:

for (; _First != _Last; ++_First)

If you look at the standard itself this shouldn't come as a surprise. After the function declaration there are just seven lines detailing what both generate and generate_n need to conform to. Many other useful STL functions have a similar level of complexity.

This won't come as a surprise to everyone, but for me it was encouraging to know that although I may never understand the language completely, there are many parts that I can get to the bottom of quite easily. It's also encouraged me to step into standard library functions as a first port of call, and shown me a less intimidating side of the C++ Standard itself. That said, I think I need some more experience before I'm ready to tackle:

std::cout << "Hello, world\n";

Tuesday, 12 November 2013

Understanding Debug Iterators in VC++ 2012

In debug mode Visual C++ STL iterators can provide useful information at run-time via failed assertions, but there’s a performance cost and in certain scenarios they can seem to grind your program to a halt. This post shows why that happens, and how to fix it. Consider the following code:

vector<int> myInts;
myInts.push_bach(1);
cout << myInts[4] << endl;

It will compile O.K., and when done so in release mode the resultant program will happily reach into memory beyond the myInts vector and send whatever garbage it finds there to cout. However, in debug mode stepping over the third line produces the popup “Debug Assertion Failed… vector subscript out of range”. That may not be a big deal in the above case but if the index were calculated in a loop dependent on a number of variables that information could be very useful.

For another example consider what’s wrong with now adding the following lines:

auto it = myInts.begin();
myInts.push_back(2);
cout << *it << endl;

Chances are this code doesn't do what you want it to because pushing the 2 onto the vector probably requires it to grow, which requires it to move, which leaves the iterator it pointing at unallocated memory. Again, in release mode cout will get whatever it is now pointing at, but in debug mode we get “vector iterator not dereferencable”. While dereferencable may not be a real word you’ll soon figure out what it means and it will probably save some debugging time.

However, these extra checks come at a cost. Pretty much any call to an iterator involves a lot of run-time checks. For example dereferencing an iterator involves checking that the parent container isn’t null, checking the iterator’s internal pointer isn’t null, and checking that the iterator doesn’t point to memory before the parent vector’s first iterator, or after its last iterator. The parent container also has extra concerns when debug iterators are being used. All iterators that point into one container form a singly linked list. This allows the parent container to communicate with those iterators, but walking the list when it needs to is slow.

Example of slow iterators in debug mode

As a simple example, consider a list of 10,000 numbers from 1-100. I wish to find the number of pairs within that list that sum to 100. Here’s some code to achieve that:

auto noCombosSumTo100 = 0;
for (auto it1 = begin(myNums); it1 != end(myNums); it1++) {
  for (auto it2 = it1; it2 != end(myNums); it2++) {
    if (*it1 + *it2 == 100) {
      noCombosSumTo100++;
    }
  }
}

On my machine this takes about 34 milli-seconds fully optimised in release mode, and about 75 seconds in debug. That's about 2200x difference between debug and release!

Obviously this difference can't be entirely attributed to the iterators, there are lots of other opportunities for the compiler to optimise the above code in release mode. However, if I implement the same algorithm using vector indexing, like this:

for (auto index1 = 0; index1 < numTerms; index1++) {
  for (auto index2 = index1; index2 < numTerms; index2++) {
    if (mynums[index1] + mynums[index2] == 100) {
      noCombosSumTo100++;
    }
  }
}

it runs in about 4.5 seconds, so there's clearly a big cost to using debug iterators.

Ways to improve iterator performance

The improvements I'm about to show for the first code snippet could apply to any code, but it's worth understanding that any inefficiencies are exacerbated when debug iterators are involved. To improve the first code snippet while still getting the benefit of debug iterators I can:

These three modifications bring the execution time down to 21 seconds. Interestingly none of them make any difference in the second snippet where built in types are involved. Again, it's worth understanding that these improvements would apply to any user defined type (in debug mode), but it probably wouldn't be just as noticeable as when iterators are involved.

If you want a further speedup then you'll need to sacrifice some additional run-time checks by setting _ITERATOR_DEBUG_LEVEL=1 or 0 (This macro supersedes and combines the functionality of the older _SECURE_SCL and _HAS_ITERATOR_DEBUGGING macros). Doing so decreases the amount of debug help you get, but it brings the original example's execution time down to about 8.5 seconds.

At 8.5 seconds it's still about twice as slow as indexing the vectors, and this is because even with the reduced checking there are more indirections (extra function calls) associated with dereferencing a debug iterator.

In summary...

Remember, the example code used here takes 34 milli-seconds in release mode. The example using vector indexing also takes about 34 milli-seconds in release, and the same algorithm indexing a raw array of ints takes the same time to execute, about 34 milli-seconds. This indicates that the code is being optimised to similar assembly in all three cases.

So, to summarise, STL iterators provide useful debug information for development, they allow you to access STL containers as if you are using a pointer, and they create production code that is as fast as using raw pointers. However, there is a cost in debug mode as they can slow your code down a lot. But, if these costs become too much of a hindrance there are ways to reduce the effect without changing how you are using the iterators.

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.

Sunday, 6 October 2013

Computing the Collatz sequence

Problem 14 from ProjectEuler concerns the Collatz sequence. From the problem's description:

The following iterative sequence is defined for the set of positive integers:
n → n/2 (n is even)
n → 3n + 1 (n is odd)

This shows how to compute the Collatz sequence for n. As an example, the Collatz sequence for 13 is:
13 → 40 → 20 → 10 → 5 → 16 → 8 → 4 → 2 → 1

The problem is to find which number, under one million, has the longest Collatz sequence. I took the following approach:

From the above example, the length of the Collatz for 13 is equal to the length for 40's  sequence plus one. In general, the length of the Colaltz sequence for n is equal to the length for second term of that sequence plus one. This can be implemented neatly with the recursive function:
uint getCollatzSeriesLength(uint value)
{
  uint seriesLength;

  if (value == 1) {
    seriesLength = 1;
  }
  else {
    if (value%2 == 0) {
      seriesLength = getCollatzSeriesLength(value/2) + 1;
    }
    else {
      seriesLength = getCollatzSeriesLength(value*3+1) +1;
    }
  }
  
  return seriesLength;
}

The recursion is terminated by returning a length of 1 for a value 1. Using the above function for it takes about 580 milliseconds to calculate the Collatz series length for every integer from 1 to 1,000,000 on my machine.

I then made the following optimisation. Consider again the Collatz series of 13:
13 → 40 → 20 → 10 → 5 → 16 → 8 → 4 → 2 → 1
If the program starts at 1 then when it comes to 13 it has already calculated the series length for 10. If this value was cached then it could be returned without recursing all the way to the end of the series. Therefore when when the algorithm calculates a series length it should cache that, and the function should check that cache before calculating the series length for the next term. Note an upshot of this is that the series length's for 40 and 20 will be calculated in the course of calculating the series length of 13. The updated function looks like this (changes highlighted):
uint getCollatzSeriesLength(uint value)
{
  uint seriesLength;

  if ((value <= noTerms) && (lengthsCache[value] != 0)) {
    seriesLength = lengthsCache[value];
  }
  else if (value == 1) {
    seriesLength = 1;
  }
  else {
    if (value%2 == 0) {
      seriesLength = getCollatzSeriesLength(value/2) + 1;
    }
    else {
      seriesLength = getCollatzSeriesLength(value*3+1) +1;
    }
    if (value <= noTerms) {
      lengthsCache[value] = seriesLength;
    }
  }

  return seriesLength;
}

Here noTerms is the number of terms to be checked (one million) and lengthsCached is a vector initialised with one million and one zeros (so it's effectively indexed from one). This optimisation makes around a 21x difference in terms of speed, as it now takes about 28 milliseconds to run the program.

I made a final optimisation by considering that when the series length has been calculated for some value, n, then the series length for 2n will be one greater. Using this the following ammendment was made to the code:
uint getCollatzSeriesLength(uint value)
{
  uint seriesLength;

  if ((value <= noTerms) && (lengthsCache[value] != 0)) {
    seriesLength = lengthsCache[value];
  }
  else if (value == 1) {
    seriesLength = 1;
  }
  else {
    if (value%2 == 0) {
      seriesLength = getCollatzSeriesLength(value/2) + 1;
    }
    else {
      seriesLength = getCollatzSeriesLength(value*3+1) +1;
    }
    auto additionalSteps = 0;
    for (auto i=value; i<=noTerms; i*=2) {
      lengthsCache[value] = seriesLength + additionalSteps++;
    }
  }

  return seriesLength;
}

This brings the program's execution time down by a further factor of ~2.3x to around 12 milliseconds on my machine.

Wednesday, 2 October 2013

Project Euler Progress

I've been working through the problems on Project Euler and as I was making notes anyway I thought I'd put them in a log, on the web, hence this blog. If it serves no other purpose it will be somewhere to jog my memory on the things I've learned.

Project Euler (in case you haven't clicked on the link above) is a good way to exercise your programming chops. The problems are maths based but generally can't be solved with pen and paper as they involve too much work. For example, the first problem is "Find the sum of all multiples of 3 or 5 below 1000", and a simple C++ solutions is:
 int total = 0;
 for (int i=0; i<1000; i++) {
  if ((i%3==0) || (i%5==0)) {
   total += i;
  }
 }
 cout << "Answer is " << total << endl;
They get progressively more difficult from there but the emphasis stays on algorithms as opposed to mathematics.

To get started, here's my solution to Problem 24.


What is the millionth lexicographic permutation of the digits 0, 1, 2, 3, 4, 5, 6, 7, 8 and 9?
I began by thinking about the simple case of the permutations of the digits 0, 1, and 2, let's say we want to find the fourth permutation. There are 3! (i.e. 6) elements in the list, and listed in lexicographic order they are:

012, 021, 102, 120, 201, 210

We can split the list into three sub-sections: two elements beginning with 0, two beginning with 1, and two beginning with 2. It's easy to see that the fourth element in the list is in the second of these sub-sections. Stated in terms of a zero indexed array: position 3 in the array is in position 1 of the array of sub-sections. We can use the size of these sub-sections, 2!, to calculate this result as 3/2! = 1 (because integers are truncated after division this works when the numbers don't divide evenly). The elements in this sub-sections start with the digit 1 and so the first digit of the answer is a 1.

Now we descend into the sub-section 102, 120. We just need to remove the digit 1 from the list of digits we're considering, recalculate the number of elements in this sub-section, and recalculated the position of our answer in this sub-section, then we can repeat the process until we've got the answer. The array is one element shorter than before so the new size is 2!. We can essentially discount all sub-sections before the one we're now in so the new position is the previous position minus 2!*1. Repeating the process brings us to the fourth element, 120.

This 3 digit example is trivial but the algorithm can be applied to a list with an arbitrary number of digits. Here's the C++ code to find the 1,000,000th lexicographic permutation of the digits 0 to 9:
auto pos = 999999; // 1,000,000th position, zero indexed
uint numDigits = 10;
vector<uint> digits(numDigits);
emuint n = 0;
generate(digits.begin(), digits.end(), [&n](){return n++;});

cout << "Answer is ";
for (uint i=numDigits; i>1; i--) {
  uint nextDigit = pos/Factorial(i-1);
  cout << digits[nextDigit];
  digits.erase(digits.begin()+nextDigit);
  pos -= nextDigit*Factorial(i-1);
}
cout << digits[0] << endl;
Note, this uses a factorial function, like one you could find here.