IBM Ponder This, January 2017

The IBM Ponder This puzzle for January 2017 asked for a rational solution for x^y=y^x, where the last eight digits of x are different from each other, and for a bonus asked for a solution where the last eight digits of x are a permutation of 1,2,3,4,5,6,7,8.

It is possible to do better than this, and construct a solution where the last nine digits of x are are a permutation of 1,2,3,4,5,6,7,8,9 , and a solution where the last ten digits of x are are a permutation of  0,1,2,3,4,5,6,7,8,9.

Original Problem and Bonus

Goldbach and Euler proved that a set of 2 rationals {x,y} satisfies x^y=y^x    iff   $ \{x,y\} = \{ (\frac{n+1}{n})^n , (\frac{n+1}{n})^{n+1} \} $  for integer n. [1], [2]

If these numbers are to have last digits, n must be of the form $ n = 2^t 5^f $

$ \frac{n+1}{n}$ has $ d = max(t,f) $ decimal places.

$ \frac{(n+1)10^d}{n} $ is an integer whose last digits are the last decimal digits of $ \frac{n+1}{n}$

The last $D$ digits of $x,y$ are 

$ (\frac{(n+1)10^d}{n})^{n+k} $ mod $ 10^D $,  for $k = 0,1 $                  

Which can be calculated as

$ ((n+1)2^{d-t}5^{d-f})^{n+k} $ mod $ 10^D$,  for $k = 0,1 $

or

$ (10^d +2^{d-t}5^{d-f})^{n+k} $ mod $ 10^D$,  for $k = 0,1 $                            ...... (1)

For modest values of n, it is easy to calculate these residues using modular exponentiation, for example in C or Python.  This Python program searches for solutions to both problems.

The first problem has a solution for $n=2^4 \cdot 5^1 =80$, $x =( \frac{81}{80})^{80} $

The bonus problem has a solution for $n= 5^5 = 3125$, $x =( \frac{3126}{3125})^{3126} $

Extending the Problems

Now the more interesting problems: find a solution where the last nine digits of x are are a permutation of 1,2,3,4,5,6,7,8,9 , and a solution where the last ten digits of x are are a permutation of  0,1,2,3,4,5,6,7,8,9.

Here we need to examine larger values of $n$,  so expression (1) becomes impractical to use and we have to find a better way of determining the last digits of $ (\frac{(n+1)10^d}{n})^{n+k} $

If $d \ge  D$, then (1) can be reduced to

$ (2^{d-t}5^{d-f})^{n+k} $ mod $ 10^D$,  for $k = 0,1 $ 

There are two cases to consider:

$f>t$

$ (2^{d-t}5^{d-f})^{n+k} $ mod $ 10^D  = ( 2^{f-t})^{n+k} $ mod $ 10^D =  2^{(f-t)(n+k)} $ mod $ 10^D $

Now [3] shows that powers of 2 mod $10^D$ cycle with a period of $P = 4 \cdot 5^{D-1}$, starting at $2^D$

so if $ (f-t)(n+k) = (f-t)(2^t 5^f+k) \ge D $,

$2^{(f-t)(n+k)} $ mod $ 10^D $ =  $ 2^{ D + ((f-t)(2^t 5^f+k) - D) \: mod  \: P } $ mod $ 10^D $                ..... (2)

$t >f$

$ (2^{d-t}5^{d-f})^{n+k} $ mod $ 10^D = ( 5^{t-f})^{n+k} $ mod $ 10^D =  5^{(t-f)(n+k)} $ mod $ 10^D $

 [4] shows that powers of 5 mod $10^D$ cycle with a period of $P=2^{D-2}$, starting at $5^D$ 

so if $ (t-f)(n+k) = (t-f)(2^t 5^f+k) \ge D $,

$5^{(t-f)(n+k)} $ mod $ 10^D $ =  $ 5^{ D + ((t-f)(2^t 5^f+k) - D)   \: mod  \: P } $ mod $ 10^D $                ..... (3)

The residues in (2) and (3) can be calculated using modular exponentiation without the need to handle excessively large numbers. This Python program searches for solutions to both problems, finding the following solutions:

For $n = 5^{267}$, the last ten digits of $( \frac{n+1}{n})^{n+1} $ are 3628975104

