Sum of all digits from 1 to N — FAST

A job-interview problem I liked too much to leave at the brute-force answer. Ended up at roughly O(log N) — an N with 190 digits takes about 75 ms, an N with 1300 digits takes under 7 seconds. The trick is three recursive relations that peel the number apart and leave only a handful of tiny calls to brute-sum at the bottom.

role

Solo · just-for-the-fun-of-it problem

kind

Job-interview task · algorithm play

language

Python 2

best run

N with 190 digits in ≈ 75 ms · N with 1300 digits in ≈ 7 s

The problem. Find the sum of the digits of every integer from 1 to N, for both ends included. For N = 10 that's 1 + 2 + 3 + 4 + 5 + 6 + 7 + 8 + 9 + (1 + 0) = 46. Now do it for N with a thousand digits.

The task showed up in a job interview. I liked it enough that I wanted to see how far it could be pushed.

Brute force

The honest baseline.

def sum_of_digits(n):
    s = 0
    while n > 0:
        s += n % 10
        n /= 10
    return s

def brute_sum(N):
    n = N
    ds = 0
    while n > 0:
        ds += sum_of_digits(n)
        n -= 1
    return ds

Nothing fancy. For a seven-digit N it spends just over 14 seconds; anything past that and you can start making coffee.

Three relations that do the actual work

The optimisation is three recursive relations. Each peels the problem down one level; together they leave brute-sum called only on numbers 0 through 10.

Relation 1 — split a messy number into round pieces:

foo(3456) == foo(3000)
           + foo(400) + 400 * (3)
           + foo(50)  + 50  * (3 + 4)
           + foo(6)   + 6   * (3 + 4 + 5)

This reduces an arbitrary N to a set of calls in the shape L × 10^M — a single leading digit followed by zeros.

Relation 2 — collapse L × 10^M to 1 × 10^M via triangular numbers:

triangular = [0, 1, 3, 6, 10, 15, 21, 28, 36]  # 0 unused
foo(3000) == 3 * foo(1000) + triangular[3 - 1] * 1000

Those triangular numbers fell out of the arithmetic on paper — I did not recognise them at first and was cheerfully calling them "empirical thingies" until I did.

Relation 3 — the one that's truly recursive: reduce 1 × 10^M down one order of magnitude:

foo(1000) == foo(100) * 10 + 44 * 100 + 100 - 9

The 44 and the 9 I don't fully have an intuition for. The tenth triangular number is 45, which is suspicious; I didn't dissect it further.

Result

With the three relations in place, actual brute-force summing only runs for numbers up to 10 — and only a handful of times per call. The main solver:

def round10_sum(lead_digit, power10):
    rounds = 10 ** power10
    n = lead_digit * rounds
    if n > 10:
        if lead_digit == 1:
            lower = 10 ** (power10 - 1)
            return round10_sum(1, power10 - 1) * 10 + 44 * lower + lower - 9
        return lead_digit * round10_sum(1, power10) + triangular[lead_digit - 1] * rounds
    return cached_sum(n)

def optimised_digits_sum(n):
    digs = map(int, str(n))
    total = 0
    power10 = 0
    while len(digs):
        lead = digs.pop(-1)
        total += round10_sum(lead, power10) + lead * sum(digs) * 10 ** power10
        power10 += 1
    return total

That's the whole thing. Runs at roughly O(log N). An N with 190 digits completes in about 75 ms; an N with 1300 digits in under 7 seconds. For a range N to M rather than 1 to N, compute both and subtract.

Full code

For completeness — the complete Python 2 program, driven off a big N (the first ~400 digits of pi, for fun):

import time
import sys
sys.setrecursionlimit(10000)

numberpie = 31415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679821480865132823066470938446095505822317253594081284811174502841027019385211055596446229489549303819644288109756659334461284756482337867831652712019091456485669234603486104543266482133936072602491412737245870066063155881748815209209628292540917153643678925903600113305305488204665213841469519415116094330572703657595919530921861173819326117931051185480744623799627495673518857527248912279381830119491298336733624406566430860213949463952247371907021798609437027705392171762931767523846748184676694051320005681271452635608277857713427577896091736371787214684409012249534301465495853710507922796892589235420199561121290219608640344181598136297747713099605187072113499999983729780499510597317328160963185950244594553469083026425223082533446850352619311881710100031378387528865875332083814206171776691473035982534904287554687311595628638823537875937519577818577805321712268066130019278766111959092164201989380952572010654858632788659361533818279682303019520353018529689957736225994138912497217752834791315155748572424541506959508295331168617278558890750983817546374649393192550604009277016711390098488240128583616035637076601047101819429555961989467678374494482553797747268471040475346462080466842590694912933136770289891521047521620569660241  # rounded :p 405

N = numberpie

# --- brute force ------------------------------------------------

def sum_of_digits(n):
    """Sum all digits of n."""
    s = 0
    while n > 0:
        s += n % 10
        n /= 10
    return s

def brute_sum(N):
    """Sum of digits of all numbers 1..N."""
    n = N
    ds = 0
    while n > 0:
        ds += sum_of_digits(n)
        n -= 1
    return ds

cache = {}

def cached_sum(n):
    """Trivial cache around brute_sum()."""
    ds = cache.get(n, None)
    if ds is None:
        ds = brute_sum(n)
        cache[n] = ds
    return ds

# --- optimised --------------------------------------------------

triangular_numbers = [0, 1, 3, 6, 10, 15, 21, 28, 36]  # index 0 unused
meaning_of_life = 44  # :p @ 42
the_number_ni = 9

def round10_sum(lead_digit, power10):
    """Round numbers (single digit followed by zeroes).
    Calls cached_sum() only for numbers 0..10, otherwise recursive.
    """
    rounds = 10 ** power10
    n = lead_digit * rounds
    if n > 10:
        if lead_digit == 1:
            lower_power10 = 10 ** (power10 - 1)
            digit_sum = round10_sum(1, power10 - 1) * 10 \
                      + meaning_of_life * lower_power10 \
                      + lower_power10 - the_number_ni
        else:
            digit_sum = lead_digit * round10_sum(1, power10) \
                      + triangular_numbers[lead_digit - 1] * rounds
    else:
        digit_sum = cached_sum(n)
    return digit_sum

def optimised_digits_sum(n):
    """Reduce any N to round10_sum() calls on round numbers."""
    digs = map(int, str(n))
    digit_sum = 0
    power10 = 0
    while len(digs):
        lead_digit = digs.pop(-1)
        digit_sum += round10_sum(lead_digit, power10) \
                   + lead_digit * sum(digs) * 10 ** power10
        power10 += 1
    return digit_sum

# --- go ---------------------------------------------------------

start_time = time.time()
print optimised_digits_sum(N)
end_time = time.time()
print "Time: %g s" % (end_time - start_time)

# Examples of the relations used in the recursion:
#
# print optimised_digits_sum(1000)
# print optimised_digits_sum(100) * 10 + 44 * 100 + 100 - 9
#
# print optimised_digits_sum(6000)
# print 6 * optimised_digits_sum(1000) + triangular_numbers[5] * 1000
#
# print optimised_digits_sum(3456)
# print optimised_digits_sum(3000) + 3 * 0 \
#     + optimised_digits_sum(400) + 4 * (3)     * 100 \
#     + optimised_digits_sum(50)  + 5 * (3 + 4) * 10 \
#     + optimised_digits_sum(6)   + 6 * (3 + 4 + 5)