1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 |
sqrt3 = 3 ** (1/2) BL = [] # list of balls as potential answers for s in range(3, 50): # square pyramid for t in range(s, 50): # triangular pyramid # using Brian's derived condition for layer ratios if t == 2 * s: balls_s = 4 * (s * (s + 1) * (2 * s + 1)) // 6 balls_t = (t * (t + 1) * (t + 2)) // 6 # working in r^2 units (omitted from As and At below) As = 16 * s**2 At = sqrt3 * (t - 1 + sqrt3) ** 2 # find the minimum number of balls # ..for a more than 55% base area reduction if 100 * At < 45 * As: BL.append(balls_s) # show minimum number of balls for display print(f"Edward displayed {min(BL)} balls.") |

1 2 3 |
if c == n[:len(c)]: |

Walrus notation could also be used to good effect in some places to avoid recalculations.

]]>
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 |
# balls in four S layer square base pyramids: 4.S.(S + 1).(2.S + 1) / 6 # balls in one T layer triangular base pyramid: T.(T + 1).(T + 2) / 6 # # The difference of ball numbers in the two constructions: # # == [4.S.(S + 1).(2.S + 1) - T.(T + 1).(T + 2))] / 6 # == (2.S - T).(4.S^2 + 6.S + T^2 + 3.T + 2.S.T) / 6 # # Hence the two constructions have equal numbers of balls when T = 2.S # # with ball radius r: # # Area of square base = 4.(2.r.S)^2 = 16.(r.S)^2 # Area of triangle base = {sqrt(3)/4}.{2.r.(T - 1 + sqrt(3))}^2 sqrt3, s = 3 ** (1 / 2), 0 while True: s += 1 # the number of balls nb = 4 * s * (s + 1) * (2 * s + 1) // 6 # the areas of the square and triangular bases (in units of ball radii squared) bs, bt = 16 * s * s, sqrt3 * (2 * s - 1 + sqrt3) ** 2 # find the minimum number of balls for a more than 55% base area reduction if bt < 0.45 * bs: print(f"{nb} balls (4 x {s} layer squares, 1 x {2 * s} layer triangular).") break |

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 |
# build list of four digit primes p1 = (2, 3, 5, 7) p2 = p1 + tuple(x for x in range(11, 100, 2) if all(x % p for p in p1)) p4 = p2 + tuple(x for x in range(101, 10000, 2) if all(x % p for p in p2)) # four digit primes in < order as strings p4s = tuple(ps for p in p4[::-1] if len(ps := str(p)) == len(set(ps))) # cubes with four digits or less and no repeat digits cs = tuple(str(si) for i in range(1, 22) if len(si := str(i ** 3)) == len(set(si))) # make ten digit number from primes in < order with no repeat digits def solvep(i = 0, prs = '', ixs = ()): if len(prs) == 10: yield ixs else: for j, p in enumerate(p4s[i:]): if len(p) <= 10 - len(prs): if len(set(prs + p)) == len(p) + len(prs): yield from solvep(i + j + 1, prs + p, ixs + (i + j,)) # decompose number into cubes of 4 digits or less def solvec(n, ixs = ()): if len(n) == 0: yield(ixs) else: for i, c in enumerate(cs): if len(c) <= len(n): if all(x[0] == x[1] for x in zip(c, n)): yield from solvec(n[len(c):], ixs + (i,)) for n in solvep(): # make large number from all primes but the smallest prs = ''.join((p4s[i] for i in(n[:-1]))) # decompose into cubes for x in solvec(prs): print('primes:',*[p4s[i] for i in(n[:-1])]) print('discarded prime:', p4s[n[-1]]) print('cubes:',*[cs[i] for i in x]) |

The diagrams above illustrate that if the disc is to fall entirely within a single hexagonal (triangular) tile, it centre has to fall within the inner green areas shown. Hence the probability of this happening is simply the ratio of the areas of these inner and outer shapes.

If the radius of the circular disc is \(r\) and the sides of the regular hexagon and equilateral triangle are \(h\) and \(t\) respectively, inspection of the small lower right triangles in each shape shows that the side lengths of the inner hexagon and the triangle are respectively \(h\,-\, 2r/\sqrt{3}\) and \(t\,-\,2\sqrt{3}r\). With \(p_h\) and \(p_t\) as the probabilities of the disc falling entirely within a tile, these become:\[p_h = \left(\frac{h \;- 2r/\sqrt{3}}{h}\right)^2=\left(1-\frac{2r}{\sqrt{3}h}\right)^2\]

\[p_t = \left(\frac{t\; – 2\sqrt{3}r}{t}\right)^2=\left(1-\frac{2\sqrt{3}r}{t}\right)^2\]

Equating these shows that equal probabilities requires that \(t=3h\). (see [1])

Since the sides are triangular numbers, let the hexagon’s and the triangle’s numbers be \(h(h+1)/2\) and \(t(t+1)/2\) respectively (this is a a change of notation from that above). We then have: \[t(t+1) = 3h(h+1)\] which, after rearrangement, can be put in the form:\[(2t+1)^2\; -\; 3(2h+1)^2=-2\]This is a quadratic Diophantine equation for which there are standard techniques for solution which can be used to show that, if we can find one solution, an infinite series of solutions can be produced using two recursions:\[t_{n+1}=2 t_{n}+3h_{n}+2\]\[h_{n+1}=t_{n}+2h_{n}+1\] Since the equation has the solution \((t, h) = (0, 0)\) the sequence of solutions can easily be generated \[(t,h) = (0,0),(2,1),(9,5),(35,20),(132,76), …\] from which the possible side lengths are \[(0,0),(3,1),(45,15),(630,210),(8778,2926), …\]The solution is hence equilateral triangle tiles with sides of 630 mm and regular hexagon tiles with sides of 210 mm.

[1] More generally for regular polyhedra with \(m\) and \(n\) sides \((s)\) equal probabilities requires that \[\tan(180/m)/s_m = \tan(180/n)/s_n\]Here is a programmed solution using the quadratic Diophantine equation solver in my number theory library (available here).

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 |
from number_theory import Qde # triangle_side = 3 x hexagon_side # # triangle side: t.(t + 1) / 2 # hexagon side: h.(h + 1) / 2 # # t.(t + 1) / 2 = 3.h.(h + 1) / 2 ==> # # (2.t + 1)^2 - 3.(2.h + 1)^2 == -2 # Hence solve the quadratic diophantine equation # generate solutions for 2.t + 1 and 2.h + 1 t = Qde(3, -2) for ttp1, hhp1 in t.solve(limit=4): if ttp1 % 2 and hhp1 % 2: t, h = ttp1 // 2, hhp1 // 2 t_side, h_side = t * (t + 1) // 2, h * (h + 1) // 2 # both sides must be even (and non zero) if t and h and not (t_side % 2 or h_side % 2): print(f"Triangle side = {t_side} mm, Hexagon side = {h_side} mm") |

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 71 72 73 74 75 |
from collections import defaultdict from itertools import product from number_theory import Primes # test for no duplicate digits diff_dgts = lambda n: len((s := str(n))) == len(set(s)) # cubes with no repeated digits up to four digits in length cubes = set(c for n in range(22) if diff_dgts((c := n ** 3))) # decompose a string of digits in <s> into integer values # in setf (that have four or less digits) def decompose(s, setf, p=()): if s == '': yield p else: for i in range(1, min(4, len(s)) + 1): v = int(s[:i]) if s[0] != '0' and v in setf: yield from decompose(s[i:], setf, p + (v,)) # collect primes indexed on their length and # then on their sorted digits n2d = defaultdict(lambda: defaultdict(list)) for p in Primes(10000): str_p = str(p) set_p = set(str_p) if len(str_p) == len(set_p): dgts = ''.join(sorted(str_p)) n2d[len(dgts)][dgts].append(p) # combine digit sequences from <n2d> with lengths # in <ns> and return those combinations of disjoint # sequences that contain all ten digits def compose(ns, n2d, setd=set(), dseq=()): if not ns: if len(setd) == 10: yield dseq else: for ds in n2d[ns[0]]: if not setd.intersection(ds) and ds not in dseq: yield from compose(ns[1:], n2d, setd.union(ds), dseq + (ds, )) # partition the integer <n> into a sequence of # non-increasing integers of maximum length <mx> def partition(n, mx, seq=()): if n == 0: yield seq else: for i in range(1, min(mx, n) + 1): if not seq or i <= seq[-1]: yield from partition(n - i, mx, seq + (i, )) # consider all lengths (in decreasing order) for primes with # lengths that sum to ten for s in partition(10, 4): # consider all digit sequences that have all ten digits for dgs in compose(s, n2d): # now convert the digit sequences into the primes that # can be made from them for pd in product(*(n2d[len(d)][d] for d in dgs)): # check that the primes are in decreasing order if pd == tuple(sorted(pd, reverse=True)): # concatenate the digits of the primes in order (and # without the smallest) to create a single integer pcbs = ''.join(str(p) for p in pd[:-1]) # find resulting integers that can be split # into digit sequences that are all cubes for q in decompose(pcbs, cubes): print(f"{pd} -> {pcbs} -> {q}") |

This runs in about 300ms. But tackling the cubes first is much faster, running in less than 30ms on my PC:

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 |
from collections import defaultdict from number_theory import Primes # test for no duplicate digits diff_dgts = lambda n: len((s := str(n))) == len(set(s)) # primes and cubes with up to four digits primes = set(p for p in Primes(10000) if diff_dgts(p)) cubes = set(c for n in range(22) if diff_dgts((c := n ** 3))) # collect primes indexed on their length and # then on their sorted digits n2d = defaultdict(lambda: defaultdict(list)) for p in Primes(10000): str_p = str(p) set_p = set(str_p) if len(str_p) == len(set_p): dgts = ''.join(sorted(str_p)) n2d[len(dgts)][dgts].append(p) # generate sequences of cubes from <cubes> with disjoint digits for # which the unused digits are capable of making single primes; if so # return the cube sequence and the list of makeable single primes def cube_seq(n2d, cubes, dgts=set(), seq=()): if seq and dgts: # sort unused digits so far and check if they can form a single prime k = ''.join(sorted(set('0123456789').difference(''.join(sorted(dgts))))) # if they can return the cube sequence and the list of primes that can # be made from the unused digits if (t := n2d[len(k)][k]): yield seq, t # add another cube to the sequence for c in cubes: # if it does not use previously used digits if not dgts.intersection(str(c)): yield from cube_seq(n2d, cubes.difference(str(c)), dgts.union(str(c)), seq + (c,)) # split a sequence of primes from the digit sequence <cseq> by # splitting it into consecutive sub-sequences that form primes # of between one and four digits in decreasing order def prime_seq(cseq, n2d, pseq=()): if not cseq: yield pseq else: # consider the length of the next prime for i in range(1, min(4, len(cseq)) + 1): # avoid a leading zero if cseq[0] == '0': continue p = int(cseq[:i]) # is it a prime in a decreasing size order if p in primes and (not pseq or p < pseq[-1]): yield from prime_seq(cseq[i:], n2d, pseq + (p, )) # consider possible sequences of cubes and last primes for cseq, t in cube_seq(n2d, cubes): # form the digit sequence from the digits of the cubes ss = ''.join(str(x) for x in cseq) # can this sequence make a decreasing prime sequence for ps in prime_seq(ss, n2d): pr = tuple(p for p in t if p < ps[-1]) if pr: print(f"{ps + pr} -> {ss} -> {cseq}") |

Forwarned that there is no solution for ‘7/96’, a solution is considered found if the denominator of the probability = 96.

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 71 72 73 74 75 76 |
from itertools import product, combinations from fractions import Fraction from math import gcd # Original couples are A B C D, new couple in game 2 is E; in game 3 is F cp = set('ABCD') # win or lose wl = 'WL' # number of times one of 'ABCD' wins second game; no win in first game: £2K pLWX = 0 # times one of 'ABCD' wins third game after a win in game 1, no win in 2: £2K pWLW = 0 # times one of 'ABCD' wins third game after no wins in game 1 and 2: £3K pLLW = 0 # consider possible first game outcomes: couple in fourth round; W/L for c1, g1 in product(cp, wl): # replace couple who was in fourth round of game 1 cp2 = cp - {c1} | {'E'} # consider possible second game outcomes: couple in fourth round; W/L for c2, g2 in product(cp2, wl): if g1 == "L" and g2 == "W": if c2 in cp: # c2 has won £2000 pLWX += 1 elif g1 == g2 == "W": # no possibility of £2K or £3K win, only £1K wins continue else: # go on to third game if first 2 are WL or LL # replace couple who was in fourth round of game 2 cp3 = cp2 - {c2} | {'F'} # consider possible third game outcomes: couple in fourth round; # W/L is indifferent c3 wins half, loses half for c3 in cp3: if g1 == "W" and g2 == "L" and c3 in cp: # c3 wins £2000 pWLW += 1 elif g1 == g2 == "L" and c3 in cp: # c3 wins £3000 pLLW += 1 print(f'Chances of £2000 win in games 2 & 3: {pLWX}, {pWLW}; £3000: {pLLW}') # Second stage: try different values of p # consider numerator and denominator of p for n, d in combinations(range(1, 100), 2): # p must be < 1/2 if n * 2 >= d or gcd(n, d) > 1: continue # prob of winning in the fourth round of a game p = Fraction(n, d) # calculate probabilty of one of the starting couples winning £2000 prob2 = (Fraction(pLWX * (1 - p) * p, 64) + Fraction(pWLW * p * (1 - p) * p, 512)) # or £2000 or £3000 prob23 = prob2 + Fraction(pLLW * p * (1 - p) * p, 512) # found solution with correct denominator? if prob2.denominator == 96 or prob23.denominator == 96: print(f'p = {p}, probability of winning £2000: {prob2}') print(f' probability of winning £2K or £3K: {prob23}') |

Consequently I don’t have a full Python solution. But here is a program that paves the way to a discussion of possible solutions by giving the equations for the probability of winning the various possible amounts in terms of the probability \(p\) described in the teaser.

Please note that this code uses the latest version of my polynomial library which I have just updated to provide features that help with this teaser.

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 |
from polynomial import Poly from fractions import Fraction as RF from gmpy2 import mpq as RF # probability of winning the jackpot in the final round (= p) win = Poly([0, 1]) # credit: Jim Randell # generate the polynomials P[p] that determine the probabilities of # all possible outcome amounts (in units of the base prize) with a # maximum of <n> games; return: <amount>, <probability polynomial> def outcomes(n, t=1, cp=1): if n == 0: # finished without a win yield 0, cp else: # there is a 1/4 chance of reaching the final and leaving q = cp * RF(1, 4) # with a win yield t, q * win # or a loss yield 0, q * (1 - win) # there is a 3/4 chance of getting knocked out # and then entering more gaames after q = cp * RF(3, 4) # a win and no rollover yield from outcomes(n - 1, 1, q * win) # or a loss and a rollover yield from outcomes(n - 1, t + 1, q * (1 - win)) # accumulate probabilities of winning amounts over three games cp, ps = dict(), Poly([0]) for v, p in outcomes(3): ps += p if v not in cp: cp[v] = Poly([0]) cp[v] += p # output the amounts and their probability polynomials for v, p in sorted(cp.items()): print(f" £{v}000 ==> P[p] = {p:p}") print(f" Sum ==> P[p] = {ps}") cp[23] = cp[2] + cp[3] cp[123] = cp[1] + cp[23] print(f">= £1000 ==> P[p] = {cp[123]:p}") print(f">= £2000 ==> P[p] = {cp[23]:p}") print("\nPossible Solutions:") # consider solutions for other amounts and other # probabilities that produce rational results for c, f in ((cp[2], ' £2000'), (cp[3], ' £3000'), (cp[23], '>= £2000')): for n in range(-20, 21): # for the polynomial to solve r = (c - RF(n, 96)).zq_roots()[1] if r: r = r[0][0] if 0 <= r <= 1: print(f" {f} (p = {n}/96) -> {r}") |