For $n = 2^2 \cdot 5^{957}$, the last nine digits of $( \frac{n+1}{n})^{n+1} $ are 247315968

[1] Cut The Knot, Rational Solutions to x^y = y^x

[2] Sved, Marta. “On the Rational Solutions of x^y = y^x.” Mathematics Magazine, vol. 63, no. 1, 1990, pp. 30–33., http://www.jstor.org/stable/2691508.

[3] Exploring Binary, Cycle Length of Powers of Two Mod Powers of Ten, Rick Regan, November 6th, 2009

[4]  Exploring Binary, Cycle Length of Powers of Five Mod Powers of Ten,  Rick Regan, December 22nd, 2009

Cheryl's Birthday

The puzzle, originally set in the 2015 Singapore and Asian Schools Math Olympiad, know as Cheryl's Birthday, has generated a lot of interest over the past few weeks.

A discussion of the puzzle here queried the meaning of the term "day" in the statement "Cheryl then tells Albert and Bernard separtely the month and day of her birthday respectively".

The generlly accepted interpretation is that it is the day of the month, but I wondered whether the puzzle could be solved if "day" were interpreted as the day of the week, as suggested in the discussion.

The solution, assuming the "day of the month" interpretation, looks like this:

The days of the week associated with each date in 2015 are:

 

Now suppose Albert and Bernard are told all these dates, Albert is told the month of Cheryl's birthday and Bernard is told the day of the week. A solution can still be derived, as follows:

A solution with this structure would work for any other year, as the days of the week would all be offset by the same amount.

Factoring numbers in Python and C++

While trying to solve the "Squares on a Cube" puzzle on Contest Center, I became sidetracked into writing an efficient program to produce the prime factors of a number, an exercise that most programmers will have attempted at some time.

The hybrid Python / C++ code presented here combines well known techniques, Miller Rabin and Brent's variation of Pollard Rho, together with a modern fast technique for modular reduction developed by Rutten and Eekelen.

Python implementation

After a bit of experimenting, I settled on this mixture for a Python implementation:

• Lookup tables for numbers below 10⁶
• Trial division of primes below 300
• Miller-Rabin primality test (deterministic up to the published limit)
• Brent's variant of the Pollard Rho algorithm, by default running for up to 1 second
• Elliptic Curve Method (ECM) if the Brent algorithm fails to find a factor within its time limit. 

C++ implementation

While the Python version is reasonably fast, a C++ version  is much faster. My C++ code factors numbers below 2⁶⁴, where Brent Pollard Rho can be used without having to use ECM or Multiple Polynomial Quadratic Sieve.

The performance of Brent Pollard Rho and Miller Rabin is ultimately dominated by modular multiplication, so any improvement of modular multiplication time is reflected in the overall performance of the factorisation algorithm. Here are my attempts to speed up modular multiplication:

I. Shift and add

To calculate $a \times b$ mod $m$  I initially used an intricate, and relatively slow, shift and add method. Several implementations of this can be found on the internet. The version I used was stolen from Craig McQueen's post on stackexchange.  This version is valid for any $ m < 2^{64} $

Function modmulSA implements this technique.

II, Repeated reduction of most significant 64 bits

  My first attempt to speed up modular multiplication, valid for a modulus $m < 2^{63}$,  uses the following techniques:

a) Division by multiplication

If a modulus of the same denominator needs to be obtained many times, as in the Brent algorithm and the Miller Rabin test, an efficient  technique is to construct an inverse $M=2^k/m$ of the denominator $m$ once (M is often termed a "magic number"), then perform division by $m$ by multiplying the numerator by M and bit shifting to divide by $2^k$.  Of course it's a bit more complicated than that, and is described comprehensively by Henry S. Warren Jr in Hackers Delight, (Chapter 10, Integer Division by Constants), and explained in a readable form by Ridiculous Fish in Labor of Division (Episode I).

C++ functions magicu2 and magicdiv (stolen from Hackers Delight) and magicmod implement these techniques.

b) Modular multiplication

1. Reduce a and b modulo m (using the division by multiplication technique)
2. Construct the product $c = a \times b$. If $a$ and $b$ are 64 bit numbers, $c$ is a 128 bit number. I used the MS VC++ intrinsic function for this calculation, but if this is not available, Hackers Delight function mulul64 can be used.
3. Repeatedly find the most significant 64 bits of $c$ and reduce them modulo $m$ (again using division by multiplication), until c is reduced to a number less than $m$

