Problem 56 — Powerful Digit Sum

Problem Analysis

This problem asks for the maximum digit sum - the sum of the decimal digits of a number - over all values of the form a^b with a and b both below 100. The digit sum of 1234, for instance, is 1+2+3+4 = 10. There is no closed-form shortcut here: the digit sum is a property of the decimal representation, and the decimal representation of a power has no simple analytic structure. The problem is therefore inherently computational. Its one genuine obstacle is size: a^b for a and b near 99 has roughly b * log₁₀(a) ~= 99 * 2 ~= 196 decimal digits, far beyond what a 64-bit integer can hold. We need arbitrary-precision arithmetic (commonly called "bignum" arithmetic) to compute these numbers exactly, because a single rounding error anywhere in the high digits would corrupt the digit sum.

The observation that makes the problem tractable is that the search space is minuscule. There are at most 99 * 99 = 9801 pairs (a, b), and even after generalising the bound to 10^num_digits the count stays modest. So a brute-force search suffices: enumerate every pair, compute a^b exactly, sum its digits, keep a running maximum. Letting N be the range limit and D ~= N * log₁₀(N) the digit count of the largest value, computing one power by repeated multiplication costs O(N * D) digit-operations and summing its digits costs O(D), so the whole search is roughly O(N³ log N) in the worst case. For the canonical N = 99 this is trivially fast.

Solution Approach

Brute-force big-integer digit-sum search (Solution 0)

Both the Python and C files implement the same algorithm: iterate over every (base, exp) pair in the search range, compute base**exp exactly, sum its digits, and track the maximum. The search is parameterised by num_digits, with stop_at = 10**num_digits as the exclusive upper bound, so num_digits=2 reproduces the original Project Euler range 1..99.

A deliberate, pragmatic choice in both files is that they do not scan the full range. They restrict the base to the top 100 values and the exponent to the top 10 values below stop_at (for num_digits=2: bases 1..99, exponents 90..99). The reasoning is empirical. A digit sum tends to be largest when a number has many large digits, and the most promising candidates are large bases raised to high exponents - small bases or small exponents produce shorter numbers with fewer digits to sum. This pruning yields the correct answer for every test case here, but it is a heuristic shortcut rather than a proven bound; a fully rigorous solution would scan all pairs. The Python --show mode exists partly to sanity-check exactly this kind of assumption by visualising where the maximum actually lands in the grid.

In the Python solution, the built-in int is already arbitrary-precision, so base**exp directly produces the exact value with no library and no overflow. I extract the digit sum with the idiom sum(int(d) for d in str(base**exp)): convert the number to its decimal string and sum the integer value of each character. The entire search collapses into a single comprehension passed to max(). The optional --show branch builds NumPy matrices of digit sums, digit lengths, and average digit values across the grid and renders three Matplotlib heatmaps - a diagnostic aid entirely separate from the answer computation.

In the C solution, there is no native bignum, so I implement a minimal one. A BigInt is a heap-allocated struct holding a little-endian array of single decimal digits (index 0 = units place), a length, and a capacity. I chose base-10, one-digit-per-element storage on purpose: it makes the digit sum a trivial loop over the array with no extraction at all - exactly the operation this problem performs constantly. The only multiplication I need is bigint_mul_int (multiply by a scalar), a left-to-right carry-propagation loop. To compute the powers efficiently, for each base I build base^exp_start once via repeated scalar multiplies, then for each later exponent I multiply by the base just once. This incremental scheme means every exponent step after the first costs exactly one scalar multiply instead of recomputing the power from scratch. I pre-allocate the digit buffer to a conservative upper bound (max_exp * (num_digits+1) + 100) and keep a doubling-realloc fallback in the multiply routine as a safety net.