Sunday Times Teaser 2630 – River Crossing

by Des MacHale

Two walkers were returning to their tent, a short distance south east of them on the other side of a five metre wide river running east to west.

They both minimised their distance while in the water, one taking the shortest possible route with this constraint and the other simply going south and then east.

The shorter of these two routes was three quarters of the longer one.

What was its length?

4 Comments Leave one →

brian gladman permalink

It is easy to solve this manually since it involves a simple quadratic equation. It is also trivial to solve in Python, but to make things more interesting, here is a Python solution using a Golden Section search for a function minimum.


# golden section function minimiser 
def gs_search(f, a, b, tol=1.0e-14):
  R = 0.5 * (5.0 ** 0.5 - 1.0)
  C = 1.0 - R
  x1, x2 = R * a + C * b, C * a + R * b
  f1, f2 = f(x1), f(x2)
  while b - a > tol:
    if f1 > f2:
      a, x1, x2 = x1, x2, C * x1 + R * b
      f1, f2 = f2, f(x2)
    else:
      b, x1, x2 = x2, R * a + C * x2, x1
      f1, f2 = f(x1), f1
  return x1 if f1 < f2 else x2

# Since the position of the river doesn't matter, put the
# starting point on its northern edge so that both routes 
# immediately cross it.   Let the distances due south and
# due east on the longest route be x.

def dif(x):
  # the distance going due south and then due east
  d1 = 2 * x
  # the distance going diagonally after crossing the river 
  d2 = 5 + (x ** 2 + (x - 5) ** 2) ** 0.5
  # return the squared difference of (3/4) d1 and d2
  return (0.75 * d1 - d2) ** 2

x = gs_search(dif, 0, 30)
print('shorter distance = {:.2f} metres'.format(1.5 * x))

# golden section function minimiser

def gs_search(f, a, b, tol=1.0e-14):

R = 0.5 * (5.0 ** 0.5 - 1.0)

C = 1.0 - R

x1, x2 = R * a + C * b, C * a + R * b

f1, f2 = f(x1), f(x2)

while b - a > tol:

if f1 > f2:

a, x1, x2 = x1, x2, C * x1 + R * b

f1, f2 = f2, f(x2)

else:

b, x1, x2 = x2, R * a + C * x2, x1

f1, f2 = f(x1), f1

return x1 if f1 < f2 else x2

# Since the position of the river doesn't matter, put the

# starting point on its northern edge so that both routes

# immediately cross it. Let the distances due south and

# due east on the longest route be x.

def dif(x):

# the distance going due south and then due east

d1 = 2 * x

# the distance going diagonally after crossing the river

d2 = 5 + (x ** 2 + (x - 5) ** 2) ** 0.5

# return the squared difference of (3/4) d1 and d2

return (0.75 * d1 - d2) ** 2

x = gs_search(dif, 0, 30)

print('shorter distance = {:.2f} metres'.format(1.5 * x))

The version above finds the minimum of the squared difference between the the shortest and three quaters of the longest route. But we can instead simply look for the root (i.e. without the squaring) using a Brent method root solver.

In this example the time difference is minimal but in more complex examples the root finding method will generally be much faster.


# Brent method root finder (based on
# the description on Wikipedia).
def bm_solve(fn, a, b, tol=1.0e-14, max_it=100):

  fa, fb = fn(a), fn(b) 
  if fa * fb >= 0:
    raise ValueError
  if abs(fa) < abs(fb):
    a, b = b, a
 
  c, fc, mf = a, fa, True
  for i in range(max_it):
    
    # found a root within tolerance
    if not fb or abs(a - b) <= tol:
      return b
    
    if fa != fc and fb != fc:
      # inverse quadratic interpolation
      s = ( a * fb * fc / (fa - fb) / (fa - fc) 
          + b * fa * fc / (fb - fa) / (fb - fc) 
          + c * fa * fb / (fc - fa) / (fc - fb) )
    else:
      # secant rule
      s = b - fb * (b - a) / (fb - fa)
      
    t = (3 * a + b) / 4
    mf = ( not (t < s < b or b < s < t) 
           or mf and 2 * abs(s - b) >= abs(b - c) 
           or not mf and 2 * abs(s - b) >= abs(c - d)
           or mf and abs(b - c) < tol
	   or not mf and abs(c - d) < tol )
    if mf:
      s = (a + b) / 2
    
    d, c, fc, fs = c, b, fb, fn(s)
    if fa * fs < 0:
      b, fb = s, fs
    else:
      a, fa = s, fs
    if abs(fa) < abs(fb):
      a, b, fa, fb = b, a, fb, fa
  else:
    raise ValueError

# Since the position of the river doesn't matter, put the
# starting point on its northern edge so that both routes 
# immediately cross it.   Let the distances due south and
# due east on the longest route be x.