Function modmulRR implements this technique.

III, Assembler

On 64 bit Intel processors, the multiply instruction multiplies two 64 bit numbers and stores the 128 bit result in two registers, rdx:rax. The divide instruction divides a 128 bit number stored in rdx:rax by a 64 bit number and stores the 64 bit quotient in rax and the 64 bit remainder in rdx.

So these instruction are ideal for modular multiplication. However, the intrinsic function to perform division of a 128 bit number is not available in MS Visual Studio, so a few lines of assembler (mod64.asm) have to be written.

The assembler version has the advantage of being valid for a modulus up to $2^{64}-1$, and is noticeably faster than repeated reduction, and of course much faster than the shift and add method.

Function modmulAS is the C interface to this function.

IV. Rutten-Eekelen

An even faster method of modular multiplication than assembler, for a modulus $m<2^{63}$, uses the algorithm described in Efficient and formally proven reduction of large integers by small moduli by Luc M. W. J. Rutten and Marko C. J. D. van Eekelen. This method requires the pre-calculation of a similar, but slightly different, magic number to the one described above, and then makes use of this number in an intricate algorithm whose correctness has been verified using a computer generated proof.

Function modmulRE implements this function.

Performance

This graph shows the average time to factor numbers in the range $10^n$ to $10^n+10000$ (giving a representative mix of easy and hard factorisations) for all my C++ versions, and for reference, a Python implementation of Brent Pollard and Miller Rabin. The graph shows that the shift and add version has performance no better than the Python version, and that the Rutten-Eekelen version beats the Assembler version, for numbers below $2^{63}$

C++ / Python integration

The ctypes module integrates Python with a DLL of the C++ implementation. Calling a DLL from Python  incurs an overhead of about 1.5µs, so it is faster to use the Python implementation for small numbers that are factored using a lookup table. For larger numbers, the overhead soon becomes negligible compared to the time to factorise the number.

Rutten-Eekelen has the best performance for a modulus $m<2^{63}$, even beating the assembler version. For a modulus $2^{63} \le m \lt 2^{64}$ the assembler version is fastest, so the final C++ version called from Python uses a hybrid of these two methods.

To verify that the Brent Pollard Rho algorithm copes with difficult cases, this graph shows the average and worst case times to factorise numbers within $ \pm 5000$ of various baselines. Note the jump in worst case factorisation times for numbers $ \ge 2^{63}$, where the algorithm switches from Rutten-Eekelen to Assembler.

The worst cases are listed here.

Source Code 

All the source code is located here, except for the optional ECM algorithm, for which I used the Sourceforge implementation pyecm.py.

A stand-alone Python module, factors.py can be run without the C++ DLL or the ECM module being present. 
This module has the following principal functions:

Description	Example call	Return value
all_factors(n) returns a list of all the factors of n The function chooses whether to use the Python or C++ version depending on the size of n.	all_factors(12)	[1,2,3,4,6,12]
totient(n) returns the Euler totient function	totient(30)	8
factorise(n, maxBrentTime=1.0) returns a list of the prime factors of n maxBrentTime determines how many seconds the Brent algorithm should run before reverting to ECM  (if available). The function chooses whether to use the Python or C++ version depending on the size of n.	factorise(12)	[2,2,3]

 
A DLL, FactorsDLL.dll,  built on Windows 7 64 bit, is included, together with a Python module factorsDLL.py to call the DLL from Python. I'm not sure whether the DLL will work on other machines. In case it doesn't, the source code to build the DLL is also included.

 Implementation notes

The Miller Rabin test can occasionally derive a factor of a composite number, when testing witness $a$ for the primality of $n$ encouters the situation

$a^{2^rd} \ne \pm 1$ mod $n$  and  $a^{2^{r+1}d} =1$ mod $n$

In this situation  $gcd(a^{2^rd}-1,n)$ produces a factor of $n$. I have included this test, even though the condition rarely occurs in practice.

The Brent Pollard Rho algorithm occasionally fails to find a factor of a composite, so I simply run the algorithm again. The random selection of parameters seems to ensure that the algorithm does not have to be run too many times to obtain a factor.

