A 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, then we want n = aA + b = b2 - a2 where A = 10k.

So, every 2k digit excellent number gives rise to divisors i,j of N where ij = N and i <= j

This process can be reversed: if i is a divisor of N, with j = N/i and i <= j, we have X = (j+i)/2, Y = (j-i)/2, then a = (X-A)/2 and b = (Y+1)/2. If all the divisions by 2 are exact (and in this case they are – N is odd, so i and j are too, also writing i = 2i'+1 and j = 2j'+1, we can show that i' and j' must have different parities) then we have a potentially excellent number – all we need to check is that a has exactly k digits and that b has at most k (otherwise a and b are not the upper and lower halves of a 2k digit number).

Now we have a nice algorithm: find all divisors i of N = 10k-1, with i <= sqrt(N), find a and b as above and check if they are in the appropriate range, if so, we have an excellent number and it should be clear that all excellent numbers can be generated in this way.

For small N, we can find all divisors just by a linear scan, but for larger N something better is needed: given a prime factorization we can generate all possible combinations of the factors to get the divisors, so now all we need to do is factorize 102k-1. This of course is a hard problem but we can use, for example, Python’s primefac library, and give it some help by observing that 102k-1 = (10k-1)(10k+1). The factorization is harder for some values of k, particularly if k is prime, but we can always have a look at:

if we run in to trouble. My Pi 2 gets stuck at k = 71 where 11, 290249, 241573142393627673576957439049, 45994811347886846310221728895223034301839 and 31321069464181068355415209323405389541706979493156189716729115659 are the factors needed, so it’s not surprising it is struggling. Also, the number of divisors to check is approximately 2n-1 where n is the number of prime factors, of which, for example 1090-1 has 35 so just generating all potential 180 digit numbers will take a while.

So, after all that, here’s some code. Using Python generators keeps the memory usage down – we can process each divisor as it is constructed, (though it does mean that results for a particular size don’t come out in order) – after running for around 24 hours on a Pi 2, we are up to 180 digits and around 2000000 numbers but top reports less than 1% of memory in use.

There is, of course more: notably the matrix based solutions – we observe that (an+2,an+1) is M(an+1,an) for 2×2 matrix M = (1 1 1 0). We then have (Fn+1,Fn) = Mn(1 0) and we can use conventional matrix operations, including inversion, to compute Fn. We can apply this to any linear recurrence, eg. an+2 = Aan+1+Ban corresponds to the matrix (A B 1 0).

It’s notable that M only depends on the recurrence and not on the initial values, so we can compute, for example, the Lucas numbers with (Ln+1,Ln) = Mn(1 2), and there is a regularity to the structure of M that can be used to optimise the computation, eg. for Fibonacci, M is of the form (a+b a a b) (cf. HAKMEM: http://www.inwap.com/pdp10/hbaker/hakmem/recurrence.html). But that probably needs another post.

Often, I’d like to embed a reasonably capable command interpreter in a C++ application. Python seems a likely candidate, so here’s some investigative code using separate processes (the next step will be to use threads, if that’s possible, so the interpreter can live in the same memory space as our application, that can wait for part II though). As well as the mechanics of embedding Python, we have a pleasant excursion through the sometimes murky worlds of signal handling and pseudo-terminals.

The server structure is conventional (though not necessarily suitable for a serious production server), on each incoming connection we fork a handler process, this in turn splits into two processes, which form their own process group under the control of a pseudo-terminal (pty). One forwarding process copies data between the socket and the master side of the pty, the other process runs the interpreter itself on the slave side. Simple enough, with a few subtleties. To get signal handling right, we have to ignore SIGINT in the forwarding process (otherwise it will terminate on interrupt, taking the interpreter with it), but leave the default handler in the interpreter process – Python sets up its own signal handler, but it only seems to do this if the handler hasn’t been redefined already. Also, Python seems to insist that it uses fds 0,1 and 2 so we need to rebind them, and, finally, to get Python to do line editing, we need to import readline in the interpreter.

My main interest here is in getting external access to the interpreter, rather than the mechanics of calling between C and Python, so we just have a couple of simple functions init() and func() defined in the embedded interpreter as examples. At this simple level I don’t think we need to worry about reference counts etc.

Of course, all this will depend on your exact Python version and where it is installed. Embedding has changed somewhat in Python 3, but most of this will still apply.

To connect to the interpreter, we can use our good friend netcat, with some extra tty mangling (we want eg. control-C to be handled by the pty defined above in the server code, not the user terminal, so we put that into raw mode).