Overview

A few short algorithms play essential roles in multiple Project Euler problems.
Instead of writing them from scratch over and over again, I copy'n'paste'n'modify where appropriate.

My goal is that all solutions include only standard C++ library headers (or C++11/14/17/you-know-the-game).
They are "single-file-only" programs without external dependencies.
Therefore the snippets shown below are NOT put into a dedicated toolbox.h file.

List of algorithms:

Greatest Common Divisor

Least Common Multiple
Coprime Test

Prime Sieves

Sieve Of Eratosthenes
Trial Division
Precomputed Tables

Primality Tests

Trial Division
Wheel Factorization
Miller-Rabin Test
Fermat Test

Prime Counting

Binary Search
Legendre Algorithm
Meissel-Lehmer Algorithm

Modular Arithmetic

Modular Multiplication
Modular Exponentiation
Modular Inverse

Combinatorics

Binomial Coefficient

Pseudo-Random Numbers

Linear Congruential Generator

Working With Digits

Digital Sum
Fingerprint

Big Integers

Fractions

Greatest Common Divisor

       // greatest common divisor
       template <typename T>
       T gcd(T a, T b)
       {
         while (a != 0)
         {
           T c = a;
           a = b % a;
           b = c;
         }
         return b;
       }

The Euclidean algorithm is a very efficient, simple (and surprisingly old !) algorithm. It's younger sibling, the binary GCD algorithm is supposed to be a tad bit faster, though.