Postscript

Even though the program vastly reduced the time to factorise a number, I still haven't solved the Contest Center problem that motivated me to write it.

Calculating square roots on a pinwheel calculator

Before electronic calculators became affordable, complex calculations were often done using pinwheel calculators. These were manufactured by several companies, for example, Odhner, Brunsviga, Thales, Schubert, all based on designs of Frank Baldwin and Willgodt Odhner in the 19th century.

Original Odhner 127

Original Odhner 227

Brunsviga 20

Having recently acquired three of these calculators, I set about investigating what they were capable of. As well as the four basic arithmetic operations, addition, subtraction, multiplication and division, square roots can also be calculated. The instruction manual for the Original Odhner 227 describes the procedure to extract the square root of 966289. At first glance, the instructions look like an alchemist's spell, but they are performing the following calculations to derive √966289 = 983:

From Odhner 227 instructions (click to enlarge)

c)	966289  - 10000 - 30000 - 50000 - 70000 - 90000 -110000 - 130000 - 150000 - 170000 = 156289	( 9 subtractions)
d) & e)	156289  - 18100 - 18300 - 18500 - 18700 - 18900 - 19100 - 19300 - 19500 = 5889	( 8 subtractions)
f)	5889 -  1961 - 1963 -1965 = 0	( 3 subtractions)

 The algorithm is based on the fact that the square of n is the sum of the first n odd numbers:

12 = 1     =     1 22 = 4     =     1+3 32 = 9     =     1+3+5 42 = 16    =     1+3+5+7 52 = 25    =     1+3+5+7+9 etc.

The square root of a number can therefore be calculated by subtracting successive odd numbers until the sum is zero. The square root is the number of subtractions performed. For example, to find the square root of 64, successively subtract 1,3,5,... until zero is reached. The number of subtractions required, in this case 8, is the square root of the number. For large numbers, e.g. 966289, this would obviously be impractical, but we can attack large numbers one pair of digits at a time, to produce one digit of the square root at a time:

64 -  1 = 63 63 -  3 = 60 60 -  5 = 55 55 -  7 = 48 48 -  9 = 39 39 - 11 = 28 28 - 13 = 15 15 - 15 =  0

To find the first digit of the square root of 966289, subtract odd ten thousands: 10000, 30000, 50000, ...  until the point just before the result would turn negative, so subtract 10000+30000+...+170000.  Now, 1+3+...+17 = 9² , so 10000+30000+...+170000 = 900² . The first digit of the square root is 9.

To find the first two digits of the square root, we could start with 966289 and subtract odd hundreds: 100, 300, 500,... up to 19500. However, we have already subtracted 900², which is equal to 100+300+500+...+17700+17900, so we can start subtracting at 18100, subtracting only 18100+18300+...+19500.  We have now subtracted a total of 980².

Similarly to find the first three digits of the square root (i.e. the complete square root), we could start with 966289 and subtract 1, 3, 5,...up to 1965. But we have already subtracted 980² = 1+3+5+...+1959, so we can start subtracting at 1961, only needing to subtract 1961+1963+1965

10000+30000+...+170000		= 9002
18100+18300+...+19500	= (100+300+...+19500) - (100+300+...+17900)	= 9802 - 9002
1961+1963+1965	= (1+3+...+1965) - (1+3+...+1959)	= 9832 - 9802

^{The algorithm can also be used to calculate approximations to the square root of numbers that are not perfect squares. Here is a video by James Grime (singingbanana) and Hugh Hunt that demonstrates an Original Odhner being used to calculate √2.}

Feeler Gauges and Magic Circles

I recently came across Sunday Times Teaser puzzle, 777 on Brian Gladman's Puzzling In Python site:

I have half-a-dozen feeler gauges on a ring like keys on a key ring. They are attached in a certain order and by selecting a single gauge or combinations of adjacent gauges it is possible to measure from 1 to 31 thousandths of an inch (thou) inclusive in steps of 1 thou. If I remove two adjacent gauges, replace them with a single gauge and rearrange the five on the ring I can now measure from 1 to 21 thou inclusive in steps of 1 thou with these five gauges by again selecting a single gauge or combinations of adjacent gauges.

What are the thicknesses of the gauges and the order of arrangement on the ring in each case? (start with gauge of unit thickness and then its thinner neighbour for the sets).

The puzzle is based on some interesting but relatively obscure mathematics, Perfect Difference Sets, described in A Theorem in Finite Projective Geometry and Some Applications to Number Theory, James Singer, Transactions of the American Mathematical Society Vol. 43, No. 3 (May, 1938), pp. 377-385, but initially analysed by the splendidly titled Rev. Thomas P. Kirkman, A.M., F.R.S., Rector of Croft with Southworth, and Honorary Member of the Literary and Philosophical Society of Manchester, in a paper entitled On the Perfect r-Partitions of $r^2-r+1$, Transactions of the Historical Society of Lancashire and Cheshire, vol 9 (1857), pp 127-142.

The puzzle is slightly contrived, as in reality feeler gauges allow non-adjacent gauges to be combined. The restriction of using only adjacent gauges makes the puzzle equivalent to finding Perfect Difference Sets of length 5 and 6.

Perfect Difference Sets

A Perfect Difference Set (PDS)  is a set of $k+1$ residues ${a_1,...,a_{k+1}}$ modulo $n=k^2+k+1$ such that every nonzero residue can be uniquely expressed in the form $a_i - a_j$

For example, the set {0,1,4,6} is a perfect difference set for $n=3^2+3+1=13$.

A PDS can also be expressed as the sequence of gaps between successive numbers in the sorted PDS with $n=k^2+k+1$ added to the set. For the PDS {0,1,4,6} the sequence of gaps is 1,3,2,7. If the sequence is wrapped round (as in feeler gauges on a ring), sums of adjacent numbers produce all the numbers from 1 to 13. These sequences will be referred to as magic circles.

Magic Circles

The article Magic Circles, David A. James, Mathematics Magazine, Vol. 54, No. 3 (May, 1981), pp. 122-125, describes a related puzzle, where Evariste Galois challenges 18 soldiers to arrange themselves around a circular pond with circumference 307 (=$17^2+17+1$) paces so that the number of paces between any two soldiers is unique. The soldiers would of course form a PDS to achieve this.

The article describes an algorithm to construct a PDS for $p^s+1$ numbers (prime p), using Galois fields (a.k.a. finite fields), and explains the mathematics behind the algorithm. The steps are:

1. Construct the Galois Field $GF(p^{3s})$, by finding an irreducible polynomial over GF(p) of degree $3s$.
   
   For s=1 this is straightforward, as a reducible polynomial of degree 3 always has a factor of degree one, so a polynomial $I(x)$ of degree 3 is irreducible if and only if $ I(x) \ne 0$ mod p for $x = 0,...,p-1$.
   
   For $s >  1$, finding an irreducible polynomial of degree 3s is more involved, so I cheated and pre-programed some that I obtained from the University of Victoria's Polynomial Generator
2. Find a generator, $g$ of the related cyclic group $GF^*(p^{3s})$.
   
   This is best done by making random guesses. I tested guesses by brute force, checking that $g^i  \ne 1$ for $1 \leq i \leq p^{3s}/2$. Typically a generator is found in one or two guesses.


3. Find elements of the subgroup $GF^*(p^s)$ in $GF^*(p^{3s})$
   
   These are the elements $g^i$ for values of $i$ that are multiples of $p^{2s}+p^s+1$ = $\frac{p^{3s}-1}{p^s-1}$
4. Find integers $i$ that satisfy $g^i - g \in GF^*(p^s)$ .
These integers, modulo $p^{2s}+p^s+1$ , together with zero, form a PDS for $p^s+1$

This algorithm will find a single PDS. Other PDS's of the same order can be found by multiplying this PDS by each number between 1 and $ p^{2s}+p^s+1$ that is coprime to $p^{2s}+p^s+1$.

Each PDS found can be be coverted to a magic circle, which can then be arranged to start with 1 followed by its smaller neighbour. Different PDS's can result in the same magic circle when this rearrangement is applied. For all the cases that I have tested, each magic circle is created by $3s$ different PDSs.

In fact, Golay, in Notes on the Representation of 1,2, ……, N by differences, J. London Math. Soc. (1972) s2-4 (4): 729-734, claims that there are $ \frac{\phi(p^{2s}+p^s+1)}{3s}$ magic circles.

