Enigmatic Code

Programming Enigma Puzzles

Tantalizer 451: Death rates

From New Scientist #1002, 27th May 1976 [link]

A minor problem which has long troubled medical historians is why two seemingly identical Edwardian TB sanataria in Mercia should have strikingly different death rates. The answer turns out to be simple enough — St. Bede’s in fact had 3,185 more patients than St. Crispin’s.

How did this come to be overlooked? Well, the probable reason is that both were built on the same formula. Each had as many wings as it had wards in each wing and had as many wards in each wing as it had patients in each ward. This was so plainly the whim of some mad bureaucrat that no historian troubled to check whether the number of patients per ward was the same in each place.

So what is the figure in each case?


7 responses to “Tantalizer 451: Death rates

  1. Jim Randell 6 June 2018 at 8:04 am

    If the characteristic numbers of St. Bede’s and St. Crispin’s are b and c respectively, then we are looking to solve:

    b³ – c³ = 3185

    This Python program looks at increasing values of c and the corresponding value of ³√(c³ + 3185).

    It runs in 68ms.

    Run: [ @repl.it ]

    from itertools import count
    from enigma import arg, iroot, printf
    N = arg(3185, 0, int)
    for c in count(1):
      c3 = c ** 3
      b3 = c3 + 3185
      b = iroot(b3, 3)
      if b == c: break
      if b ** 3 == b3:
        printf("{N} = {b3} - {c3} = {b}^3 - {c}^3")

    Solution: C has 1728 patients (= 12³), B has 4913 patients (= 17³).

    Most of the reported run time is taken up in housekeeping of the Python interpreter. But we can use the Timer() class from the enigma.py library to find the actual elapsed time of a chunk of code, and doing that we find the internal run time of this program is about 107µs. I came with a couple of other approaches, but they both clocked in with a similar, but slightly slower, internal run time of around 120µs.

  2. Hugh Casement 6 June 2018 at 9:08 am

    This puzzle strikes me as more trivial than most.
    b must be at least 15 because 14³ < 3185 and cannot be more than 33 because 34³ – 33³ > 3185.
    Even I could write a very short program in Basic, that would find the solution very quickly (without the need to take cube roots).
    Though I have to say that run time is quite unimportant for something like this that is executed only once and is not a real-time control application. The time it takes to type in the program is inevitably very much greater!

  3. Tessa Fullwood 6 June 2018 at 10:23 am

    Minor problem… so I agree, rather more trivial than most.

  4. Hugh Casement 6 June 2018 at 12:25 pm

    Thank you for fixing the > and <, Jim.

    It’s since occurred to me that b³ – c³ = (b – c)(b² + bc + c²); and 3185 has prime factors 5, 7, and 13, so (b – c) must be one of those, which reduces the possibilities.  If we tentatively take 31 as the upper limit for b, we can use single-length integer arithmetic, insert an inner loop for c, and avoid taking cube roots.

  5. arthurvause 6 June 2018 at 4:41 pm

    Using Hugh’s factorisation of b^3-c^3:

    Suppose the number of patients of St Crispins and St Bedes are x, x+n
    then (x+n)^3-x^3 = 3185 = 3n*x^2 +3n^2*x +n^3 = n(3x^2 +3nx +n^2)
    3185 = 5*7^2*13, so choosing a value of n to be a (small) divisor of 3185: 5,7,13,35,49
    3x^2 + 3nx + n^2-3185/n = 0
    x = (-3n + sqrt( 9n^2 - 12*(n*n- 3185/n))/6
      = (-3n + sqrt(12*3185/n - 3*n*n))/6
    squares = {n*n:n for n in xrange(196)}
    for n in (5,7,13,35,49):
      d = 12*3185/n - 3*n*n
      if d in squares:
        xnum = -3*n + squares[d]
        if xnum%6==0:
          x = xnum/6
          print "patients per ward is ",x,x+n
  6. Jim Randell 6 June 2018 at 11:36 pm

    Here’s a solution that uses factorisation to solve the problem for a general difference in cubes:

    b³ – c³ = N


    b³ – c³ = (b – c)(b² + bc + c²) = N

    We can consider divisor pairs p and q of N such that p×q = N, and p < q.


    p = b – c
    q = b² + bc + c²


    q – p² = (b² + bc + c²) – (b² – 2bc + c²) = 3bc

    so, if (q – p²) is positive and divisible by 3 we can determine candidate values for b and c directly by looking at the divisor pairs of (q – p²) / 3.

    from enigma import arg, divisors_pairs, printf
    N = arg(3185, 0, int)
    for (p, q) in divisors_pairs(N):
      (d, r) = divmod(q - p ** 2, 3)
      if d < 0: break
      if r != 0: continue
      for (c, b) in divisors_pairs(d):
        if b - c == p:
          printf("{N} = {b3} - {c3} = {b}^3 - {c}^3", c3=c ** 3, b3=b ** 3)
    • Jim Randell 7 June 2018 at 7:08 am

      In fact if we have an integer d:

      d = (q – p²) / 3

      Then we have equations:

      bc = d
      b – c = p

      Solving for c gives the following quadratic:

      c² + pc – d = 0

      Which, using the quadratic formula has a solution at:

      c = (√(p² + 4d) – p) / 2

      and substituting for d:

      c = (√((4q – p²) / 3) – p) / 2
      b = p + c

      Which gives a similar approach to Arthur’s method above – factorising N and then solving a quadratic equation.

      from enigma import arg, divisors_pairs, is_square, printf
      N = arg(3185, 0, int)
      for (p, q) in divisors_pairs(N):
        (d, r) = divmod(4 * q - p * p, 3)
        if d < 0: break
        if r != 0: continue
        x = is_square(d)
        if x is None: continue
        (c, r) = divmod(x - p, 2)
        if c < 1 or r != 0: continue
        b = c + p
        printf("{N} = {b3} - {c3} = {b}^3 - {c}^3", c3=c ** 3, b3=b ** 3)

Leave a Comment

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 )

Google photo

You are commenting using your Google 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

This site uses Akismet to reduce spam. Learn how your comment data is processed.

%d bloggers like this: