Excellent Numbers
March 24, 2015
Today’s exercise channels our inner Project Euler:
An excellent number n has an even number of digits and, if you split the number into the front half a and the back half b, then b 2 − a 2 = n. For example, 3468 = 682 − 342 = 4624 − 1156 = 3468, so 3468 is an excellent number. The only two-digit excellent number is 48 and the only four-digit excellent number is 3468. There are eight six-digit excellent numbers, 140400, 190476, 216513, 300625, 334668, 416768, 484848, and 530901, and their sum is 2615199. What is the sum of the 10-digit excellent numbers?
Your task is to compute the sum of the 10-digit excellent numbers; in the spirit of Project Euler, your solution should take no more than one minute of computation time. When you are finished, you are welcome to read or run a suggested solution, or to post your own solution or discuss the exercise in the comments below.
Programming Praxis 
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:
Or in one line (rewriting b):
Both result in 13171037992
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))) #13171037992Strictly 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.
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 else: n = a*N + b total += n print(n) a += 1; b += 1 print(total)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).
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.
@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
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 print(n) total += n print(total)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?
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.
I’ll double you: 617972744134287544096427144466999931525969228137094506182988
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:
3016693653810692443471528771919068734283766879102218756224920866230260139596321662663687362816936566086878697389529438778910186253288857357094019024639046138113
along with another 153843 most excellent 160-digit numbers.
Thanks to this excellent website, http://stdkmd.com/nrr/repunit/, we can now go to 2016 digits, for example:
[…] 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, […]
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:
618033865994856576190462466045494275576755186660999999862755693037993862789961848095236722795025
Looks like we are converging on the Golden Ratio!