Euler problem #16 has been discussed many times here, but I could not find an answer that gives a good overview of possible solution approaches, the lay of the land as it were. Here's my attempt at rectifying that.
This overview is intended for people who have already found a solution and want to get a more complete picture. It is basically language-agnostic even though the sample code is C#. There are some usages of features that are not available in C# 2.0 but they are not essential - their purpose is only to get boring stuff out of the way with a minimum of fuss.
Apart from using a ready-made BigInteger library (which doesn't count), straightforward solutions for Euler #16 fall into two fundamental categories: performing calculations natively - i.e. in a base that is a power of two - and converting to decimal in order to get at the digits, or performing the computations directly in a decimal base so that the digits are available without any conversion.
For the latter there are two reasonably simple options:
- repeated doubling
- powering by repeated squaring
Native Computation + Radix Conversion
This approach is the simplest and its performance exceeds that of naive solutions using .Net's builtin BigInteger
type.
The actual computation is trivially achieved: just perform the moral equivalent of 1 << 1000
, by storing 1000 binary zeroes and appending a single lone binary 1.
The conversion is also quite simple and can be done by coding the pencil-and-paper division method, with a suitably large choice of 'digit' for efficiency. Variables for intermediate results need to be able to hold two 'digits'; dividing the number of decimal digits that fit in a long
by 2 gives 9 decimal digits for the maximum meta-digit (or 'limb', as it is usually called in bignum lore).
class E16_RadixConversion
{
const int BITS_PER_WORD = sizeof(uint) * 8;
const uint RADIX = 1000000000; // == 10^9
public static int digit_sum_for_power_of_2 (int exponent)
{
var dec = new List<int>();
var bin = new uint[(exponent + BITS_PER_WORD) / BITS_PER_WORD];
int top = bin.Length - 1;
bin[top] = 1u << (exponent % BITS_PER_WORD);
while (top >= 0)
{
ulong rest = 0;
for (int i = top; i >= 0; --i)
{
ulong temp = (rest << BITS_PER_WORD) | bin[i];
ulong quot = temp / RADIX; // x64 uses MUL (sometimes), x86 calls a helper function
rest = temp - quot * RADIX;
bin[i] = (uint)quot;
}
dec.Add((int)rest);
if (bin[top] == 0)
--top;
}
return E16_Common.digit_sum(dec);
}
}
I wrote (rest << BITS_PER_WORD) | big[i]
instead of using operator + because that is precisely what is needed here; no 64-bit addition with carry propagation needs to take place. This means that the two operands could be written directly to their separate registers in a register pair, or to fields in an equivalent struct like LARGE_INTEGER
.
On 32-bit systems the 64-bit division cannot be inlined as a few CPU instructions, because the compiler cannot know that the algorithm guarantees quotient and remainder to fit into 32-bit registers. Hence the compiler calls a helper function that can handle all eventualities.
These systems may profit from using a smaller limb, i.e. RADIX = 10000
and uint
instead of ulong
for holding intermediate (double-limb) results. An alternative for languages like C/C++ would be to call a suitable compiler intrinsic that wraps the raw 32-bit by 32-bit to 64-bit multiply (assuming that division by the constant radix is to be implemented by multiplication with the inverse). Conversely, on 64-bit systems the limb size can be increased to 19 digits if the compiler offers a suitable 64-by-64-to-128 bit multiply primitive or allows inline assembler.
Decimal Doubling
Repeated doubling seems to be everyone's favourite, so let's do that next. Variables for intermediate results need to hold one 'digit' plus one carry bit, which gives 18 digits per limb for long
. Going to ulong
cannot improve things (there's 0.04 bit missing to 19 digits plus carry), and so we might as well stick with long
.
On a binary computer, decimal limbs do not coincide with computer word boundaries. That makes it necessary to perform a modulo operation on the limbs during each step of the calculation. Here, this modulo op can be reduced to a subtraction of the modulus in the event of carry, which is faster than performing a division. The branching in the inner loop can be eliminated by bit twiddling but that would be needlessly obscure for a demonstration of the basic algorithm.
class E16_DecimalDoubling
{
const int DIGITS_PER_LIMB = 18; // == floor(log10(2) * (63 - 1)), b/o carry
const long LIMB_MODULUS = 1000000000000000000L; // == 10^18
public static int digit_sum_for_power_of_2 (int power_of_2)
{
Trace.Assert(power_of_2 > 0);
int total_digits = (int)Math.Ceiling(Math.Log10(2) * power_of_2);
int total_limbs = (total_digits + DIGITS_PER_LIMB - 1) / DIGITS_PER_LIMB;
var a = new long[total_limbs];
int limbs = 1;
a[0] = 2;
for (int i = 1; i < power_of_2; ++i)
{
int carry = 0;
for (int j = 0; j < limbs; ++j)
{
long new_limb = (a[j] << 1) | carry;
carry = 0;
if (new_limb >= LIMB_MODULUS)
{
new_limb -= LIMB_MODULUS;
carry = 1;
}
a[j] = new_limb;
}
if (carry != 0)
{
a[limbs++] = carry;
}
}
return E16_Common.digit_sum(a);
}
}
This is just as simple as radix conversion, but except for very small exponents it does not perform anywhere near as well (despite its huge meta-digits of 18 decimal places). The reason is that the code must perform (exponent - 1) doublings, and the work done in each pass corresponds to about half the total number of digits (limbs).
Repeated Squaring
The idea behind powering by repeated squaring is to replace a large number of doublings with a small number of multiplications.
1000 = 2^3 + 2^5 + 2^6 + 2^7 + 2^8 + 2^9
x^1000 = x^(2^3 + 2^5 + 2^6 + 2^7 + 2^8 + 2^9)
x^1000 = x^2^3 * x^2^5 * x^2^6 * x^2^7 * x^2*8 * x^2^9
x^2^3 can be obtained by squaring x three times, x^2^5 by squaring five times, and so on. On a binary computer the decomposition of the exponent into powers of two is readily available because it is the bit pattern representing that number. However, even non-binary computers should be able to test whether a number is odd or even, or to divide a number by two.
The multiplication can be done by coding the pencil-and-paper method; here I'm using a helper function that computes one row of a product and adds it into the result at a suitably shifted position, so that the rows of partial products do not need to be stored for a separate addition step later. Intermediate values during computation can be up to two 'digits' in size, so that the limbs can be only half as wide as for repeated doubling (where only one extra bit had to fit in addition to a 'digit').
Note: the radix of the computations is not a power of 2, and so the squarings of 2 cannot be computed by simple shifting here. On the positive side, the code can be used for computing powers of bases other than 2.
class E16_DecimalSquaring
{
const int DIGITS_PER_LIMB = 9; // language limit 18, half needed for holding the carry
const int LIMB_MODULUS = 1000000000;
public static int digit_sum_for_power_of_2 (int e)
{
Trace.Assert(e > 0);
int total_digits = (int)Math.Ceiling(Math.Log10(2) * e);
int total_limbs = (total_digits + DIGITS_PER_LIMB - 1) / DIGITS_PER_LIMB;
var squared_power = new List<int>(total_limbs) { 2 };
var result = new List<int>(total_limbs);
result.Add((e & 1) == 0 ? 1 : 2);
while ((e >>= 1) != 0)
{
squared_power = multiply(squared_power, squared_power);
if ((e & 1) == 1)
result = multiply(result, squared_power);
}
return E16_Common.digit_sum(result);
}
static List<int> multiply (List<int> lhs, List<int> rhs)
{
var result = new List<int>(lhs.Count + rhs.Count);
resize_to_capacity(result);
for (int i = 0; i < rhs.Count; ++i)
addmul_1(result, i, lhs, rhs[i]);
trim_leading_zero_limbs(result);
return result;
}
static void addmul_1 (List<int> result, int offset, List<int> multiplicand, int multiplier)
{
// it is assumed that the caller has sized `result` appropriately before calling this primitive
Trace.Assert(result.Count >= offset + multiplicand.Count + 1);
long carry = 0;
foreach (long limb in multiplicand)
{
long temp = result[offset] + limb * multiplier + carry;
carry = temp / LIMB_MODULUS;
result[offset++] = (int)(temp - carry * LIMB_MODULUS);
}
while (carry != 0)
{
long final_temp = result[offset] + carry;
carry = final_temp / LIMB_MODULUS;
result[offset++] = (int)(final_temp - carry * LIMB_MODULUS);
}
}
static void resize_to_capacity (List<int> operand)
{
operand.AddRange(Enumerable.Repeat(0, operand.Capacity - operand.Count));
}
static void trim_leading_zero_limbs (List<int> operand)
{
int i = operand.Count;
while (i > 1 && operand[i - 1] == 0)
--i;
operand.RemoveRange(i, operand.Count - i);
}
}
The efficiency of this approach is roughly on par with radix conversion but there are specific improvements that apply here. Efficiency of the squaring can be doubled by writing a special squaring routine that utilises the fact that ai*bj == aj*bi
if a == b
, which cuts the number of multiplications in half.
Also, there are methods for computing addition chains that involve fewer operations overall than using the exponent bits for determining the squaring/multiplication schedule.
Helper Code and Benchmarks
The helper code for summing decimal digits in the meta-digits (decimal limbs) produced by the sample code is trivial, but I'm posting it here anyway for your convenience:
internal class E16_Common
{
internal static int digit_sum (int limb)
{
int sum = 0;
for ( ; limb > 0; limb /= 10)
sum += limb % 10;
return sum;
}
internal static int digit_sum (long limb)
{
const int M1E9 = 1000000000;
return digit_sum((int)(limb / M1E9)) + digit_sum((int)(limb % M1E9));
}
internal static int digit_sum (IEnumerable<int> limbs)
{
return limbs.Aggregate(0, (sum, limb) => sum + digit_sum(limb));
}
internal static int digit_sum (IEnumerable<long> limbs)
{
return limbs.Select((limb) => digit_sum(limb)).Sum();
}
}
This can be made more efficient in various ways but overall it is not critical.
All three solutions take O(n^2) time where n is the exponent. In other words, they will take a hundred times as long when the exponent grows by a factor of ten. Radix conversion and repeated squaring can both be improved to roughly O(n log n) by employing divide-and-conquer strategies; I doubt whether the doubling scheme can be improved in a similar fastion but then it was never competitive to begin with.
All three solutions presented here can be used to print the actual results, by stringifying the meta-digits with suitable padding and concatenating them. I've coded the functions as returning the digit sum instead of the arrays/lists with decimal limbs only in order to keep the sample code simple and to ensure that all functions have the same signature, for benchmarking.
In these benchmarks, the .Net BigInteger type was wrapped like this:
static int digit_sum_via_BigInteger (int power_of_2)
{
return System.Numerics.BigInteger.Pow(2, power_of_2)
.ToString()
.ToCharArray()
.Select((c) => (int)c - '0')
.Sum();
}
Finally, the benchmarks for the C# code:
# testing decimal doubling ...
1000: 1366 in 0,052 ms
10000: 13561 in 3,485 ms
100000: 135178 in 339,530 ms
1000000: 1351546 in 33.505,348 ms
# testing decimal squaring ...
1000: 1366 in 0,023 ms
10000: 13561 in 0,299 ms
100000: 135178 in 24,610 ms
1000000: 1351546 in 2.612,480 ms
# testing radix conversion ...
1000: 1366 in 0,018 ms
10000: 13561 in 0,619 ms
100000: 135178 in 60,618 ms
1000000: 1351546 in 5.944,242 ms
# testing BigInteger + LINQ ...
1000: 1366 in 0,021 ms
10000: 13561 in 0,737 ms
100000: 135178 in 69,331 ms
1000000: 1351546 in 6.723,880 ms
As you can see, the radix conversion is almost as slow as the solution using the builtin BigInteger class. The reason is that the runtime is of the newer type that does performs certain standard optimisations only for signed integer types but not for unsigned ones (here: implementing division by a constant as multiplication with the inverse).
I haven't found an easy means of inspecting the native code for existing .Net assemblies, so I decided on a different path of investigation: I coded a variant of E16_RadixConversion
for comparison where ulong
and uint
were replaced by long
and int
respectively, and BITS_PER_WORD
decreased by 1 accordingly. Here are the timings:
# testing radix conv Int63 ...
1000: 1366 in 0,004 ms
10000: 13561 in 0,202 ms
100000: 135178 in 18,414 ms
1000000: 1351546 in 1.834,305 ms
More than three times as fast as the version that uses unsigned types! Clear evidence of numbskullery in the compiler...
In order to showcase the effect of different limb sizes I templated the solutions in C++ on the unsigned integer types used as limbs. The timings are prefixed with the byte size of a limb and the number of decimal digits in a limb, separated by a colon. There is no timing for the often-seen case of manipulating digit characters in strings, but it is safe to say that such code will take at least twice as long as the code that uses double digits in byte-sized limbs.
# E16_DecimalDoubling
[1:02] e = 1000 -> 1366 0.308 ms
[2:04] e = 1000 -> 1366 0.152 ms
[4:09] e = 1000 -> 1366 0.070 ms
[8:18] e = 1000 -> 1366 0.071 ms
[1:02] e = 10000 -> 13561 30.533 ms
[2:04] e = 10000 -> 13561 13.791 ms
[4:09] e = 10000 -> 13561 6.436 ms
[8:18] e = 10000 -> 13561 2.996 ms
[1:02] e = 100000 -> 135178 2719.600 ms
[2:04] e = 100000 -> 135178 1340.050 ms
[4:09] e = 100000 -> 135178 588.878 ms
[8:18] e = 100000 -> 135178 290.721 ms
[8:18] e = 1000000 -> 1351546 28823.330 ms
For the exponent of 10^6 there is only the timing with 64-bit limbs, since I didn't have the patience to wait many minutes for full results. The picture is similar for radix conversion, except that there is no row for 64-bit limbs because my compiler does not have a native 128-bit integral type.
# E16_RadixConversion
[1:02] e = 1000 -> 1366 0.080 ms
[2:04] e = 1000 -> 1366 0.026 ms
[4:09] e = 1000 -> 1366 0.048 ms
[1:02] e = 10000 -> 13561 4.537 ms
[2:04] e = 10000 -> 13561 0.746 ms
[4:09] e = 10000 -> 13561 0.243 ms
[1:02] e = 100000 -> 135178 445.092 ms
[2:04] e = 100000 -> 135178 68.600 ms
[4:09] e = 100000 -> 135178 19.344 ms
[4:09] e = 1000000 -> 1351546 1925.564 ms
The interesting thing is that simply compiling the code as C++ doesn't make it any faster - i.e., the optimiser couldn't find any low-hanging fruit that the C# jitter missed, apart from not toeing the line with regard to penalising unsigned integers. That's the reason why I like prototyping in C# - performance in the same ballpark as (unoptimised) C++ and none of the hassle.
Here's the meat of the C++ version (sans reams of boring stuff like helper templates and so on) so that you can see that I didn't cheat to make C# look better:
template<typename W>
struct E16_RadixConversion
{
typedef W limb_t;
typedef typename detail::E16_traits<W>::long_t long_t;
static unsigned const BITS_PER_WORD = sizeof(limb_t) * CHAR_BIT;
static unsigned const RADIX_DIGITS = std::numeric_limits<limb_t>::digits10;
static limb_t const RADIX = detail::pow10_t<limb_t, RADIX_DIGITS>::RESULT;
static unsigned digit_sum_for_power_of_2 (unsigned e)
{
std::vector<limb_t> digits;
compute_digits_for_power_of_2(e, digits);
return digit_sum(digits);
}
static void compute_digits_for_power_of_2 (unsigned e, std::vector<limb_t> &result)
{
assert(e > 0);
unsigned total_digits = unsigned(std::ceil(std::log10(2) * e));
unsigned total_limbs = (total_digits + RADIX_DIGITS - 1) / RADIX_DIGITS;
result.resize(0);
result.reserve(total_limbs);
std::vector<limb_t> bin((e + BITS_PER_WORD) / BITS_PER_WORD);
bin.back() = limb_t(limb_t(1) << (e % BITS_PER_WORD));
while (!bin.empty())
{
long_t rest = 0;
for (std::size_t i = bin.size(); i-- > 0; )
{
long_t temp = (rest << BITS_PER_WORD) | bin[i];
long_t quot = temp / RADIX;
rest = temp - quot * RADIX;
bin[i] = limb_t(quot);
}
result.push_back(limb_t(rest));
if (bin.back() == 0)
bin.pop_back();
}
}
};
Conclusion
These benchmarks also show that this Euler task - like many others - seems designed to be solved on a ZX81 or an Apple ][, not on our modern toys that are a million times as powerful. There's no challenge involved here unless the limits are increased drastically (an exponent of 10^5 or 10^6 would be much more adequate).
A good overview of the practical state of the art can be got from GMP's overview of algorithms. Another excellent overview of the algorithms is chapter 1 of "Modern Computer Arithmetic" by Richard Brent and Paul Zimmermann. It contains exactly what one needs to know for coding challenges and competitions, but unfortunately the depth is not equal to that of Donald Knuth's treatment in "The Art of Computer Programming".
The radix conversion solution adds a useful technique to one's code challenge toolchest, since the given code can be trivially extended for converting any old big integer instead of only the bit pattern 1 << exponent
. The repeated squaring solutiono can be similarly useful since changing the sample code to power something other than 2 is again trivial.
The approach of performing computations directly in powers of 10 can be useful for challenges where decimal results are required, because performance is in the same ballpark as native computation but there is no need for a separate conversion step (which can require similar amounts of time as the actual computation).