def dif(x):
  # the distance going due south and then due east
  d1 = 2 * x
  # the distance going diagonally after crossing the river 
  d2 = 5 + (x ** 2 + (x - 5) ** 2) ** 0.5
  # return the squared difference of (3/4) d1 and d2
  return (0.75 * d1 - d2)

x = bm_solve(dif, 0, 30)
print('shorter distance = {:.2f} metres'.format(1.5 * x))

# Brent method root finder (based on

# the description on Wikipedia).

def bm_solve(fn, a, b, tol=1.0e-14, max_it=100):

fa, fb = fn(a), fn(b)

if fa * fb >= 0:

raise ValueError

if abs(fa) < abs(fb):

a, b = b, a

c, fc, mf = a, fa, True

for i in range(max_it):

# found a root within tolerance

if not fb or abs(a - b) <= tol:

return b

if fa != fc and fb != fc:

# inverse quadratic interpolation

s = ( a * fb * fc / (fa - fb) / (fa - fc)

+ b * fa * fc / (fb - fa) / (fb - fc)

+ c * fa * fb / (fc - fa) / (fc - fb) )

else:

# secant rule

s = b - fb * (b - a) / (fb - fa)

t = (3 * a + b) / 4

mf = ( not (t < s < b or b < s < t)

or mf and 2 * abs(s - b) >= abs(b - c)

or not mf and 2 * abs(s - b) >= abs(c - d)

or mf and abs(b - c) < tol

or not mf and abs(c - d) < tol )

if mf:

s = (a + b) / 2

d, c, fc, fs = c, b, fb, fn(s)

if fa * fs < 0:

b, fb = s, fs

else:

a, fa = s, fs

if abs(fa) < abs(fb):

a, b, fa, fb = b, a, fb, fa

else:

raise ValueError

# Since the position of the river doesn't matter, put the

# starting point on its northern edge so that both routes

# immediately cross it. Let the distances due south and

# due east on the longest route be x.

def dif(x):

# the distance going due south and then due east

d1 = 2 * x

# the distance going diagonally after crossing the river

d2 = 5 + (x ** 2 + (x - 5) ** 2) ** 0.5

# return the squared difference of (3/4) d1 and d2

return (0.75 * d1 - d2)

x = bm_solve(dif, 0, 30)

print('shorter distance = {:.2f} metres'.format(1.5 * x))

Here is a manual solution.

Since the position of the river doesn’t matter, put the starting point its northern edge so that both routes immediately cross it. Let the southerly and easterly distances on the longest route be $x$ metres.

So the distances on the longer and shorter routes are $2x$ and $5+\sqrt{x^2+(x-5)^2}$ respectively. Setting the shorter route equal to three quarters of the longer one now gives:

$$\begin{array} ((3x/2-5)^2=&x^2+(x-5)^2 \\ 9x^2/4-15x+25=&2x^2-10x+25 \\ x^2/4 – 5x=&0 \\ x(x-20)=&0 \end{array}$$

This gives x as 20 metres and the shorter of the two routes as 30 metres.

Frits permalink

bm_solve() will not find the root for some complex functions.

Python

def dif(x): return (x - 1.7)**17 - 2.45 + 3.7 * (3.1 - x)**9 x, i = bm_solve(dif, -1, 100)

1
2
3
4
5

def dif(x): return (x - 1.7)**17 - 2.45 + 3.7 * (3.1 - x)**9

x, i = bm_solve(dif, -1, 100)

Chandrupatla gives the correct x = 2.7541183874033877

Reply

Frits permalink

Another root finding method (code from the internet)


sign = lambda x: -1 if x < 0 else 0 if x == 0 else 1

def tracker(f, a=0):
  ''' decorates calls to f(x) to track each call '''
  i = [0]
  def g(x):
    i[0] += 1
    y = f(x)
    print(f"i={i[0]}, x={x}, y={y}, err={y-a}")
    return y-a
  return g

def cm_solve(f, x0, x1,verbose=False, eps_m=None, eps_a=None):
  # credit: Jason Sachs
  # https://www.embeddedrelated.com/showarticle/855.php

  # adapted from Chandrupatla's algorithm as described in Scherer
  # https://books.google.com/books?id=cC-8BAAAQBAJ&pg=PA95

  # Initialization
  b = x0
  a = x1
  c = x1
  fa = f(a)
  fb = f(b)
  fc = fa
  # Make sure we are bracketing a root
  assert sign(fa) * sign(fb) <= 0
  t = 0.5
  do_iqi = False

  # jms: some guesses for default values of the eps_m and eps_a settings
  # based on machine precision... not sure exactly what to do here
  #eps = np.finfo(float).eps
  eps = 1.0e-14
  if eps_m is None:
    eps_m = eps
  if eps_a is None:
    eps_a = 2*eps

  while True:
    # use t to linearly interpolate between a and b,
    # and evaluate this function as our newest estimate xt
    xt = a + t * (b - a)
    ft = f(xt)
    if verbose:
        print('%s t=%f xt=%f ft=%g a=%f b=%f c=%f' % (
           'IQI' if do_iqi else 'BIS', t, xt, ft, a, b, c))
    # update our history of the last few points so that
    # - a is the newest estimate (we're going to update it from xt)
    # - c and b get the preceding two estimates
    # - a and b maintain opposite signs for f(a) and f(b)
    if sign(ft) == sign(fa):
      c = a
      fc = fa
    else:
      c = b
      b = a
      fc = fb
      fb = fa
    a = xt
    fa = ft
    # set xm so that f(xm) is the minimum magnitude of f(a) and f(b)
    if abs(fa) < abs(fb):
      xm = a
      fm = fa
    else:
      xm = b
      fm = fb
    if fm == 0:
      return xm
    # Figure out values xi and phi
    # to determine which method we should use next
    tol = 2 * eps_m * abs(xm) + eps_a
    tlim = tol / abs(b - c)
    if tlim > 0.5:
      return xm

    xi  = (a - b) / (c - b)
    phi = (fa - fb) / (fc - fb)
    do_iqi = phi**2 < xi and (1 - phi)**2 < 1 - xi

    if do_iqi:
      # inverse quadratic interpolation
      t = fa / (fb - fa) * fc / (fb - fc) + \
          (c - a) / (b - a) * fa / (fc - fa) * fb / (fc - fb)
    else:
      # bisection
      t = 0.5

    # limit to the range (tlim, 1-tlim)
    t = min(1 - tlim, max(tlim, t))

# Since the position of the river doesn't matter, put the
# starting point on its northern edge so that both routes
# immediately cross it.   Let the distances due south and
# due east on the longest route be x.

def dif(x):
  # the distance going due south and then due east
  d1 = 2 * x
  # the distance going diagonally after crossing the river
  d2 = 5 + (x ** 2 + (x - 5) ** 2) ** 0.5
  # return the squared difference of (3/4) d1 and d2
  return (0.75 * d1 - d2)


x = cm_solve(tracker(dif), 0, 40, verbose=True)
print('shorter distance = {:.2f} metres'.format(1.5 * x))

100

101

102

103

104

105

106

107

108

109

110

111

sign = lambda x: -1 if x < 0 else 0 if x == 0 else 1

def tracker(f, a=0):

''' decorates calls to f(x) to track each call '''

i = [0]

def g(x):

i[0] += 1

y = f(x)

print(f"i={i[0]}, x={x}, y={y}, err={y-a}")

return y-a

return g

def cm_solve(f, x0, x1,verbose=False, eps_m=None, eps_a=None):

# credit: Jason Sachs

# https://www.embeddedrelated.com/showarticle/855.php

# adapted from Chandrupatla's algorithm as described in Scherer

# https://books.google.com/books?id=cC-8BAAAQBAJ&pg=PA95

# Initialization

b = x0

a = x1

c = x1

fa = f(a)

fb = f(b)

fc = fa

# Make sure we are bracketing a root

assert sign(fa) * sign(fb) <= 0

t = 0.5

do_iqi = False

# jms: some guesses for default values of the eps_m and eps_a settings

# based on machine precision... not sure exactly what to do here

#eps = np.finfo(float).eps

eps = 1.0e-14

if eps_m is None:

eps_m = eps

if eps_a is None:

eps_a = 2*eps

while True:

# use t to linearly interpolate between a and b,

# and evaluate this function as our newest estimate xt

xt = a + t * (b - a)

ft = f(xt)

if verbose:

print('%s t=%f xt=%f ft=%g a=%f b=%f c=%f' % (

'IQI' if do_iqi else 'BIS', t, xt, ft, a, b, c))

# update our history of the last few points so that

# - a is the newest estimate (we're going to update it from xt)

# - c and b get the preceding two estimates

# - a and b maintain opposite signs for f(a) and f(b)

if sign(ft) == sign(fa):

c = a

fc = fa

else:

c = b

b = a

fc = fb

fb = fa

a = xt

fa = ft

# set xm so that f(xm) is the minimum magnitude of f(a) and f(b)

if abs(fa) < abs(fb):

xm = a

fm = fa

else:

xm = b

fm = fb

if fm == 0:

return xm

# Figure out values xi and phi

# to determine which method we should use next

tol = 2 * eps_m * abs(xm) + eps_a

tlim = tol / abs(b - c)

if tlim > 0.5:

return xm

xi = (a - b) / (c - b)

phi = (fa - fb) / (fc - fb)

do_iqi = phi**2 < xi and (1 - phi)**2 < 1 - xi

if do_iqi:

# inverse quadratic interpolation

t = fa / (fb - fa) * fc / (fb - fc) + \

(c - a) / (b - a) * fa / (fc - fa) * fb / (fc - fb)

else:

# bisection

t = 0.5

# limit to the range (tlim, 1-tlim)

t = min(1 - tlim, max(tlim, t))

# Since the position of the river doesn't matter, put the

# starting point on its northern edge so that both routes

# immediately cross it. Let the distances due south and

# due east on the longest route be x.

def dif(x):

# the distance going due south and then due east

d1 = 2 * x

# the distance going diagonally after crossing the river

d2 = 5 + (x ** 2 + (x - 5) ** 2) ** 0.5

# return the squared difference of (3/4) d1 and d2

return (0.75 * d1 - d2)

x = cm_solve(tracker(dif), 0, 40, verbose=True)

print('shorter distance = {:.2f} metres'.format(1.5 * x))

Frits permalink

And another…


from math import log
from sys import float_info

eps = float_info.epsilon

sign = lambda x: -1 if x < 0 else x > 0

# ITP "Interpolate, Truncate and Project" method root finder
# (based on the description on Wikipedia)
def itp_solve(fn, a, b, e=eps, k1=0.1, k2=2, n0=1):
  n_h = log((b - a) / (2 * e), 2)
  n_m = n_h + n0
  j = 0
  ya, yb = fn(a), fn(b)
  assert ya < 0 < yb, "incorrect values for a and b"
  assert k1 >= 0 and 1 <= k2 <= 1.5 + 5**.5 / 2, \
         "incorrect values for k1 and k2"

  while(b - a > 2 * e):
    # calculating parameters
    x_h = (a + b) / 2
    r = e * 2**(n_m - j) - (b - a) / 2
    d = k1 * (b - a)**k2
    # interpolation
    x_f = (yb * a - ya * b) / (yb - ya)
    # truncation
    s = sign(x_h - x_f)
    x_t = x_f + s * d if d <= abs(x_h - x_f) else x_h
    # projection
    x_ITP = x_t if abs(x_t - x_h) <= r else x_h - s * r
    # updating interval
    y_ITP = fn(x_ITP)
    if y_ITP > 0:
      b, yb = x_ITP, y_ITP
    elif y_ITP < 0:
      a, ya = x_ITP, y_ITP
    else:
      a, b = x_ITP, x_ITP
    j += 1

  return (a + b) / 2

def dif(x):
  # the distance going due south and then due east
  d1 = 2 * x
  # the distance going diagonally after crossing the river
  d2 = 5 + (x ** 2 + (x - 5) ** 2) ** 0.5
  # return the squared difference of (3/4) d1 and d2
  return (0.75 * d1 - d2)

x = itp_solve(dif, 0, 100)
print('shorter distance = {:.2f} metres'.format(1.5 * x))

from math import log

from sys import float_info

eps = float_info.epsilon

sign = lambda x: -1 if x < 0 else x > 0

# ITP "Interpolate, Truncate and Project" method root finder

# (based on the description on Wikipedia)

def itp_solve(fn, a, b, e=eps, k1=0.1, k2=2, n0=1):

n_h = log((b - a) / (2 * e), 2)

n_m = n_h + n0

j = 0

ya, yb = fn(a), fn(b)

assert ya < 0 < yb, "incorrect values for a and b"

assert k1 >= 0 and 1 <= k2 <= 1.5 + 5**.5 / 2, \

"incorrect values for k1 and k2"

while(b - a > 2 * e):

# calculating parameters

x_h = (a + b) / 2

r = e * 2**(n_m - j) - (b - a) / 2

d = k1 * (b - a)**k2

# interpolation

x_f = (yb * a - ya * b) / (yb - ya)

# truncation

s = sign(x_h - x_f)

x_t = x_f + s * d if d <= abs(x_h - x_f) else x_h

# projection

x_ITP = x_t if abs(x_t - x_h) <= r else x_h - s * r

# updating interval

y_ITP = fn(x_ITP)

if y_ITP > 0:

b, yb = x_ITP, y_ITP

elif y_ITP < 0:

a, ya = x_ITP, y_ITP

else:

a, b = x_ITP, x_ITP

j += 1

return (a + b) / 2

def dif(x):

# the distance going due south and then due east

d1 = 2 * x

# the distance going diagonally after crossing the river

d2 = 5 + (x ** 2 + (x - 5) ** 2) ** 0.5

# return the squared difference of (3/4) d1 and d2

return (0.75 * d1 - d2)

x = itp_solve(dif, 0, 100)

print('shorter distance = {:.2f} metres'.format(1.5 * x))

Sunday Times Teaser 2630 – River Crossing

by Des MacHale

Leave a Reply Cancel reply

Recent Posts

Recent Comments

Archives

Categories

Meta