Python Code

Python code is located in this folder. The folder contains:

finitefield.py -  implements Galois fields, representing polynomials $a_0+a_1x+...+a_nx^n$ by a python list $[a_0,a_1,...,a_n]$. The module contains procedures to multiply, divide, add, subtract polynomials modulo a prime and to multiply polynomials modulo an irreducible. The module also contains procedures to find an irreducible polynomial (for a limited range of $p^n$) and to find a generator of  $GF^*(p^n)$

magiccircles.py - implements the algorithm described above and uses the algorithm to solve the Sunday Times Teaser, and to generate magic circles for a range of values including $17^1+1 =18$, Evarise Galois' soldiers.  My algorithm eliminates reflections, producing $ \frac{\phi(p^{2s}+p^s+1)}{6s}$ magic circles.

Three angles

A puzzle recently posted in Enigmatic Code,  Enigma 1329: Height of ignorance, reminded me of a related puzzle:  For three squares laid side by side as shown in the diagram, prove that  α = β + γ

Clearly  α = 45⁰. The obvious approach is to use multiple angle formulae, but there is a more elegant geometrical solution.

Click here to see the geometrical solution
The solution of this puzzle can be used to solve Enigma 1329.

There is also a simple solution using complex numbers which consists of 17 symbols: 

$(3+i)(2+i) = 5(1+i)$

Does this count as a haiku?

Countdown Numbers Game

Several people have already written programs based on the Numbers Game from the TV program Countdown, but anyway here is another one, written in Python.

As has been noted elsewhere, a set of 2 large + 4 small or 1 large + 5 small typically covers a large proportion of  target numbers (101 to 999), and there are typically many ways of expressing each target using some or all of the 6 numbers.

These sets each cover all targets from 101 to 999. The lists of solutions () contain all targets from 0 until a gap above 999 is reached:

 5   6   7   8   9   75  (five consecutive numbers) 

 2   5   7   8   9   10  (six small numbers) 

 4   7   7   10   25  100  (a repeated number) 

 6   7   7   9   9  100  is perhaps more noteworthy,  covering every target number except the number 367 

 8   9   25   50   75  100  is the best set containing all four large numbers, covering 885 targets  


The code constructs solutions for a few games that have appeared on the TV show:

•  3   6   25   50   75  100 with a target of  952 , the James Martin 952 game, which has 2 solutions, (3×75×(6+100)-50)÷25 which was found by the contestant, and 25+(6×75×(3+100))÷50.
  
•  1   4   8   9   10   50  with a target of  500  ,possibly the easiset game ever?  for which my program finds 21 distinct solutions .
  
•  Some games solved by George Ford , one of the best countdown contestants.

Most of the complexity of the program is to do with minimising brackets,
e.g. (((((6)+(7))-(9))×(7))×(9))-(100)=(6+7-9)×7×9-100
and avoiding solutions that are essentially the same, e.g. 152=(6+7-9)×7×9-100=7×(7-9+6)×9-100.

The program also excludes "frivolous" solutions that contain an expression that equates to zero, e.g. 500 = (9-8-1)×4+50×100, solutions that involve multiplying or dividing by 1, e.g. 500 = (9-8)×50×100 or 500 = 50×100÷1 and solutions that have other redundancies, e.g. 500 = 9+10×50-1-8. 

# -*- coding: cp1252 -*-
from itertools import combinations

def wrap(t, condition=True):
    return ''.join(["(",t,")"]) if condition else t

def expressions( A, processed=set() ):
  """ initially, A is an array of numbers

   In recursive calls, A is an array of tuples (n,t,o,x,y,h) where
   n is a number,
   t is a text expression of n,
   o is the last operation applied to x,y to produce n = xoy
   h is the history of the expression, containing all the numbers produced in the construction
     of the expression
  """
  if type(A[0])==int: # inital call, so convert to array of tuples
    A = [ (A[i],str(A[i]),"",0,0, frozenset((A[i],)) ) for i in xrange(len(A))]
    processed=set()
  else: # recursive call, so return the expression formed in the parent
    n,t,o,x,y,h = A[0]
    if (n,t) not in processed:
      yield A[0]
      processed.add( (n,t) )