Note: On newer versions of the GCC compiler there is a __gcd(a, b) function in #include <stl_algo.h> which might be faster on your computer.
(it's very likely that std::gcd will be added to the C++17 standard)

Modern CPU's have a few instructions that can accelerate the GCD computation by about 30%. However, accessing these instructions isn't portable.
An implementation for the GCC compiler looks like this (with automatic fallback to default algorithm on other compilers):

       // greatest common divisor
       template <typename T>
       T gcd(T a, T b)
       {
       #ifdef __GNUC__
         // adopted from https://lemire.me/blog/2013/12/26/fastest-way-to-compute-the-greatest-common-divisor/
         // and          https://github.com/lemire/Code-used-on-Daniel-Lemire-s-blog/blob/master/2013/12/26/gcd.cpp
         if (a == 0)
           return b;
         if (b == 0)
           return a;
        
         // MSVC++: _BitScanForward intrinsic instead
         auto shift = __builtin_ctz(a | b);
         a >>= __builtin_ctz(a);
         do
         {
           b >>= __builtin_ctz(b);
           if (a > b)
             std::swap(a, b);
        
           b -= a;
         } while (b != 0);
         return a << shift;
        
       #else
        
         // standard GCD
         while (a != 0)
         {
           T c = a;
           a = b % a;
           b = c;
         }
         return b;
       #endif
       }

Least Common Multiple

       // least common multiple
       template <typename T>
       T lcm(T a, T b)
       {
         // parentheses avoid overflows for certain input values
         return a * (b / gcd(a, b));
       }

The Least Common Multiple builds upon the Greatest Common Divisor.

Coprime Test

       // return true if a and b are coprime
       // note: needs an implementation of gcd()
       template <typename T>
       bool isCoprime(T a, T b)
       {
         // fast reject if both are even (=> gcd(a,b) >= 2)
         if (((a|b) & 1) == 0)
           return false;
        
         return gcd(a, b) == 1;
       }

Prime Sieves

Sieve Of Eratosthenes

Advantages:

faster than other simple sieves
ok-ish memory consumption

Disadvantages:

optimized version needs a dedicated lookup-function

       #include <vector>
        
       // odd prime numbers are marked as "true" in a bitvector
       std::vector<bool> sieve;
        
       // return true, if x is a prime number
       bool isPrime(unsigned int x)
       {
         // handle even numbers
         if ((x & 1) == 0)
           return x == 2;
        
         // lookup for odd numbers
         return sieve[x >> 1];
       }
        
       // find all prime numbers from 2 to size
       void fillSieve(unsigned int size)
       {
         // store only odd numbers
         const unsigned int half = (size >> 1) + 1;
        
         // allocate memory
         sieve.resize(half, true);
         // 1 is not a prime number
         sieve[0] = false;
        
         // process all relevant prime factors
         for (unsigned int i = 1; 2*i*i < half; i++)
           // do we have a prime factor ?
           if (sieve[i])
           {
             // mark all its multiples as false
             unsigned int current = 3*i+1;
             while (current < half)
             {
               sieve[current] = false;
               current += 2*i+1;
             }
           }
       }

My optimized code needs 0.5 bits per number which is less than what is required for the sieve based on trial division for most use cases.
All prime numbers from 2 to 100 million can be found in about one third of a second and need about 6 MByte.
A multi-threaded version is much faster, though.

Note: I wrote a segmented prime sieve in problem 196.

Based On Trial Division

Advantages:

simple code
can be computed in an incremental way ("on demand")

Disadvantages:

slow division
not very memory-efficient (allocation is slow, too !)

       std::set<unsigned int> primes = { 2 };
       for (unsigned int i = 3; i <= limit; i += 2)
       {
         bool isPrime = true;
        
         // test against all prime numbers we have so far (in ascending order)
         for (auto p : primes)
         {
           // next prime is too large to be a divisor ?
           if (p*p > i)
             break;
        
           // divisible ? => not prime
           if (i % p == 0)
           {
             isPrime = false;
             break;
           }
         }
        
         // yes, we have a prime number
         if (isPrime)
           primes.insert(i);
       }

While std::set offers some nice features "for free", such as looking up whether it contains an element, its memory efficiency can be sub-prime. In those cases I replace it by std::vector and change both inserts to push_back.

Sometimes std::unordered_set might be an option, too.

Precomputed Tables

If code size and/or compilation times are no issue, then a precomputed table of prime numbers is certainly the fastest approach:

       const unsigned int primes[] = {
       #include "primes1000000.txt"
       };

You can download these files from my server:

all primes from 2 to 100 (0.1k)
all primes from 2 to 1000 (0.6k)
all primes from 2 to 10000 (5.8k)
all primes from 2 to 100000 (54.8k)
all primes from 2 to 1000000 (525.8k)

All files contain a comma-separated list of integers, e.g. primes100.txt looks like this:

       2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79,83,89,97

Primality Tests

Trial Division

This test is basically the same used in the Trial Division Prime Sieve, except that we don't know all smaller primes and therefore have to check against all odd numbers.

       bool isPrime(unsigned int x)
       {
         // even number ? only even prime is 2
         if ((x & 1) == 0)
           return x == 2;
        
         // check all even numbers from 3 to sqrt(x)
         for (unsigned int i = 3; i*i <= x; i += 2)
           if (x % i == 0)
             return false;
        
         // passed all tests, must be a prime number
         return x > 1;
       }

Wheel Factorization

Wheel Factorization reduces the number of divisions by observing that all prime factors can be expressed as:
2, 3, 5 or any m: 30m+1, 30m+7, 30m+11, 30m+13, 30m+17, 30m+19, 30m+23, 30m+29 (wheel size 30=2*3*5).
Some of these numbers are actually composite numbers (and not true prime factors) and therefore cause redundant work - but at least we don't forget about any prime factor.
More formally: this set of numbers contains all prime factors but is not minimal.

Note: running a Wheel-based primality test on all numbers from 1 to 10 million is roughly two times faster than with Trial Division
(1.8 vs. 3.3 seconds ... and only 1.4 seconds for the Miller-Rabin test).

       bool isPrime(unsigned long long x)
       {
         // prime test for 2, 3 and 5 and their multiples
         if (x % 2 == 0 || x % 3 == 0 || x % 5 == 0)
           return x == 2 || x == 3 || x == 5;
        
         // wheel with size 30 (=2*3*5):
         // test against 30m+1, 30m+7, 30m+11, 30m+13, 30m+17, 30m+19, 30m+23, 30m+29
         // their deltas/increments are:
         const unsigned int Delta[] = { 6, 4, 2, 4, 2, 4, 6, 2 };
         // start with 7, which is 30*0+7
         unsigned long long i = 7;
         // 7 belongs to the second test group
         auto pos = 1;
        
         // check numbers up to sqrt(x)
         while (i*i <= x)
         {
           // not prime ?
           if (x % i == 0)
             return false;
        
           // skip forward to next test divisor
           i += Delta[pos];
           // next delta/increment
           pos = (pos + 1) & 7;
         }
        
         // passed all tests, must be a prime number
         return x > 1;
       }

Miller-Rabin Test

The Miller-Rabin test is a probabilistic primality test. Using specific witnesses it can be made a deterministic test for all 64 bit numbers.
My code needs the mulmod and powmod function also found on this page (see here).

       // Miller-Rabin-test
       bool isPrime(unsigned long long p)
       {
         // IMPORTANT: requires mulmod(a, b, modulo) and powmod(base, exponent, modulo)
        
         // some code from             https://ronzii.wordpress.com/2012/03/04/miller-rabin-primality-test/
         // with optimizations from    http://ceur-ws.org/Vol-1326/020-Forisek.pdf
         // good bases can be found at http://miller-rabin.appspot.com/
        
         // trivial cases
         const unsigned int bitmaskPrimes2to31 = (1 <<  2) | (1 <<  3) | (1 <<  5) | (1 <<  7) |
                                                 (1 << 11) | (1 << 13) | (1 << 17) | (1 << 19) |
                                                 (1 << 23) | (1 << 29); // = 0x208A28Ac
         if (p < 31)
           return (bitmaskPrimes2to31 & (1 << p)) != 0;
        
         if (p %  2 == 0 || p %  3 == 0 || p %  5 == 0 || p % 7 == 0 || // divisible by a small prime
             p % 11 == 0 || p % 13 == 0 || p % 17 == 0)
           return false;
        
         if (p < 17*19) // we filtered all composite numbers < 17*19, all others below 17*19 must be prime
           return true;
        
         // test p against those numbers ("witnesses")
         // good bases can be found at http://miller-rabin.appspot.com/
         const unsigned int STOP = 0;
         const unsigned int TestAgainst1[] = { 377687, STOP };
         const unsigned int TestAgainst2[] = { 31, 73, STOP };
         const unsigned int TestAgainst3[] = { 2, 7, 61, STOP };
         // first three sequences are good up to 2^32
         const unsigned int TestAgainst4[] = { 2, 13, 23, 1662803, STOP };
         const unsigned int TestAgainst7[] = { 2, 325, 9375, 28178, 450775, 9780504, 1795265022, STOP };
        
         // good up to 2^64
         const unsigned int* testAgainst = TestAgainst7;
         // use less tests if feasible
         if (p < 5329)
           testAgainst = TestAgainst1;
         else if (p < 9080191)
           testAgainst = TestAgainst2;
         else if (p < 4759123141ULL)
           testAgainst = TestAgainst3;
         else if (p < 1122004669633ULL)
           testAgainst = TestAgainst4;
        
         // find p - 1 = d * 2^j
         auto d = p - 1;
         d >>= 1;
         unsigned int shift = 0;
         while ((d & 1) == 0)
         {
           shift++;
           d >>= 1;
         }
        
         // test p against all bases
         do
         {
           auto x = powmod(*testAgainst++, d, p);
           // is test^d % p == 1 or -1 ?
           if (x == 1 || x == p - 1)
             continue;
        
           // now either prime or a strong pseudo-prime
           // check test^(d*2^r) for 0 <= r < shift
           bool maybePrime = false;
           for (unsigned int r = 0; r < shift; r++)
           {
             // x = x^2 % p
             // (initial x was test^d)
             x = mulmod(x, x, p);
             // x % p == 1 => not prime
             if (x == 1)
               return false;
        
             // x % p == -1 => prime or an even stronger pseudo-prime
             if (x == p - 1)
             {
               // next iteration
               maybePrime = true;
               break;
             }
           }
        
           // not prime
           if (!maybePrime)
             return false;
         } while (*testAgainst != STOP);
        
         // prime
         return true;
       }

Fermat Test

The Fermat primality test is neither fast nor accurate - but its implementation is quite short.
My code needs the mulmod and powmod function also found on this page (see here).

       // return true if p is probably prime, accuracy depends on the number of iterations
       // unfortunately, so-called Carmichael numbers pass all tests, see https://en.wikipedia.org/wiki/Carmichael_number
       bool isPrime(unsigned long long p, unsigned int iterations = 90)
       {
         // fast check against small primes, eliminates a few Carmichael numbers, too
         if (p < 2 || p % 2 == 0 || p % 3 == 0 || p % 5 == 0 || p % 7 == 0 || p % 11 == 0)
           return p == 2 || p == 3 || p == 5 || p == 7 || p == 11;
        
         // just a pseudo-random seed
         unsigned long long a = p / 2 + 1;
        
         // check p against several pseudo-random numbers
         for (unsigned int i = 0; i < iterations; i++)
         {
           // we test against random numbers between 1 and p - 1 (inclusive)
           // let's compute a simple hash: (i * 2654435761) % 2^32
           // (Knuth's multiplicative hash)
           a *= 2654435761;
           if (a == 0 || a >= p)
             a = a % (p - 1) + 1;
        
           // Fermat found that for each prime: a^(p-1) % p == 1
           if (powmod(a, p-1, p) != 1)
             return false; // failed test, p is definitely not prime ...
         }
        
         // p is probably prime
         return true;
       }

Prime Counting

The source code for all three algorithm can be found below.

Binary Search

If the number of primes is small (a few million) then a standard prime sieve can precompute all prime numbers.
When those primes are stored in a std::vector named primes then a basic bisection is extremely fast for an existing prime x:

       std::distance(primes.begin(), std::upper_bound(primes.begin(), primes.end(), x));

This code is actually part of the next two algorithms, too (see the first lines of the pi function).

Legendre Algorithm

Two functions pi and phi are essential for the Legendre algorithm:
pi(x) counts all prime numbers up to (and including) x.
phi(x, a) counts all numbers up to (and including) x that are coprime to the first a primes.
The only difference between the Legendre and the faster Meissel-Lehmer algorithm is the implementation of pi.

Note: my version of phi contains the optimizations of Mapes and I added a cache for small values (phicache).
The code works perfectly fine without the neat tricks - just slower.

Meissel-Lehmer Algorithm

The Meissel-Lehmer algorithm has a more complex pi function which puts less work on phi.

       #include <vector>
       #include <algorithm>
       #include <cmath>
        
       // uncomment to run (slower) Legendre algorithm
       //#define LEGENDRE
        
       // cache small value of phi(x,a)
       static std::vector<std::vector<unsigned short>> phicache;
        
       // forward declaration
       unsigned long long pi(unsigned long long x);
        
       // return how many numbers are coprime to the first a primes
       long long phi(unsigned long long x, unsigned int a, int sign = +1)
       {
         // Mapes equations for small a
         if (a == 0)
           return sign * x;
         if (a == 1)
           return sign * (x - x / 2);
         if (a == 2)
           return sign * (x - x / 2 - x / 3 + x / 6);
        
         // precomputed tables for small x
         if (x < phicache.size() && a < phicache[x].size())
           return sign * phicache[x][a];
        
       #ifndef LEGENDRE
         // http://www.cnblogs.com/boceng/p/7222701.html
         if (x < primes[a] * primes[a])
           return sign * ((int)pi(x) - (int)a + 1);
       #endif
        
         // partially iterative algorithm
         // most people use the slower version: phi(x, a-1) - phi(x/p[a], a-1)
         unsigned long long result = 0;
         while (a > 3) // a = 0..3 handled by Mapes equations
         {
           if (a < primes.size() && x < primes[a]) // check for primes.size() only needed for Legendre
             return result + sign;
        
           a--;
           result += phi(x / primes[a], a, -sign);
         }
        
         // Mapes a = 3
         return result + sign * (x - x / 2 - x / 3 + x / 6 - x / 5 + x / 10 + x / 15 - x / 30);
       }
        
       // two implementations:
       // - Legendre: simpler but slower ==> pi(10^11) in 11 seconds
       // - Lehmer:   more complex but faster => pi(10^11) in about 0.7 seconds
       // both numbers including the overhead of the prime sieve (up to sqrt(10^11)) and caching pi up to 100 and phi up to (1000,168)
       unsigned long long pi(unsigned long long x)
       {
         // small x: fast binary search in list of primes
         if (x < primes.back())
         {
           auto nextPrime = std::upper_bound(primes.begin(), primes.end(), x);
           return std::distance(primes.begin(), nextPrime);
         }
        
       #ifdef  LEGENDRE
         // Legendre
         auto a = pi(sqrt(x));
         return phi(x, a) + a - 1;
       #else
         // Meissel-Lehmer
         // good overview: http://acganesh.com/blog/2016/12/23/prime-counting
         // Lehmer paper:  https://projecteuclid.org/download/pdf_1/euclid.ijm/1255455259
        
         auto a = pi(pow(x, 0.25));
         auto b = pi(sqrt(x));
         auto c = pi(cbrt(x));
        
         auto result = phi(x, a) + (b + a - 2) * (b - a + 1) / 2;
         for (auto i = a; i < b; i++)
         {
           auto xDivP = x / primes[i];
           result -= pi(xDivP);
        
           if (i < c)
           {
             auto b_i = pi(sqrt(xDivP));
             for (auto j = i; j < b_i; j++)
               result -= pi(xDivP / primes[j]) - j;
           }
         }
         return result;
       #endif
       }

Benchmark

The function countPrimes in problem 501 does the same job - and is much faster.
However, I didn't write it (only modified it a little bit) and still have no complete understanding of its inner workings despite Muthu Veerappan's explanation on this blog.

pi(10¹¹) = 4118054813
Legendre 11.2 seconds
Meissel-Lehmer 0.7 seconds
PE501 0.25 seconds

pi(10¹²) = 37607912018
Legendre 103 seconds,
Meissel-Lehmer 4.3 seconds
PE501 1.2 seconds

Modular Arithmetic

Modular Multiplication

There are multiple versions of my modular arithmetic multiplication code (mulmod) which solves (a*b)%modulo.
The 32 bit version is extremely short:

       // return (a*b) % modulo
       unsigned int mulmod(unsigned int a, unsigned int b, unsigned int modulo)
       {
         return ((unsigned long long)a * b) % modulo;
       }

The story is completely different for 64 bit because of a lack of standardized 128 bit integers in C/C++.
Bitwise multiplication is very similar to multiplication taught in school (and works only up to 63 bits):

       // return (a*b) % modulo
       unsigned long long mulmod(unsigned long long a, unsigned long long b, unsigned long long modulo)
       {
         // (a * b) % modulo = (a % modulo) * (b % modulo) % modulo
         a %= modulo;
         b %= modulo;
        
         // fast path
         if (a <= 0xFFFFFFF && b <= 0xFFFFFFF)
           return (a * b) % modulo;
        
         // we might encounter overflows (slow path)
         // the number of loops depends on b, therefore try to minimize b
         if (b > a)
           std::swap(a, b);
        
         // bitwise multiplication
         unsigned long long result = 0;
         while (a > 0 && b > 0)
         {
           // b is odd ? a*b = a + a*(b-1)
           if (b & 1)
           {
             result += a;
             if (result >= modulo)
               result -= modulo;
             // skip b-- because the bit-shift at the end will remove the lowest bit anyway
           }
        
           // b is even ? a*b = (2*a)*(b/2)
           a <<= 1;
           if (a >= modulo)
             a -= modulo;
        
           // next bit
           b >>= 1;
         }
        
         return result;
       }

The blockwise algorithm is an extension of the bitwise algorithm and multiplies multiple bits at once.
The 63-bit limit applies as well:

       // return (a*b) % modulo
       unsigned long long mulmod(unsigned long long a, unsigned long long b, unsigned long long modulo)
       {
         // (a * b) % modulo = (a % modulo) * (b % modulo) % modulo
         a %= modulo;
         b %= modulo;
        
         // fast path
         if (a <= 0xFFFFFFF && b <= 0xFFFFFFF)
           return (a * b) % modulo;
        
         // count leading zero bits of modulo
         unsigned int leadingZeroes = 0;
         unsigned long long m = modulo;
         while ((m & 0x8000000000000000ULL) == 0)
         {
           leadingZeroes++;
           m <<= 1;
         }
        
         // cover all bits of modulo
         unsigned long long mask = (1 << leadingZeroes) - 1;
        
         // bitwise multiplication
         unsigned long long result = 0;
         while (a > 0 && b > 0)
         {
           result += (b & mask) * a;
           result %= modulo;
        
           // next bits
           b >>= leadingZeroes;
           a <<= leadingZeroes;
           a  %= modulo;
         }
         return result;
       }

A GCC-only version of mulmod can be considerably faster:

       unsigned long long mulmod(unsigned long long a, unsigned long long b, unsigned long long modulo)
       {
         // note: you can remove the following two lines
         // but if (a*b)/c > 2^64 then the Assembler code will trigger a floating-point exception
         a %= modulo;
         b %= modulo;
        
         // fast path
         if (a <= 0xFFFFFFF && b <= 0xFFFFFFF)
           return (a * b) % modulo;
        
         // perform 128 bit arithmetic
        
       #ifdef __x86_64__
         // use GCC inline assembler
         unsigned long long asmResult;
         __asm__
         (
           "mulq %2\n" // result in rdx:rax
           "divq %3"   // quotient in rax, remainder in rdx
           : "=&d" (asmResult), "+%a" (a)
           : "rm" (b), "rm" (modulo)
           : "cc"      // clear conditions
         );
         return asmResult;
        
       #else
        
         // based on GCC's 128 bit implementation
         return ((unsigned __int128)a * b) % modulo;
       #endif
       }

Modular Exponentiation

powmod solves (base^exponent)%modulo. My algorithm is based on Exponentiation by squaring.
The 32-bit version doesn't have any dependencies:

       // return (base^exponent) % modulo for 32-bit values, no need for mulmod
       unsigned int powmod(unsigned int base, unsigned int exponent, unsigned int modulo)
       {
         unsigned int result = 1;
         while (exponent > 0)
         {
           // fast exponentation:
           // odd exponent ? a^b = a*a^(b-1)
           if (exponent & 1)
             result = (result * (unsigned long long)base) % modulo;
        
           // even exponent ? a^b = (a*a)^(b/2)
           base = (base * (unsigned long long)base) % modulo;
           exponent >>= 1;
         }
         return result;
       }

The 64-bit version needs mulmod (see above):

       // return (base^exponent) % modulo
       unsigned long long powmod(unsigned long long base, unsigned long long exponent, unsigned long long modulo)
       {
         unsigned long long result = 1;
         while (exponent > 0)
         {
           // fast exponentation:
           // odd exponent ? a^b = a*a^(b-1)
           if (exponent & 1)
             result = mulmod(result, base, modulo);
        
           // even exponent ? a^b = (a*a)^(b/2)
           base = mulmod(base, base, modulo);
           exponent >>= 1;
         }
         return result;
       }

Modular Inverse

The Extended Euclidean algorithm is the fastest way and comes with the least restrictions. It's code looks a bit messy, though.
If the modulo is a prime number then thanks to Fermat's little theorem I can find the modular inverse with powmod:

       // extended Euclidean algorithm, see https://en.wikipedia.org/wiki/Extended_Euclidean_algorithm
       unsigned int modularInverse(unsigned int a, unsigned int modulo)
       {
         // pretty much the same code can be found on https://rosettacode.org/wiki/Modular_inverse
         auto originalModulo = modulo;
         // note: s and t can be negative inside the loop
         int s = 0;
         int t = 1;
         while (a > 1)
         {
           auto tmp      =     modulo;
           auto quotient = a / modulo;
           modulo        = a % modulo;
           a = tmp;
        
           auto tmp2 = s;
           s = t - quotient * s;
           t = tmp2;
         }
        
         // avoid negative result
         return t < 0 ? t + originalModulo : t;
       }
        
       // Fermat's little theorem, see https://en.wikipedia.org/wiki/Fermat%27s_little_theorem
       // modulo needs to be a prime number
       unsigned int modularInverseOfPrime(unsigned int a, unsigned int modulo)
       {
         // Fermat's Little Theorem: a^p     == a mod p
         // divide both sides by a:  a^(p-1) == 1 mod p
         // once more:               a^(p-2) == a^-1 mod p
         // where a^-1 is the inverse of a => equal to a^(p-2) % p
         return powmod(a, modulo - 2, modulo);
       }

Pseudo-Random Numbers

For most of the statistics / probability problems I write a simple Monte-Carlo simulation in order to estimate the result (to avoid that my "proper" algorithm doesn't go wild).
The implementation of the rand() function from the C standard library varies among compilers / platforms.
Therefore I prefer this simple code to reproduce exactly the same result on different machines:

Linear Congruential Generator

The constants were taken from the Wikipedia page.

       // a simple pseudo-random number generator
       // (produces the same result no matter what compiler you have - unlike rand() from math.h)
       unsigned int myrand()
       {
         static unsigned long long seed = 0;
         seed = 6364136223846793005ULL * seed + 1;
         return (unsigned int)(seed >> 30);
       }

The output isn't on a cryptographically safe level but more than sufficient for the Project Euler problems. And pretty fast !

Working With Digits

Digital Sum

Efficient routines to analyze the digits of a number are part of many problems.
The most frequent is the digital sum (sum of all digits):

       // sum of all digits
       unsigned int digitSum(unsigned long long x)
       {
         unsigned int result = 0;
         while (x > 0)
         {
           result += x % 10;
           x      /= 10;
         }
         return result;
       }
        
       // slightly more efficient (saves one iteration)
       unsigned int digitSumFaster(unsigned long long x)
       {
         unsigned int result = 0;
         while (x >= 10)
         {
           result += x % 10;
           x      /= 10;
         }
         return result + x;
       }

The second version is slightly more efficient but in most cases you won't notice a significant difference.
I tried std::lldiv (from the C++11 standard) instead of modulo/division but it's much slower with my G++ compiler.

Fingerprint

In order to find out whether two numbers are permutations of each other, I generate their "fingerprint":
I count how often each digit can be found and return a number where the lowest 5 bits represent the amount of zeros,
the next 5 bits the amount of ones, then twos, etc.
Note: the term "fingerprint" isn't something official - I just like the word.

       // count digits, two numbers share the same fingerprint if they are permutations of each other
       unsigned long long fingerprint(unsigned long long x)
       {
         unsigned long long result = 0;
         while (x > 0)
         {
           auto digit = x % 10;
           x /= 10;
        
           result += 1ULL << (5 * digit);
         }
         return result;
       }

Combinatorics

Binomial Coefficient

The formula for a binomial coefficient is pretty simple but easily causes overflows due to the factorial function.
The following code can compute most of the values if the result fits in 64 bits (even though the factorials may exceed 64 bits):

       // number of ways to choose n elements from k available
       unsigned long long choose(unsigned int n, unsigned int k)
       {
         // n! / (n-k)!k!
         unsigned long long result = 1;
         // choose(n, 0) = 1
         // choose(n, k) = choose(n, k - 1) * (n + 1 - k) / k
         // reduce overflow by dividing as soon as possible to keep numbers small
         for (unsigned int invK = 1; invK <= k; invK++)
         {
           result *= n;
           result /= invK;
           n--;
         }
         return result;
       }

Big Integers

Working with very large numbers can be hard in C++ (whereas it's extremely simple in Java and Python).
Of course, there are dedicated C/C++ libraries, such as GMP, but I don't want to link/include external libraries.

Addition, subtraction and multiplication are straightforward implementations of basic school algorithms.
MaxDigit can be set to other powers-of-10 in order to speed up the code but then the std::string functions will fail.

Fractions

The C++ STL class std::ratio supports only compile-time arithmetic.
Therefore I had to write my class Fraction which was tuned for performance:

it doesn't automatically reduce fractions (to become proper fractions), you have to explicitly call reduced()
there are no checks whether the denominator is zero nor are the signs of numerator / denominator normalized
overflows do not cause exceptions and will produce undefined results
adjust the data type T to fit your needs, e.g. reduced() is much slower for long long then int

       // simple class to representation a fraction
       // note: no checks for a zero denominator
       //       signs are not normalized
       struct Fraction
       {
         // change to long long if you need to support larger values
         typedef int T;
        
         // numerator
         T num;
         // denominator
         T den;
        
         Fraction(T numerator, T denominator = 1)
         : num(numerator), den(denominator)
         {}
        
         // add
         Fraction operator+(const Fraction& other) const
         {
           if (den == other.den)
             return Fraction(num + other.num, den);
           // n1/d1 + n2/d2 = (n1*d2 + n2*d1) / d1*d2
           return { num * other.den + other.num * den, den * other.den };
         }
         // subtract
         Fraction operator-(const Fraction& other) const
         {
           if (den == other.den)
             return Fraction(num - other.num, den);
           // n1/d1 - n2/d2 = (n1*d2 - n2*d1) / d1*d2
           return { num * other.den - other.num * den, den * other.den };
         }
         // multiply
         Fraction operator*(const Fraction& other) const
         {
           // n1/d1 * n2/d2 = n1*n2 / d1*d2
           return Fraction(num * other.num, den * other.den);
         }
         // divide
         Fraction operator/(const Fraction& other) const
         {
           // n1/d1 / n2/d2 = n1*d2 / d1*n2
           return Fraction(num * other.den, den * other.num);
         }
         // sort
         bool operator< (const Fraction& other) const
         {
           // n1/d1 < n2/d2 => n1*d2 < n2*d2
           return num * other.den <  other.num * den;
         }
         // compare
         bool operator==(const Fraction& other) const
         {
           // n1/d1 < n2/d2 => n1*d2 < n2*d2
           return num * other.den == other.num * den;
         }
        
         // return Fraction with swapped numerator and denominator
         Fraction inverse() const
         {
           return { den, num };
         }
        
         // convert to a proper fraction
         Fraction reduced() const
         {
           // catch simple cases: 0/d and n/1
           if (num == 0 || den == 1)
             return { num, 1 };
        
           // divide numerator and denominator by their gcd()
           // code taken from my toolbox / gcd()
           auto a = num;
           auto b = den;
           while (a != 0)
           {
             auto c = a;
             a = b % a;
             b = c;
           }
           // now b contains the gcd()
           if (b == 1)
             return *this;
           return { num/b, den/b };
         }
       };

Heatmap

Please click on a problem's number to open my solution to that problem:

green		solutions solve the original Project Euler problem and have a perfect score of 100% at Hackerrank, too
yellow		solutions score less than 100% at Hackerrank (but still solve the original problem easily)
gray		problems are already solved but I haven't published my solution yet
blue		solutions are relevant for Project Euler only: there wasn't a Hackerrank version of it (at the time I solved it) or it differed too much
orange		problems are solved but exceed the time limit of one minute or the memory limit of 256 MByte
red		problems are not solved yet but I wrote a simulation to approximate the result or verified at least the given example - usually I sketched a few ideas, too
black		problems are solved but access to the solution is blocked for a few days until the next problem is published
[new]		the flashing problem is the one I solved most recently

I stopped working on Project Euler problems around the time they released 617.

1	2	3	4	5	6	7	8	9	10	11	12	13	14	15	16	17	18	19	20	21	22	23	24	25
26	27	28	29	30	31	32	33	34	35	36	37	38	39	40	41	42	43	44	45	46	47	48	49	50
51	52	53	54	55	56	57	58	59	60	61	62	63	64	65	66	67	68	69	70	71	72	73	74	75
76	77	78	79	80	81	82	83	84	85	86	87	88	89	90	91	92	93	94	95	96	97	98	99	100

101	102	103	104	105	106	107	108	109	110	111	112	113	114	115	116	117	118	119	120	121	122	123	124	125
126	127	128	129	130	131	132	133	134	135	136	137	138	139	140	141	142	143	144	145	146	147	148	149	150
151	152	153	154	155	156	157	158	159	160	161	162	163	164	165	166	167	168	169	170	171	172	173	174	175
176	177	178	179	180	181	182	183	184	185	186	187	188	189	190	191	192	193	194	195	196	197	198	199	200

201	202	203	204	205	206	207	208	209	210	211	212	213	214	215	216	217	218	219	220	221	222	223	224	225
226	227	228	229	230	231	232	233	234	235	236	237	238	239	240	241	242	243	244	245	246	247	248	249	250
251	252	253	254	255	256	257	258	259	260	261	262	263	264	265	266	267	268	269	270	271	272	273	274	275
276	277	278	279	280	281	282	283	284	285	286	287	288	289	290	291	292	293	294	295	296	297	298	299	300

301	302	303	304	305	306	307	308	309	310	311	312	313	314	315	316	317	318	319	320	321	322	323	324	325
326	327	328	329	330	331	332	333	334	335	336	337	338	339	340	341	342	343	344	345	346	347	348	349	350
351	352	353	354	355	356	357	358	359	360	361	362	363	364	365	366	367	368	369	370	371	372	373	374	375
376	377	378	379	380	381	382	383	384	385	386	387	388	389	390	391	392	393	394	395	396	397	398	399	400

401	402	403	404	405	406	407	408	409	410	411	412	413	414	415	416	417	418	419	420	421	422	423	424	425
426	427	428	429	430	431	432	433	434	435	436	437	438	439	440	441	442	443	444	445	446	447	448	449	450
451	452	453	454	455	456	457	458	459	460	461	462	463	464	465	466	467	468	469	470	471	472	473	474	475
476	477	478	479	480	481	482	483	484	485	486	487	488	489	490	491	492	493	494	495	496	497	498	499	500

501	502	503	504	505	506	507	508	509	510	511	512	513	514	515	516	517	518	519	520	521	522	523	524	525
526	527	528	529	530	531	532	533	534	535	536	537	538	539	540	541	542	543	544	545	546	547	548	549	550
551	552	553	554	555	556	557	558	559	560	561	562	563	564	565	566	567	568	569	570	571	572	573	574	575
576	577	578	579	580	581	582	583	584	585	586	587	588	589	590	591	592	593	594	595	596	597	598	599	600

601	602	603	604	605	606	607	608	609	610	611	612	613	614	615	616	617	618	619	620	621	622	623	624	625
626	627	628	629	630	631	632	633	634	635	636	637	638	639	640	641	642	643	644	645	646	647	648	649	650
651	652	653	654	655	656	657	658	659	660	661	662	663	664	665	666	667	668	669	670	671	672	673	674	675
676	677	678	679	680	681	682	683	684	685	686	687	688	689	690	691	692	693	694	695	696	697	698	699	700

701	702	703	704	705	706	707	708	709	710	711	712	713	714	715	716	717	718	719	720	721	722	723	724	725
726	727	728	729	730	731	732	733	734	735	736	737	738	739	740	741	742	743	744	745	746	747	748	749	750
751	752	753	754	755	756	757	758	759	760	761	762	763	764	765	766	767	768	769	770	771	772	773	774	775
776	777	778	779	780	781	782	783	784	785	786	787	788	789	790	791	792	793	794	795	796	797	798	799	800

801

802

803

804

805

806

807

808

809

810

811

812

813

814

The 310 solved problems (that's level 12) had an average difficulty of 32.6% at Project Euler and
I scored 13526 points (out of 15700 possible points, top rank was 17 out of ≈60000 in August 2017) at Hackerrank's Project Euler+.

My username at Project Euler is stephanbrumme while it's stbrumme at Hackerrank.

Look at my progress and performance pages to get more details.

1	2	3	4	5	6	7	8	9	10	11	12	13	14	15	16	17	18	19	20	21	22	23	24	25
26	27	28	29	30	31	32	33	34	35	36	37	38	39	40	41	42	43	44	45	46	47	48	49	50
51	52	53	54	55	56	57	58	59	60	61	62	63	64	65	66	67	68	69	70	71	72	73	74	75
76	77	78	79	80	81	82	83	84	85	86	87	88	89	90	91	92	93	94	95	96	97	98	99	100

101	102	103	104	105	106	107	108	109	110	111	112	113	114	115	116	117	118	119	120	121	122	123	124	125
126	127	128	129	130	131	132	133	134	135	136	137	138	139	140	141	142	143	144	145	146	147	148	149	150
151	152	153	154	155	156	157	158	159	160	161	162	163	164	165	166	167	168	169	170	171	172	173	174	175
176	177	178	179	180	181	182	183	184	185	186	187	188	189	190	191	192	193	194	195	196	197	198	199	200

201	202	203	204	205	206	207	208	209	210	211	212	213	214	215	216	217	218	219	220	221	222	223	224	225
226	227	228	229	230	231	232	233	234	235	236	237	238	239	240	241	242	243	244	245	246	247	248	249	250
251	252	253	254	255	256	257	258	259	260	261	262	263	264	265	266	267	268	269	270	271	272	273	274	275
276	277	278	279	280	281	282	283	284	285	286	287	288	289	290	291	292	293	294	295	296	297	298	299	300

301	302	303	304	305	306	307	308	309	310	311	312	313	314	315	316	317	318	319	320	321	322	323	324	325
326	327	328	329	330	331	332	333	334	335	336	337	338	339	340	341	342	343	344	345	346	347	348	349	350
351	352	353	354	355	356	357	358	359	360	361	362	363	364	365	366	367	368	369	370	371	372	373	374	375
376	377	378	379	380	381	382	383	384	385	386	387	388	389	390	391	392	393	394	395	396	397	398	399	400

401	402	403	404	405	406	407	408	409	410	411	412	413	414	415	416	417	418	419	420	421	422	423	424	425
426	427	428	429	430	431	432	433	434	435	436	437	438	439	440	441	442	443	444	445	446	447	448	449	450
451	452	453	454	455	456	457	458	459	460	461	462	463	464	465	466	467	468	469	470	471	472	473	474	475
476	477	478	479	480	481	482	483	484	485	486	487	488	489	490	491	492	493	494	495	496	497	498	499	500

501	502	503	504	505	506	507	508	509	510	511	512	513	514	515	516	517	518	519	520	521	522	523	524	525
526	527	528	529	530	531	532	533	534	535	536	537	538	539	540	541	542	543	544	545	546	547	548	549	550
551	552	553	554	555	556	557	558	559	560	561	562	563	564	565	566	567	568	569	570	571	572	573	574	575
576	577	578	579	580	581	582	583	584	585	586	587	588	589	590	591	592	593	594	595	596	597	598	599	600

601	602	603	604	605	606	607	608	609	610	611	612	613	614	615	616	617	618	619	620	621	622	623	624	625
626	627	628	629	630	631	632	633	634	635	636	637	638	639	640	641	642	643	644	645	646	647	648	649	650
651	652	653	654	655	656	657	658	659	660	661	662	663	664	665	666	667	668	669	670	671	672	673	674	675
676	677	678	679	680	681	682	683	684	685	686	687	688	689	690	691	692	693	694	695	696	697	698	699	700

701	702	703	704	705	706	707	708	709	710	711	712	713	714	715	716	717	718	719	720	721	722	723	724	725
726	727	728	729	730	731	732	733	734	735	736	737	738	739	740	741	742	743	744	745	746	747	748	749	750
751	752	753	754	755	756	757	758	759	760	761	762	763	764	765	766	767	768	769	770	771	772	773	774	775
776	777	778	779	780	781	782	783	784	785	786	787	788	789	790	791	792	793	794	795	796	797	798	799	800

1	2	3	4	5	6	7	8	9	10	11	12	13	14	15	16	17	18	19	20	21	22	23	24	25
26	27	28	29	30	31	32	33	34	35	36	37	38	39	40	41	42	43	44	45	46	47	48	49	50
51	52	53	54	55	56	57	58	59	60	61	62	63	64	65	66	67	68	69	70	71	72	73	74	75
76	77	78	79	80	81	82	83	84	85	86	87	88	89	90	91	92	93	94	95	96	97	98	99	100

101	102	103	104	105	106	107	108	109	110	111	112	113	114	115	116	117	118	119	120	121	122	123	124	125
126	127	128	129	130	131	132	133	134	135	136	137	138	139	140	141	142	143	144	145	146	147	148	149	150
151	152	153	154	155	156	157	158	159	160	161	162	163	164	165	166	167	168	169	170	171	172	173	174	175
176	177	178	179	180	181	182	183	184	185	186	187	188	189	190	191	192	193	194	195	196	197	198	199	200

201	202	203	204	205	206	207	208	209	210	211	212	213	214	215	216	217	218	219	220	221	222	223	224	225
226	227	228	229	230	231	232	233	234	235	236	237	238	239	240	241	242	243	244	245	246	247	248	249	250
251	252	253	254	255	256	257	258	259	260	261	262	263	264	265	266	267	268	269	270	271	272	273	274	275
276	277	278	279	280	281	282	283	284	285	286	287	288	289	290	291	292	293	294	295	296	297	298	299	300

301	302	303	304	305	306	307	308	309	310	311	312	313	314	315	316	317	318	319	320	321	322	323	324	325
326	327	328	329	330	331	332	333	334	335	336	337	338	339	340	341	342	343	344	345	346	347	348	349	350
351	352	353	354	355	356	357	358	359	360	361	362	363	364	365	366	367	368	369	370	371	372	373	374	375
376	377	378	379	380	381	382	383	384	385	386	387	388	389	390	391	392	393	394	395	396	397	398	399	400

401	402	403	404	405	406	407	408	409	410	411	412	413	414	415	416	417	418	419	420	421	422	423	424	425
426	427	428	429	430	431	432	433	434	435	436	437	438	439	440	441	442	443	444	445	446	447	448	449	450
451	452	453	454	455	456	457	458	459	460	461	462	463	464	465	466	467	468	469	470	471	472	473	474	475
476	477	478	479	480	481	482	483	484	485	486	487	488	489	490	491	492	493	494	495	496	497	498	499	500

501	502	503	504	505	506	507	508	509	510	511	512	513	514	515	516	517	518	519	520	521	522	523	524	525
526	527	528	529	530	531	532	533	534	535	536	537	538	539	540	541	542	543	544	545	546	547	548	549	550
551	552	553	554	555	556	557	558	559	560	561	562	563	564	565	566	567	568	569	570	571	572	573	574	575
576	577	578	579	580	581	582	583	584	585	586	587	588	589	590	591	592	593	594	595	596	597	598	599	600

601	602	603	604	605	606	607	608	609	610	611	612	613	614	615	616	617	618	619	620	621	622	623	624	625
626	627	628	629	630	631	632	633	634	635	636	637	638	639	640	641	642	643	644	645	646	647	648	649	650
651	652	653	654	655	656	657	658	659	660	661	662	663	664	665	666	667	668	669	670	671	672	673	674	675
676	677	678	679	680	681	682	683	684	685	686	687	688	689	690	691	692	693	694	695	696	697	698	699	700

701	702	703	704	705	706	707	708	709	710	711	712	713	714	715	716	717	718	719	720	721	722	723	724	725
726	727	728	729	730	731	732	733	734	735	736	737	738	739	740	741	742	743	744	745	746	747	748	749	750
751	752	753	754	755	756	757	758	759	760	761	762	763	764	765	766	767	768	769	770	771	772	773	774	775
776	777	778	779	780	781	782	783	784	785	786	787	788	789	790	791	792	793	794	795	796	797	798	799	800

1	2	3	4	5	6	7	8	9	10	11	12	13	14	15	16	17	18	19	20	21	22	23	24	25
26	27	28	29	30	31	32	33	34	35	36	37	38	39	40	41	42	43	44	45	46	47	48	49	50
51	52	53	54	55	56	57	58	59	60	61	62	63	64	65	66	67	68	69	70	71	72	73	74	75
76	77	78	79	80	81	82	83	84	85	86	87	88	89	90	91	92	93	94	95	96	97	98	99	100

101	102	103	104	105	106	107	108	109	110	111	112	113	114	115	116	117	118	119	120	121	122	123	124	125
126	127	128	129	130	131	132	133	134	135	136	137	138	139	140	141	142	143	144	145	146	147	148	149	150
151	152	153	154	155	156	157	158	159	160	161	162	163	164	165	166	167	168	169	170	171	172	173	174	175
176	177	178	179	180	181	182	183	184	185	186	187	188	189	190	191	192	193	194	195	196	197	198	199	200

201	202	203	204	205	206	207	208	209	210	211	212	213	214	215	216	217	218	219	220	221	222	223	224	225
226	227	228	229	230	231	232	233	234	235	236	237	238	239	240	241	242	243	244	245	246	247	248	249	250
251	252	253	254	255	256	257	258	259	260	261	262	263	264	265	266	267	268	269	270	271	272	273	274	275
276	277	278	279	280	281	282	283	284	285	286	287	288	289	290	291	292	293	294	295	296	297	298	299	300

301	302	303	304	305	306	307	308	309	310	311	312	313	314	315	316	317	318	319	320	321	322	323	324	325
326	327	328	329	330	331	332	333	334	335	336	337	338	339	340	341	342	343	344	345	346	347	348	349	350
351	352	353	354	355	356	357	358	359	360	361	362	363	364	365	366	367	368	369	370	371	372	373	374	375
376	377	378	379	380	381	382	383	384	385	386	387	388	389	390	391	392	393	394	395	396	397	398	399	400

401	402	403	404	405	406	407	408	409	410	411	412	413	414	415	416	417	418	419	420	421	422	423	424	425
426	427	428	429	430	431	432	433	434	435	436	437	438	439	440	441	442	443	444	445	446	447	448	449	450
451	452	453	454	455	456	457	458	459	460	461	462	463	464	465	466	467	468	469	470	471	472	473	474	475
476	477	478	479	480	481	482	483	484	485	486	487	488	489	490	491	492	493	494	495	496	497	498	499	500

501	502	503	504	505	506	507	508	509	510	511	512	513	514	515	516	517	518	519	520	521	522	523	524	525
526	527	528	529	530	531	532	533	534	535	536	537	538	539	540	541	542	543	544	545	546	547	548	549	550
551	552	553	554	555	556	557	558	559	560	561	562	563	564	565	566	567	568	569	570	571	572	573	574	575
576	577	578	579	580	581	582	583	584	585	586	587	588	589	590	591	592	593	594	595	596	597	598	599	600

601	602	603	604	605	606	607	608	609	610	611	612	613	614	615	616	617	618	619	620	621	622	623	624	625
626	627	628	629	630	631	632	633	634	635	636	637	638	639	640	641	642	643	644	645	646	647	648	649	650
651	652	653	654	655	656	657	658	659	660	661	662	663	664	665	666	667	668	669	670	671	672	673	674	675
676	677	678	679	680	681	682	683	684	685	686	687	688	689	690	691	692	693	694	695	696	697	698	699	700

701	702	703	704	705	706	707	708	709	710	711	712	713	714	715	716	717	718	719	720	721	722	723	724	725
726	727	728	729	730	731	732	733	734	735	736	737	738	739	740	741	742	743	744	745	746	747	748	749	750
751	752	753	754	755	756	757	758	759	760	761	762	763	764	765	766	767	768	769	770	771	772	773	774	775
776	777	778	779	780	781	782	783	784	785	786	787	788	789	790	791	792	793	794	795	796	797	798	799	800