Here is its output (edited):

1 2 3 4 5 6 7 8 9 |
£0000 ==> P[p] = -(37/64)p + 1 £1000 ==> P[p] = (21/64)p^2 + (1/4)p £2000 ==> P[p] = -(9/64)p^3 - (3/64)p^2 + (3/16)p £3000 ==> P[p] = (9/64)p^3 - (9/32)p^2 + (9/64)p Sum ==> P[p] = 1 >= £1000 ==> P[p] = (37/64)p >= £2000 ==> P[p] = -(21/64)p^2 + (21/64)p |

If we set \(P[p]\) to be \(q/96\) for a win of £2000 and simplify the appropriate polynomial, we obtain the cubic equation \[27p^3+9p^2-36p+2q=0\]

For the teaser \(q\) is 7, for which the polynomial has one negative real root and two complex ones, none of which can solve the teaser. In practice it seems likely that the probability \(p\) set by the author in the teaser is a rational fraction so it makes sense to consider other \(q\) values to see if any provide solutions of this form. When this is done two possible solutions emerge: \(p=1/3\) when \(q = 5\) and \(p=2/3\) when \(q=6\). The teaser hence has solution when the teaser’s probability of \(7/96\) is replaced by \(5/96\).

Jim Randell (on his S2T2 site) has suggested that “£2000” could be interpreted as “£2000 or more”. When we use the polynomial for this value and set it to \(7/96\) it simplifies to give a quadratic equation:\[9p^2-9p+2=0\]with solutions of \(p = 1/3\) and \(p = 2/3\) the first of which meets the teaser constraints.

Of these I prefer the former since it is nice to see a solution involving a cubic instead of the much more common quadratics.

]]>
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 |
from math import isqrt # b = fraction of notes that are forgeries (bad); f = fraction # of notes wrongly classified; t = number of tests run # # 1. genuine notes wrongly tested as forgeries = f.(1 - b).t # 2. genuine notes correctly tested as genuine = (1 - f)(1 - b).t # 3. forged notes wrongly tested as genuine = f.b.t # 4. forged notes correctly tested as forged = (1 - f).b.t # # Test error rate (1 + 3): b.t = t / N ==> f = 1 / N # # Of tests indicating forgeries (1 + 4) only 1 was correct: # f.(1- b).t = N - 1 ; (1 - f).b.t = 1 # ==> f.(1 - b) = (N - 1)(1 - f).b # # Eliminating f and simpifying: 1 - b = (N - 1)^2.b # # ==> N = sqrt((1 - b) / b) + 1 # ==> N = sqrt(100 / d - 1) + 1 (with b = d / 100) # for d in range(1, 10): if 100 % d == 0: t = 100 // d - 1 if (r := isqrt(t)) ** 2 == t: print(f"N = {r + 1} ({d = })") |

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 |
from math import isqrt # the set of unit digits of squares with odd penultimate digits units = set(i ** 2 % 10 for i in range(10) if (i ** 2 // 10) & 1) # consider the possible lengths of the phone number for sq_len in units: # digits at positions 0, 1, 2, -2, -1 all non zero and different # plus one zero digit gives a minimum of six digits for r in range(317, isqrt(10 ** sq_len - 1) + 1): nbr = r ** 2 # its penultimat digit is odd if (nbr // 10) & 1: # the digits of the number d = [int(d) for d in str(nbr)] # its digits are all different with one zero ddigit if len(set(d)) == sq_len and d.count(0) == 1: # its second, third and last but two digits are exact # proper multiples of its first digit if all(n > d[0] and n % d[0] == 0 for n in (d[1], d[2], d[-2])): print(f"She dialled {nbr}.") |