Haskell Prime Number Generator (part 1 of 3)

Project Euler Problem 7 reads "What is the 10 001st prime number?". Solving this problem requires a prime number generator. I solved this using Haskell, which I'm currently learning, but after I'd done so I learned a lot more by looking at others' solutions. Here I present some of the Haskell prime number generators being discussed in the Project Euler forum, problem 7 thread (I can't link to the forum as you need to have an account and have solved Problem 7 to access that thread):

As I'm just learning Haskell my aim is to annotate these examples at my own level, i.e. for beginners. Note, the examples I'm referring to were meant to solve a particular problem, so in some cases I've tweaked them slightly to make them general prime number generators.

Prime Number Generator from user gauchopuro:

primes (x:xs) = x : primes (filter ((/= 0) . (`mod` x)) xs)

This function must be supplied with a list of the positive integers starting at 2 and it uses the Sieve of Eratosthenes to calculate the list of the primes. For example a function that uses the above to get the first n prime numbers would look like this:

getFirstNPrimes n = take n $ primes [2..]

Or to get the nth prime:

getNthPrime n = primes [2..] !! (n-1)

Although very concise, once you understand the notation involved, primes is also very expressive. Here's what you need to know to parse it:

1. primes uses pattern matching. It expects to be supplied with a list of integers and the(x:xs)to the left of the equals sign means that within the function (i.e. to the right of the equals sign) x will refer to the first element in the list and xs (pronounced exes, as in the plural of x) refers to the remainder of the list.
2. primes is recursive. The list that it generates consists of the first element of the list passed into it and the output of primes when passed the remainder of the list after some filtering.
3. The filter function expects two arguments. The first argument is a predicate and the second is a list. filter returns a new list consisting of the elements from the initial list that yield True when passed to the predicate.
4. The . operator composes two functions i.e. (f . g) x == f (g x).
5. The cons operator ':' appends a new element to the start of a list.

Haskell Lazy Evaluation

Although not strictly required by the primes function shown, it's also useful to understand Haskell's lazy evaluation. It means that Haskell will do just enough work to get to the result it needs. For example, consider the following function:

fib :: Int -> Int -> [Int]
fib a b = a : fib b (a+b)

The first line tells us that this function takes two integers and returns a list of integers. The second line defines the function and says that the list to return is made up of the first integer argument followed by the list produced by calling fib again recursively, with the original second argument as the new first argument, and the sum of the original arguments as the new second argument. This function can be called with 1 1 to generate the Fibonacci sequence as follows:

fib 1 1 = 1 : fib 1 2
        = 1 : 1 : fib 2 3
        = 1 : 1 : 2 : fib 3 5
        = 1 : 1 : 2 : 3 : fib 5 8 etc...

You can see how this generates the Fibonacci sequence, but if you're not familiar with lazy evaluation you'd expect this function to call itself recursively indefinitely and never return. However, lazy evaluation ensures that the function does just enough work to get the result it needs. So if you just want the first 10 Fibonacci numbers you call the function with:

take 10 (fib 1 1)

and the result is:

[1,1,2,3,5,8,13,21,34,55]

After the first 10 elements are generated Haskell has the result it needs and so fib returns without continuing to call recursively.

Understanding primes

Here is primes again:

primes (x:xs) = x : primes (filter ((/= 0) . (`mod` x)) xs)

and, again, it is to be called with the list of integers starting at two, for example like:

primes [2..] !! 101          -- To get the 100th prime

With the above knowledge you can now understand primes as implementing the Sieve of Eratosthenes as follows:

1. It's originally passed a list of integers starting at 2 (2, 3, 4, 5...). The first element of that list, x, will originally be 2. This becomes the first element in the list of prime numbers to be returned.
2. The remainder of the list, xs, is the integers from three upwards. They are filtered with the predicate ((/= 0) . (`mod` x)) to create a new list. Remember x is still 2, so this predicate first finds the remainder when the input is divided by 2, and returns True if that remainder is non-zero. As True is returned for all odd elements in xs, the even numbers are filtered out of the list. This filtered list (3, 5, 7, 9...) is passed into primes recursively.
3. This time the first element in the list is 3, so x  is 3, and this becomes the next prime number that is returned. Now xs is the list of odd integers starting at 5. They are filtered as before with x=3 so that the filtered list is the integers from 5 that don't divide evenly by 2 or 3 (5, 7, 11, 13, 17...). This list is passed recursively into primes.
4. This process continues adding another prime number to the list each time.

This process can be illustrated with the following equivalences:

primes [2, 3, 4, 5, 6...] = 2 : primes [3, 5, 7, 9, 11...]

                          = 2 : 3 : primes [5, 7, 11, 13, 17...]

                          = 2 : 3 : 5 : primes [7, 11, 13, 17, 19...]

                          = 2 : 3 : 5 : 7 : primes [11, 13, 17, 19, 23...]

                          = 2 : 3 : 5 : 7 : 11 : primes [13, 17, 19, 23, 29...]

Although not the most efficient prime number generator, I did find this one to be expressive, and a good learning aid.

I'll follow up in future posts with more prime number generators.

Comparison of C++ and Haskell for Project Euler

I've already solved quite a few of the problems on Project Euler using C++, but as I've recently been learning Haskell, I decided to redo them all using it. Solving real problems is a good way to learn a new language and as I've already done them in C++ I can make some comparisons between the two languages.

Haskell for Maths Problems

Without using Haskell for very long it's easy to see that it is well suited to maths problems. This is because the syntax is optimised for expressing mathematical statements and it often resembles maths notation. For example, consider the following definition for the absolute function, absolute:

absolute n
  | n >= 0    = n
  | otherwise = -n

You don't need to know Haskell to figure out what this function does. Another frequently used example is the factorial function. It can be defined, for example, in the following three ways:

factorial n
  | n == 0    = 1
  | otherwise = n * factorial (n-1)

factorial2 0 = 1
factorial2 n = n * factorial2 (n-1)

factorial3 n = product [1..n]

Not only are these more succinct than in a procedural language, but they also read more naturally.

Haskell for Project Euler, Problem 5

Given Haskell's suitability for maths problems, it should be well suited to Project Euler. I'll demonstrate that here using Problem 5, "What is the smallest positive number that is evenly divisible by all of the numbers from 1 to 20?"

Firstly, here is an explanation of my approach. The solution must be a multiple of all the prime numbers below 20, and as they don't share any common factors the solution must be a multiple of the product of the primes below 20. Therefore I calculate this product and test multiples of it until I find the first one that the non primes below 20 also divide evenly into.

Here is the solution implemented in C++

bool isValidAnswer(int numberToCheck) {
  vector<u64> nonPrimesTo20{ 4, 6, 8, 9, 10, 12, 14, 15, 16, 18 };
  auto isMultiple = true;
  for_each(nonPrimesTo20.begin(), nonPrimesTo20.end(), [&](u64 n) {
    if (numberToCheck%n != 0) {
      isMultiple = false;
    }
  });
  return isMultiple;
}

int main(void) {
  const auto productOfPrimesTo20 = 2*3*5*7*11*13*17*19;
  auto numberToCheck = productOfPrimesTo20;

  while (!isValidAnswer(numberToCheck)) {
    numberToCheck += productOfPrimesTo20;
  }
  cout << "Answer is " << numberToCheck << endl;
  return 0;
}

and here is the same algorithm implemented in Haskell:

productOfPrimesTo20 = 2 * 3 * 5 * 7 * 11 * 13 * 17 * 19
nonPrimesTo20 = [4, 6, 8, 9, 10, 12, 14, 15, 16, 18]
isValidAnswer x = and [mod x n == 0 | n <- nonPrimesTo20]
answer = head [x | x <- [productOfPrimesTo20, productOfPrimesTo20*2..], isValidAnswer x]

Again, it is clear that the Haskell solution is much more succinct than in C++. Not only that, but if you understand the Haskell notation, the Haskell solution is also more expressive, and it's easier to read as you don't need to keep track of a loop or any variables.

If you're not familiar with Haskell then this code will look pretty obfuscated, but I've got posts coming up on Haskell prime number generation that explain many of the techniques used above.

Comparison of tuples using variadic templates

In a previous post I explained how the tuple has been implemented in VC++ 2013 using variadic templates. That post gave some explanation of the definition and construction of tuples. You'll should understand that post before you read this one.

The comparison of tuples is interesting and it also uses variadic templates. The STL provides comparison operators for tuples even though they can be composed of many different types, which aren't known, and often don't exist, until the client code is written. This post shows how this is achieved.

Comparison of tuples

The include file <tuple> contains the (non-member) helper operator:

template<class... _TypesLhs, class... _TypesRhs> inline
bool operator==(const tuple<_TypesLhs...>& _Left, const tuple<_TypesRhs...>& _Right)
{
  return (_Left._Equals(_Right));
}

Here _TypesLhs is the types of the tuple on the left hand side and _TypesRhs the types of the tuple on the right. Note that this is not the type of the tuple itself, but the collection of types of all the elements of the tuple. This function simply calls through to the tuple member function _Equals which looks like this:

template<class... TypesRhs>
bool _Equals(const tuple<TypesRhs...>& _Right) const
{
  return (_Myfirst._Val == _Right._Myfirst._Val
  && _Mybase::_Equals(_Right._Get_rest()));
}

Remember that a tuple is composed of a base part that holds a single element of a templated type, and a derived part that holds the remaining elements and the collection of associated types. This pattern continues recursively down to a base class that doesn't hold any value and is based on a template definition with zero types.

The _Equals function above compares the two elements in the most derived class. If they don't match then it will return false. If they do match is will call recursively into the _Equals function in the base class. If these calls make it all the way to the base tuple, which is templated for zero arguments, then it always returns true and the tuples match. This is another example of how powerful variadic templates can be.

You should understand that if the corresponding elements of two tuples can't be compared with the == operator then compiler will complain if you try to compare those tuples. This reminds us how much of the work here is being done at compile time. Template classes and template functions are not classes and functions. They're templates that the compiler uses to write classes and functions during compilation, and this is just the same when variadic templates are being used.

Comparing tuples with different numbers of elements

Comparing tuples of different sizes obviously shouldn't return true, but this is in fact also dealt with at compile time. You'll need to understand a few more things to see how this is done, but they're interesting so I'm adding them into this post.

Firstly, I previously omitted from the _Equals function the following line:

static_assert(_Mysize == sizeof...(_Other), "comparing tuple to object with different size");

The first parameter to a static_assert must be a constant expression that can be converted to a bool. A constant expression is an expression that can be evaluated by the compiler. As I'll show, _Mysize is a const and sizeof can be figured out by the compiler, so _Mysize == sizeof...(_Other) is a constant expression. If the compiler evaluates this expression to true then the static_assert declaration has no effect. If it evaluates to false the compiler gives an error with the message supplied in the second argument to static_assert.

So, how does the compiler determine the size of each tuple? The tuple member variable _Mysize is defined as follows:

static const size_t _Mysize = 1 + sizeof...(_Rest);

Here, _Rest is the collection of types of the tuple from which this class is derived, i.e. all the types not held by this class. To understand how this calculation works consider sizeof(tuple<double, long, std::string>). Here:

sizeof(tuple<double, long, std::string>) = 1 + sizeof(tuple<long, std::string);
sizeof(tuple<long, std::string>) = 1 + sizeof(std::string);
sizeof(tuple<std::string>) = 1 + sizeof(tuple<>);
sizeof(tuple<>) = 0;

This should explain how comparing tuples with differing numbers of elements can give a compilation error.

More on STL Tuples - Implementation using Variadic Templates

My previous post on STL tuples mentioned that they can be implemented using variadic templates. This post looks at the VC++ 2013 tuple implementation to see how this has been achieved. Doing so gives a great example of how powerful variadic templates can be.

(Note, in order to keep the code below clean I've cut out some details, e.g. template parameters with defaults that are used. These extras have no affect on what I'm explaining and clutter up the code snippets. If you want to see these details go look in <tuple>.)

Simple class definition for a tuple in VC++ 2013

Firstly, there is the following class definition for tuple:

template<>
class tuple<>
{
public:
  tuple()
  {
  }
  // more implementation

This is a template class with no template parameters, and a constructor that takes no parameters. You could invoke this constructor with:

tuple<> temp;

but I can't think of a reason to ever do that.

Variadic template class definition for the tuple

There is then another tuple definition that uses variadic templates:

template<class T, class... _Rest>
class tuple<T, _Rest...> : private tuple<_Rest...>
{ // recursive tuple definition
public:
  _Tuple_val<T> _Myfirst;
  tuple<_Rest...> _Get_rest();

To understand this you'll need to know that the ... notation can be read as "collection of". Also, if you're not familiar with the class keyword inside the template angle brackets, it's equivalent to typename in this context.

Now you can begin to parse the above code. This definition of tuple is templated for a parameter of type T and a collection of types referred to as _Rest. It privately inherits from a tuple that is templated for this collection of types. The comment "recursive tuple definition" is not mine, it's taken from the VC++ 2013 STL code. It explains that when using the above definition the compiler will recursively create new tuple definitions for a shrinking list of types until that list is empty, at which point it will use the original tuple definition, which is templated for zero types.

To explain with an example, consider you wish to use tuple<double, long, std::string>. The compiler will first create a class definition of tuple that is templated for type double and the collection of types long and std::string. This class will inherit from another class definition of tuple that is templated for type long  and the collection of types std::string. This class will inherit from a tuple class that is templated for type std::string and inherits from the original tuple class definition that is templated for zero types.

Construction and Implementation details

Look at the two members that I've picked out. _Tuple_val<T> is a helper class and it's used to store a value of type T. _Get_rest() is a member function that returns a reference to the base class. So, from the previous example, in the most derived class, which is templated for type double and the collection of types long and std::string, _Myfirst refers to the value of the double and _Get_rest() returns a reference to its base class i.e. the class that is templated for a long  and the collection of types std::string.

Now look at the constructor that is used when many types are passed in:

explicit tuple(T&& _This_arg, _Rest&&... _Rest_arg)
: _Myfirst(_This_arg), tuple<_Rest...>(_Rest_arg...)
{
}

As you can see, the constructor simply initializes the _Myfirst variable with the corresponding value _This_arg, and the rest of the arguments are passed to the base class, which is by definition a tuple templated for the types of those arguments.

Hopefully this example gives a good idea of the power of variadic templates. I'll follow up with another post on some more VC++ 2013 tuple implementation details including how they are compared.

STL Tuples

The STL tuple class is a bit like the STL pair, except the tuple supports any number of types to be tied together in a single object. This can be achieved using variadic templates, although many tuple implementations were completed before variadic templates were available (e.g. VC++ 2012). They implemented the tuple by using different macros for different numbers of arguments.

STL pairs are useful but sometimes not sufficient

STL pairs are useful when two pieces of data are always coupled together, for example:

std::pair<std::string, float> raceWinner{ "M40", 17.20 };

However, in the past, to include any more pieces of data required some more work, for example creating a simple struct or class:

struct triathlonResult {
std::string ageCategory;
float swim;
float cycle;
float run;
};

which (if uniform initialization is supported) can then be used with:

triathlonResult triathlonWinner{ "M35", 13.45, 35.5, 17.3 };

STL tuples extend the pair to any number of arguments

The introduction of the tuple to the standard library means that a custom class is no longer needed, instead:

std::tuple<std::string, float, float, float> triathlonWinner( "F35", 13.45, 35.5, 17.3 );

can be used to create an object that holds the same information.

Initialization can be simplified using auto and make_tuple like this:

auto triathlonWinner = make_tuple("F35", 13.45, 35.5, 17.3);

An added convenience is that the comparison operators (==, !=, >, >=, < and <=) are included in the standard library. From the standard [tuple.rel] - "The elementary comparisons are performed in order from the zeroth index upwards. No comparisons or element accesses are performed after the first equality comparison that evaluates to false".

Retrieving elements from a tuple

Retrieving elements from a tuple may differ from how you're used to getting them from a pair. One way to retrieve members from a pair is to use the public members first and second. For example:

auto raceWinnersCategory = raceWinner.first;
auto raceWinnersTime = raceWinner.second;

However, there is no tuple member function to retrieve elements, instead they must be retrieved using the non-member function get (defined in <tuple>). The reason is that the element to retrieve must be indexed, and if this were to happen on a member function it would need to be invoked using the template keyword like this:

auto triWinnersCategory = triathlonWinner.template get<0>();
(for more information see [tuple.helper] note 11 from the standard)

To avoid this a non-member get function is provided, to be used as shown:
auto triWinnersCategory = get<0>(triathlonWinner);

auto triWinnersSwimTime = get<1>(triathlonWinner);

auto triWinnersCycleTime = get<2>(triathlonWinner);

auto triWinnersRunTime = get<3>(triathlonWinner);

This may not seem very natural to code at first but it is easy to read. Note that this get function can also be used for pairs, so if you suspect your pairs may one day be promoted to tuples then it's a good idea to use get with pairs too.

Uniform Initialization in VC++ 2013

Before uniform initialization was introduced to C++ there were many ways to initialize variables. For example:

int myInt(3);

int myInt = 3;

Musician guitarist; // Zero argument constructor

Musician drummer("Black", 10.01, true);

Uniform initialization provides a uniform way to initialize objects using curly braces. For example:

int myInt {4};

Musician piper {}; // Zero argument constructor

Musician jazzFlautist {"Black", 12.34, true};

Uniform initialization provides a standard method for initializing objects and for consistency alone should be the first preference. However, there are further advantages.

Further advantages of Uniform Initialization

The uniform initialization syntax can be used with the class std::initializer_list<> to pass a variable number of similar arguments to a constructor. This is useful for creating std containers, for example:

vector<int> myInts { 1, 3, 2, 4};

will create a vector of ints with four elements: 1, 2, 3, and 4. To make use of std::initializer_list<> in your own classes you'll need to write an initializer list constructor, i.e. a constructor that takes an std::initializer_list<> of the desired type and parses it as required. If you do this users will be able to use the { } notation to invoke that constructor.

Uniform initialization can also be nested, for example:

vector<Musician> band { guitarist ,

                        jazzFlautist,

                        {"Blue", 43.21, false} };

In the above case guitarist and jazzFlautist are added to the vector, as is a third anonymous object that is created and copied (or hopefully moved) into the vector.

Uniform initialization also brings a way to avoid Most Vexing Parse confusion. For a simple illustration consider the following code:

Musician harpist();

This may be considered an acceptable way to create a Musician object called harpist with the zero argument constructor, but in fact it declares a function called harpist that takes zero arguments and returns a Musician. This ambiguity can take more subtle forms, for example where a function call is intended as an input parameter to a function call, but the compiler reads this as a function declaration that takes a function type as a parameter. It's so called the most vexing parse for good reason, but it is avoided using by the curly braces uniform initialization notation, i.e.

Musician harpist{};

Another advantage of uniform initialization is that it is a convenient way to initialise Aggregate classes. For example given the following definition:

struct number {
int asInteger;
std::string asString;
bool isPrime;
};

a number can be initialized with:

number theNumberFour{ 4, string("four"), false };

Reasons to still use the old constructor syntax

There may still be times when you wish to use the old constructor syntax. If a class has an initializer list constructor then it will take precedence over other constructors. If you wish to call the other constructors you'll need to use the old constructor syntax to be explicit.

For example, if you wish to create a std::vector of ints and you wish to call the constructor that takes a single parameter specifying the vector's size then writing:

std::vector<int> myInts{4};

won't work. The above code would create a vector of ints and initialize it with the single element 4. Instead you need to use the round bracket notation :

std::vector<int> myInts(4);

However, in general you should default to using uniform initialization, and only break that pattern when you've got a specific reason.

Writing an Arbitrarily Large Integer Class in C++

Project Euler Problem 25 is to find the first 1000 digit Fibonacci number. I'm implementing a solution in C++ so the problem here is that the largest number that can be represented with a built-in type is 2⁶⁴. A few of the Project Euler problems prior to this one required handling numbers greater than 2⁶⁴ - for example Problem 13 is to find the 10 least significant digits of a sum of one-hundred 50-digit numbers - but in those cases it was possible to work around the 2⁶⁴ limitation by discarding digits that aren't needed in the final answer or using doubles.

However, this problem prompted me to write my own class that could deal with arbitrarily large integers. To this end I have created embiguint.

Main goals of EMBigUInt

The main goals of the class are as follows:

• The class should work naturally with the build in types. For example it should support both:
• embiguint = int + embiguint, and
• embiguint = embituint + int
• The class should provide implicit construction from the built in types.
• Like the built-in types, the class should work with all the standard arithmetic operators (+, -, /, *, +=, *=, etc...) and comparison operators (==, >, >=, !=, etc...). These operators should work with embiguints and between embiguints and the built-in types.
• The class should provide some additional functionality beyond what the built-in types, for example:
• Construction from a string of digits to support constructing an embiguint greater than 2⁶⁴ e.g. 
  embiguint myBigUInt("1234567890123456789012")
• Return a count of the number of digits.

There's still much more work to do to make the class I want, but here are some notes on design decisions I've taken so far:

• Efficiency is not my first concern. My initial aim it to get a working class. When I've got all the functionality I need with unit tests and performance tests in place I can work on efficiency in the confidence that I won't break existing functionality and I can measure any speed gains. Some outcomes from this are:
• I'm using standard library vectors, which are easy to use, provide lots of debugging help, and come with lots of standard library algorithms. It's worth investigating if raw arrays could improve efficiency, but for now this would slow down development.
• I'm relying on implicit conversion from built-in types to embiguint for some of the functions. I could implement more efficient functions if I overloaded specifically for the built-in types but, again, this can be done as a refactoring.
• So far, most of the functions in the library are non-member functions. This follows from the advise in Item 24 of Scott Meyers' Effective C++: Declare non-member functions when type conversion should apply to all parameters, and it allows me to do
  embiguint = int + embiguint
  Here I haven't overloaded operator+ for the type int, instead the non-member function 
  const embiguint operator+ (const embiguint& lhs, const embiguint& rhs)
  is being used with implicit type conversion converting the int to an embiguint.

Using EMBigUInt for Project Euler Problem 25

This class in its current state enables Project Euler Problem 25 to be solved very simply with the following:

embiguint lastValue = 0;
embiguint currentValue = 1;
emuint counter = 1;

while (currentValue.GetNoDigits() < 1000) {
  embiguint temp = currentValue;
  currentValue += lastValue;
  lastValue = temp;
  counter++;
}

cout << "Answer is " << counter << endl;