Sunday Times Teaser 3295 – Always as the Crow Flies


#   A
#  D0F  x = D0, y = A0
#   S

max_sa, sol = 0, None
# consider AD, the shortest distance between any castles
for ad in range(10, 100):
  # consider AF and DS
  for af in range(ad + 1, 100):
    for ds in range(ad + 1, 100):
      # calculate FS from FS^2 = AF^2 + DS^2 - AD^2
      fs2 = af * af + ds * ds - ad * ad
      if (fs := round(fs2 ** (1 / 2))) ** 2 == fs2 and ad < fs < 100:
        # now consider the horizontal diagonal
        for df in range(ad + 1, 100):
          # calculate x = (DF^2 + AD^2 - AF^2) / (2.DF)
          if 0 < (x := (df * df + ad * ad - af * af) / (2 * df)) < ad:
            # calculate y from y^2 = (AD^2 - x^2)
            y = (ad * ad - x * x) ** (1 / 2)
            # calculate SA - y from (SA - y)^2 = DS^2 - x^2  and hence
            # SA which must have the maximum possible integer value
            sa_ = y + (ds * ds - x * x) ** (1 / 2)
            if ad < sa_ < 100 and abs(sa_ - (sa := round(sa_))) < 1e-10:
              tmp = (ad, af, ds, fs), (df, sa)                
              # all distances between ranches are different 
              if len(set(sum(tmp, ()))) == 6 and sa > max_sa:
                max_sa = sa
                sol = tmp
per, diags = sol            
print(f"Trip = {sum(per)} miles [{per}, {diags}]")

1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

25

26

27

28

29

30

31

32

# A

# D0F x = D0, y = A0

# S

max_sa, sol = 0, None

# consider AD, the shortest distance between any castles

for ad in range(10, 100):

# consider AF and DS

for af in range(ad + 1, 100):

for ds in range(ad + 1, 100):

# calculate FS from FS^2 = AF^2 + DS^2 - AD^2

fs2 = af * af + ds * ds - ad * ad

if (fs := round(fs2 ** (1 / 2))) ** 2 == fs2 and ad < fs < 100:

# now consider the horizontal diagonal

for df in range(ad + 1, 100):

# calculate x = (DF^2 + AD^2 - AF^2) / (2.DF)

if 0 < (x := (df * df + ad * ad - af * af) / (2 * df)) < ad:

# calculate y from y^2 = (AD^2 - x^2)

y = (ad * ad - x * x) ** (1 / 2)

# calculate SA - y from (SA - y)^2 = DS^2 - x^2 and hence

# SA which must have the maximum possible integer value

sa_ = y + (ds * ds - x * x) ** (1 / 2)

if ad < sa_ < 100 and abs(sa_ - (sa := round(sa_))) < 1e-10:

tmp = (ad, af, ds, fs), (df, sa)

# all distances between ranches are different

if len(set(sum(tmp, ()))) == 6 and sa > max_sa:

max_sa = sa

sol = tmp

per, diags = sol

print(f"Trip = {sum(per)} miles [{per}, {diags}]")

Reply

JohnZ permalink


from itertools import combinations
sq = tuple(x * x for x in range(100))
trg = []
# find right triangles with integral sides
for i in range(10, 98):
  for j in range(i + 1, 99):
    k2 = i * i + j * j
    if k2 in sq:
      k = sq.index(k2)
      trg. append((i, j, k))
 
perim = 0
for tr4 in combinations(trg, 4):
  # can the four triangles form a quadrilateral
  if len(set().union(*((a, b) for a, b, c in tr4))) == 4:
    p = sum(c for a, b, c in tr4)
    if p > perim:
        perim = p
        tr4s = tr4
print(perim, 'miles')
print('with triangles:', *tr4s)

1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

22

23

from itertools import combinations

sq = tuple(x * x for x in range(100))

trg = []

# find right triangles with integral sides

for i in range(10, 98):

for j in range(i + 1, 99):

k2 = i * i + j * j

if k2 in sq:

k = sq.index(k2)

trg. append((i, j, k))

perim = 0

for tr4 in combinations(trg, 4):

# can the four triangles form a quadrilateral

if len(set().union(*((a, b) for a, b, c in tr4))) == 4:

p = sum(c for a, b, c in tr4)

if p > perim:

perim = p

tr4s = tr4

print(perim, 'miles')

print('with triangles:', *tr4s)

Reply

BRG permalink

Hi John,

A good solution if we assume that the distances between the ranches and the point where the AS and DF lines cross are all integer. I am hoping this is the only solution but right now I am looking at solutions without making this assumption.

Reply

Frits permalink

Using a formula for the area of the quadrilateral.


# return altitude of triangle a, b, c with base <a>
def altitude(a, b, c):
  # Heron's formula 
  s = (a + b + c) / 2
  # calculate area
  a2 = (s * (s - a) * (s - b) * (s - c))
  if a2 < 0: return None
  area = a2**.5
  return 2 * area / a
    
# return divisor pairs of <n> between <m> and <M> (including)
def divisors_pairs(n, m, M, mx):
  # largest divisor >= mx
  i = min(n // mx, int(n**.5)) if mx else int(n**.5) 
  while i >= m:
    d, r = divmod(n, i)
    if d > M: break
    if not r:
      yield (i, d)
    i -= 1    

# return 2 diagonals of quadrilateral
def diagonals(area, m, M, mx):
  for p, q in divisors_pairs(2 * area, m, M, mx):
    if all(m <= x < M for x in [p, q]):
      if q >= mx: yield p
      if p >= mx: yield q

# are values a and b almost the same?
same = lambda a, b, eps=1e-10: abs(a - b) < eps   

#        A
#
#      D x      F
#
# 
#        S

sols, mx = set(), 0
# FS > AF and FS > DS
for AD in range(97, 9, -1):
  AD2 = AD * AD
  for DS in range(AD + 1, 99):
    DS2 = DS * DS
    # DS^2 - AD^2 = (FS + AF).(FS - AF) 
    for x, y in divisors_pairs(DS2 - AD2, 1, 99 + 98, 0):
      if x % 2 != y % 2: continue
      AF, FS = sorted([(y - x) // 2, (x + y) // 2])
      if AF <= AD or FS > 99: continue
      
      # calculate area of quadrilateral using outcome of Brahmagupta's formula
      area = ((AF + DS)**2 - (FS - AD)**2) // 4
      # 2 * area = DF * AS 
      for DF in set(diagonals(area, AD + 1, 99, mx)):
        fAx = altitude(DF, AD, AF)
        if fAx is None or fAx in {0, AD}: continue
        fSx = altitude(DF, DS, FS)
        if fSx is None: continue
        fAS = fAx + fSx
        if same(fAS, (AS := round(fAS))):
          if AS > 99: continue
          if AS >= mx:
            if AS > mx:
              sols = set()
              mx = AS
            sols.add(AD + AF + DS + FS)  

print("answer:", ' or '.join(str(x) for x in sols))

1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

25

26

27

28

29

30

31

32

33

34

35

36

37

38

39

40

41

42

43

44

45

46

47

48

49

50

51

52

53

54

55

56

57

58

59

60

61

62

63

64

65

66

67

68

69

70

# return altitude of triangle a, b, c with base <a>

def altitude(a, b, c):

# Heron's formula

s = (a + b + c) / 2

# calculate area

a2 = (s * (s - a) * (s - b) * (s - c))

if a2 < 0: return None

area = a2**.5

return 2 * area / a

# return divisor pairs of <n> between <m> and <M> (including)

def divisors_pairs(n, m, M, mx):

# largest divisor >= mx

i = min(n // mx, int(n**.5)) if mx else int(n**.5)

while i >= m:

d, r = divmod(n, i)

if d > M: break

if not r:

yield (i, d)

i -= 1

# return 2 diagonals of quadrilateral

def diagonals(area, m, M, mx):

for p, q in divisors_pairs(2 * area, m, M, mx):

if all(m <= x < M for x in [p, q]):

if q >= mx: yield p

if p >= mx: yield q

# are values a and b almost the same?

same = lambda a, b, eps=1e-10: abs(a - b) < eps

# A

#

# D x F

#

# S

sols, mx = set(), 0

# FS > AF and FS > DS

for AD in range(97, 9, -1):

AD2 = AD * AD

for DS in range(AD + 1, 99):

DS2 = DS * DS

# DS^2 - AD^2 = (FS + AF).(FS - AF)

for x, y in divisors_pairs(DS2 - AD2, 1, 99 + 98, 0):

if x % 2 != y % 2: continue

AF, FS = sorted([(y - x) // 2, (x + y) // 2])

if AF <= AD or FS > 99: continue

# calculate area of quadrilateral using outcome of Brahmagupta's formula

area = ((AF + DS)**2 - (FS - AD)**2) // 4

# 2 * area = DF * AS

for DF in set(diagonals(area, AD + 1, 99, mx)):

fAx = altitude(DF, AD, AF)

if fAx is None or fAx in {0, AD}: continue

fSx = altitude(DF, DS, FS)

if fSx is None: continue

fAS = fAx + fSx

if same(fAS, (AS := round(fAS))):

if AS > 99: continue

if AS >= mx:

if AS > mx:

sols = set()

mx = AS

sols.add(AD + AF + DS + FS)

print("answer:", ' or '.join(str(x) for x in sols))

Reply

Frits permalink

A Google AI Overview gave me the following info:

An orthodiagonal quadrilateral with integer diagonal lengths (and typically integer side lengths and area) is known as a Brahmagupta quadrilateral.

Apparently this is not correct as Brahmagupta quadrilaterals have to be cyclic as well.
So the program above doesn’t cover the full solution space.

Reply

BRG permalink

It is interesting that this issue has turned up here at the same time that Sundar Pichai (CEO of Alphabet, Google’s parent company) warned about over-reliance on AI because it is prone to making errors.

I am not surprised that it got this one wrong because there is a lot of confused internet coverage on the relationship between orthodiagonal and Brahmagupta quadrilaterals

Reply

John Z permalink

A different approach gives a different answer


#For any orthodiagonal quadrilateral, the sum of the squares of two opposite
#sides equals that of the other two opposite sides:
#for successive sides a, b, c, and d, we have

#    a^2 + c^2 = b^2 + d^2 .

from itertools import combinations

for a, b, d, c in combinations(range(99, 0, -1), 4):

  if a * a + c * c == b * b + d * d:
    print("perimeter is",  a + b + c + d, "with sides", a, b, d, c)
    break

1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

#For any orthodiagonal quadrilateral, the sum of the squares of two opposite

#sides equals that of the other two opposite sides:

#for successive sides a, b, c, and d, we have

# a^2 + c^2 = b^2 + d^2 .

from itertools import combinations

for a, b, d, c in combinations(range(99, 0, -1), 4):

if a * a + c * c == b * b + d * d:

print("perimeter is", a + b + c + d, "with sides", a, b, d, c)

break

Reply

BRG permalink

Hi John, But in this solution you are not applying the teaser constraint that the AS distance must be the largest possible.

Reply

Sunday Times Teaser 3295 – Always as the Crow Flies

by Andrew Skidmore

Published Sunday November 16 2025 (link)

Leave a comment to BRG Cancel reply

Recent Posts

Recent Comments

Archives

Categories

Meta