Excellent Numbers

March 24, 2015

We begin by writing a function that recognizes excellent numbers:

(define (xl? n)
  (let* ((k (/ (+ (ilog 10 n) 1) 2))
         (ten (expt 10 k))
         (a (quotient n ten))
         (b (modulo n ten)))
    (= (- (* b b) (* a a)) n)))

Then we compute the sum of the six-digit excellent numbers:

> (time (sum (filter xl? (range 100000 1000000))))
(time (sum (filter xl? ...)))
    175 collections
    4908 ms elapsed cpu time, including 187 ms collecting
    4913 ms elapsed real time, including 177 ms collecting
    1468845552 bytes allocated, including 1463474816 bytes reclaimed

That’s slow because it repeatedly calculates the number of digits. Moving that computation out of the inner loop gives the following program:

(define (xl-sum k)
  (let ((ten (expt 10 k))
        (start (expt 10 (+ k k -1)))
        (stop (expt 10 (+ k k))))
    (let loop ((n start) (sum 0))
      (if (= n stop) sum
        (let ((a (quotient n ten))
              (b (modulo n ten)))
          (if (= (- (* b b) (* a a)) n)
              (loop (+ n 1) (+ sum n))
              (loop (+ n 1) sum)))))))

And now it’s much faster:

> (time (xl-sum 3))
(time (xl-sum 3))
    7 collections
    293 ms elapsed cpu time, including 7 ms collecting
    293 ms elapsed real time, including 7 ms collecting
    57602000 bytes allocated, including 58954160 bytes reclaimed

But still not fast enough. Here is the timing for eight-digit excellent numbers, which takes almost half a minute:

> (time (xl-sum 4))
(time (xl-sum 4))
    684 collections
    28921 ms elapsed cpu time, including 284 ms
    29133 ms elapsed real time, including 271 ms collectingcollecting
    5760175312 bytes allocated, including 5903355296 bytes reclaimed

At that rate, (xl-sum 5) will take almost an hour to process the nine billion ten-digit numbers. We need something better.

The 2k-digit numbers that we are processing are represented as a 10k + b. By definition, that is equal to b2a2 for excellent numbers. Rearranging terms to put the a on the left side and the b on the right side, we get a (10k + a) = b2b. There are 90,000 possible a among the ten-digit numbers; if we could determine quickly the corresponding b and test if it forms an excellent number, then we would have a very fast method to find the ten-digit excellent numbers.

Fortunately, it is easy to find b. Since b2b is very close to b2, just a little bit smaller, b − 1 ≤ sqrt(b2b) ≤ b and we can compute b = ⌈a (10k + a)⌉. Then all we need to do is test if a and b combine to form an excellent number:

(define (xl-list k)
  (let ((ten (expt 10 k))
        (start (expt 10 (- k 1)))
        (stop (expt 10 k)))
    (let loop ((a start) (xls (list)))
      (if (= a stop) (reverse xls)
        (let* ((b (+ (isqrt (* a (+ ten a))) 1))
               (n (+ (* a ten) b)))
          (if (= (- (* b b) (* a a)) n)
              (loop (+ a 1) (cons n xls))
              (loop (+ a 1) xls)))))))

Now it is very fast to calculate the ten-digit excellent numbers:

> (define xl-ten #f)
> (time (set! xl-ten (xl-list 5)))(time (set! xl-ten ...))
    17 collections
    441 ms elapsed cpu time, including 1 ms collecting
    442 ms elapsed real time, including 3 ms collecting
    145748000 bytes allocated, including 143206224 bytes reclaimed
> (sum xl-ten)
> xl-ten
(3333466668 4848484848 4989086476)

We used filter, range, sum, isqrt and ilog from the Standard Prelude. You can run the program at http://ideone.com/tFJBLg.

You might like to look at http://oeis.org/A162700, which gives a different definition of excellent numbers. You might also enjoy Mark Dominus’ version of the program, from which we stole both the exercise and the essence of the solution. A complete list of excellent numbers up to twenty digits appears on the next page.


Pages: 1 2 3

14 Responses to “Excellent Numbers”

  1. Rutger said

    Mm. b*b – a*a = n. Rewriting n to be a*100000 + b, yields b * b – b == a * (a + 100000). Which can be reduced to b = ((4*a**2+400000*a+1)**0.5+1) / 2. So we loop through all a values with 5 digits to find b and see if b is an integer. In Python:

    total = 0
    for a in range(10000, 100000):
    	b = ((4*a**2+400000*a+1)**0.5+1) / 2.0
    	if b == int(b):
    		result = int(str(a)+str(int(b)))
    		total += result
    print total

    Or in one line (rewriting b):

    print sum([int(str(a)+str(int(((4*a**2+400000*a+1)**0.5+1) / 2.0))) for a in range(10000, 100000) if (((4*a**2+400000*a+1)**0.5+1) / 2.0) == int(((4*a**2+400000*a+1)**0.5+1) / 2.0)])

    Both result in 13171037992

  2. Paul said

    In Python. Same as Rutger’s. The first part of the number a needs to be even. I think, there should be a check that b is not too large. For n=10 this never happens, but for n=12 it happens.

    import math
    def excellent(size):
        half = size // 2
        factor = 10 ** half
        for a in range(10**(half-1), factor, 2):
            b = (1 + math.sqrt(1 + 4 * (a ** 2 + factor * a))) / 2
            if b.is_integer() and b < factor:
                yield int(b) + a * factor
    print(sum(excellent(10))) #13171037992
  3. namako said

    Strictly speaking, you probably ought to check that b isn’t an integer because it’s missing low bits.
    That is, if b is an integer, make sure it’s actually a solution.

  4. matthew said

    We want solutions to a*(a+10000) = b*(b-1), both sides are monotonic, so we can just scan up through the two sequences, looking for common elements (not too far from the merge problem the other day).

    N = 100000
    def f(a): return a * (a+N)
    def g(b): return b * (b-1)
    a = N/10
    b = 1
    total = 0
    while a < N and b < N:
        if f(a) < g(b):
            a += 1
        elif f(a) > g(b):
            b += 1
            n = a*N + b
            total += n
            a += 1; b += 1

    We can optimize this of course (eg. a must be even and b can be initialized to 1000 as well, also various strength reductions are possible).

  5. Ernie Gore said

    A little bit of pencil and paper work will show that the last digit of b(b-1) can only be 0,2, or 6. Consequently,
    a must end in 0, 4, or 6. That eliminates 70% of the numbers that must be tested.

  6. matthew said

    @Ernie: good point, also b is 4k or 4k+1.

    Another approach:, we can use the quadratic Diophantine equation solver (and it’s methods) at http://www.alpertron.com.ar/QUAD.HTM to solve x^2 – y^2 +10000x +1 = 0, which reduces the problem to finding x’ and y’ that satisfy (x’ + y’)(x’ – y’) = 9999999999

  7. matthew said

    Here’s the divisors approach coded up, seems to work though I haven’t proved it correct. I don’t understand why a should be even and b should be odd, though that always seems to be the case.

    The equation is a^2 – b^2 + 100000a + b = 0, we substitute x = 2a+100000 and y = 2b-1, and as above, the equation then simplifies to (x+y)(x-y) = 9999999999. It’s very similar to the torn numbers problem of a few weeks ago and like there, a more intelligent generation of all divisors would be possible:

    A = 100000
    N = A*A-1
    total = 0
    for i in range(1,A):
        if N%i == 0:
            j = N//i
            x,y = (j+i)//2, (j-i)//2
            a,b = (x-A)//2, (y+1)//2
            if A//10 <= a < A and 0 <= b < A:
                n = a*A+b
                total += n
  8. matthew said

    Incidentally, I don’t think all the numbers on the page 3 list are right, for example 620741003025 is lacking in excellence. Maybe you need a range check for b?

  9. I’ve gone a bit crazy on this problem and have almost completed the list of 30-digit excellent numbers, along with some proofs about the interesting patterns in them. As someone pointed out, the list of numbers presented here has several errors. You can see more at excellent nums.com.

  10. matthew said

    I’ll double you: 617972744134287544096427144466999931525969228137094506182988

  11. matthew said

    In fact, with the help of http://www.numberempire.com/numberfactorizer.php and http://wpedia.goo.ne.jp/enwiki/User:Toshio_Yamaguchi/Factorizations_of_POTPO_numbers, it’s not too hard to factorize 10^160-1 as [3, 3, 11, 17, 41, 73, 101, 137, 271, 353, 449, 641, 1409, 3541, 9091, 27961, 69857, 1634881, 1676321, 5070721, 5882353, 18453761, 947147262401, 117633133855156294075039863860622161, 349954396040122577928041596214187605761], and once we have that, we can find (with a version of my python program up above) eg:


    along with another 153843 most excellent 160-digit numbers.

  12. matthew said

    Thanks to this excellent website, http://stdkmd.com/nrr/repunit/, we can now go to 2016 digits, for example:

  13. […] number n, with an even number of digits, is excellent if it can be split into two halves, a and b, such that b2 – a2 = n. Let 2k be the number of digits, […]

  14. matthew said

    I was going to add a link to my blog post working out the divisors solution properly, but WordPress seems to have done that for me.

    One curious fact: after sorting the numbers produced lexicographically, the largest that came up was:


    Looks like we are converging on the Golden Ratio!

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s

%d bloggers like this: