# Mathematician vs. Computer: A Game

A mathematician and a computer are playing a game: First, the mathematician chooses an integer from the range $2,...,1000$. Then, the computer chooses an integer uniformly at random from the same range. If the numbers chosen share a prime factor, the larger number wins. If they do not, the smaller number wins. (If the two numbers are the same, the game is a draw.)

Which number should the mathematician choose in order to maximize his chances of winning?

• Not that it will matter, but the rules of the game should specify what happens if the numbers match. Explore the probability the mathematician wins if she chooses various primes. – André Nicolas Oct 26 '14 at 19:43
• If the numbers match, the game is a draw. – user139000 Oct 26 '14 at 19:44
• Comments are not for extended discussion; this conversation has been moved to chat. – user642796 Oct 28 '14 at 5:01
• Since the top answers, while enlightening, have become a bit TLDR, with the actual solution hard to find, here are the current results: The optimal choices for the given range are 31 and 29, and for the generalized problem with a range up to $N$, it has been shown that the closest prime to $\sqrt{N}$ is optimal if $N$ is sufficiently large. – user139000 Oct 29 '14 at 14:24

## 11 Answers

First consider choosing a prime $p$ in the range $[2,N]$. You lose only if the computer chooses a multiple of $p$ or a number smaller than $p$, which occurs with probability $$\frac{(\lfloor{N/p}\rfloor-1)+(p-2)}{N-1}=\frac{\lfloor{p+N/p}\rfloor-3}{N-1}.$$ The term inside the floor function has derivative $$1-\frac{N}{p^2},$$ so it increases for $p\le \sqrt{N}$ and decreases thereafter. The floor function does not change this behavior. So the best prime to choose is always one of the two closest primes to $\sqrt{N}$ (the one on its left and one its right, unless $N$ is the square of a prime). Your chance of losing with this strategy will be $\sim 2/\sqrt{N}$.

On the other hand, consider choosing a composite $q$ whose prime factors are $$p_1 \le p_2 \le \ldots \le p_k.$$ Then the computer certainly wins if it chooses a prime number less than $q$ (other than any of the $p$'s); there are about $q / \log q$ of these by the prime number theorem. It also wins if it chooses a multiple of $p_1$ larger than $q$; there are about $(N-q)/p_1$ of these. Since $p_1 \le \sqrt{q}$ (because $q$ is composite), the computer's chance of winning here is at least about $$\frac{q}{N\log q}+\frac{N-q}{N\sqrt{q}}.$$ The first term increases with $q$, and the second term decreases. The second term is larger than $(1/3)/\sqrt{N}$ until $q \ge (19-\sqrt{37})N/18 \approx 0.72 N$, at which point the first is already $0.72 / \log{N}$, which is itself larger than $(5/3)/\sqrt{N}$ as long as $N > 124$. So the sum of these terms will always be larger than $2/\sqrt{N}$ for $N > 124$ or so, meaning that the computer has a better chance of winning than if you'd chosen the best prime.

This rough calculation shows that choosing the prime closest to $\sqrt{N}$ is the best strategy for sufficiently large $N$, where "sufficiently large" means larger than about $100$. (Other answers have listed the exceptions, the largest of which appears to be $N=30$, consistent with this calculation.)

• Although I have upvoted several of the other answers as well, I have finally chosen to accept this one because it was the first answer to give a full theoretical explanation for what is going on. – user139000 Nov 1 '14 at 6:45

For fixed range:

range = 16;
a = Table[Table[FactorInteger[y][[n, 1]], {n, 1, PrimeNu[y]}], {y, 1, range}];
b = Table[Sort@DeleteDuplicates@ Flatten@Table[
Table[Position[a, a[[y, m]]][[n, 1]],
{n, 1, Length@Position[a, a[[y, m]]]}], {m, 1, PrimeNu[y]}], {y, 1, range}];
c = Table[Complement[Range[range], b[[n]]], {n, 1, range}];
d = Table[Range[n, range], {n, 1, range}];
e = Table[Range[1, n], {n, 1, range}];
w = Table[DeleteCases[DeleteCases[Join[Intersection[c[[n]], e[[n]]],
Intersection[b[[n]], d[[n]]]], 1], n], {n, 1, range}];
l = Table[DeleteCases[DeleteCases[Complement[Range[range], w[[n]]], 1],
n], {n, 1, range}];
results = Table[Length@l[[n]], {n, 1, range}];
cf = Grid[{{Join[{"n"}, Rest@(r = Range[range])] // ColumnForm,
Join[{"win against n"}, Rest@w] // ColumnForm,
Join[{"lose against n"}, Rest@l] // ColumnForm,
Join[{"probability win for n"}, (p = Drop[Table[
results[[n]]/Total@Drop[results, 1] // N,{n, 1, range}], 1])] // ColumnForm}}]
Flatten[Position[p, Max@p] + 1]


isn't great code, but fun to play with for small ranges, gives

and perhaps more illuminating

rr = 20; Grid[{{Join[{"range"}, Rest@(r = Range[rr])] // ColumnForm,
Join[{"best n"}, (t = Rest@Table[
a = Table[Table[FactorInteger[y][[n, 1]], {n, 1, PrimeNu[y]}], {y, 1, range}];
b = Table[Sort@DeleteDuplicates@Flatten@Table[Table[
Position[a, a[[y, m]]][[n, 1]], {n, 1,Length@Position[a, a[[y, m]]]}],
{m, 1,PrimeNu[y]}], {y, 1, range}];
c = Table[Complement[Range[range], b[[n]]], {n, 1, range}];
d = Table[Range[n, range], {n, 1, range}];
e = Table[Range[1, n], {n, 1, range}];
w = Table[DeleteCases[DeleteCases[Join[Intersection[c[[n]], e[[n]]],
Intersection[b[[n]], d[[n]]]], 1], n], {n, 1, range}];
l = Table[DeleteCases[DeleteCases[Complement[Range[range], w[[n]]], 1], n],
{n,1, range}];
results = Table[Length@l[[n]], {n, 1, range}];
p = Drop[Table[results[[n]]/Total@Drop[results, 1] // N,
{n, 1, range}], 1];
{Flatten[Position[p, Max@p] + 1], Max@p}, {range, 1, rr}]/.Indeterminate-> draw);
Table[t[[n, 1]], {n, 1, rr - 1}]] // ColumnForm,
Join[{"probability for win"}, Table[t[[n, 2]], {n, 1, rr - 1}]] // ColumnForm}}]


compares ranges:

Plotting mean "best $n$" against $\sqrt{\text{range}}$ gives

For range=$1000,$ "best $n$" are $29$ and $31$, which can be seen as maxima in this plot:

# Update

In light of DanielV's comment that a "primes vs winchance" graph would probably be enlightening, I did a little bit of digging, and it turns out that it is. Looking at the "winchance" (just a weighting for $n$) of the primes in the range only, it is possible to give a fairly accurate prediction using

range = 1000;
a = Table[Table[FactorInteger[y][[n, 1]], {n, 1, PrimeNu[y]}], {y, 1, range}];
b = Table[Sort@DeleteDuplicates@Flatten@Table[
Table[Position[a, a[[y, m]]][[n, 1]], {n, 1,
Length@Position[a, a[[y, m]]]}], {m, 1, PrimeNu[y]}], {y, 1, range}];
c = Table[Complement[Range[range], b[[n]]], {n, 1, range}];
d = Table[Range[n, range], {n, 1, range}];
e = Table[Range[1, n], {n, 1, range}];
w = Table[    DeleteCases[    DeleteCases[
Join[Intersection[c[[n]], e[[n]]], Intersection[b[[n]], d[[n]]]],
1], n], {n, 1, range}];
l = Table[
DeleteCases[DeleteCases[Complement[Range[range], w[[n]]], 1],
n], {n, 1, range}];
results = Table[Length@l[[n]], {n, 1, range}];
p = Drop[Table[
results[[n]]/Total@Drop[results, 1] // N, {n, 1, range}], 1];
{Flatten[Position[p, Max@p] + 1], Max@p};
qq = Prime[Range[PrimePi[2], PrimePi[range]]] - 1;
Show[ListLinePlot[Table[p[[t]] range, {t, qq}],
DataRange -> {1, Length@qq}],
ListLinePlot[
Table[2 - 2/Prime[x] - 2/range (-E + Prime[x]), {x, 1, Length@qq + 0}],
PlotStyle -> Red], PlotRange -> All]


The plot above (there are $2$ plots here) show the values of "winchance" for primes against a plot of $$2+\frac{2 (e-p_n)}{\text{range}}-\frac{2}{p_n}$$

where $p_n$ is the $n$th prime, and "winchance" is the number of possible wins for $n$ divided by total number of possible wins in range ie $$\dfrac{\text{range}}{2}\left(\text{range}-1\right)$$ eg $499500$ for range $1000$.

Show[p // ListLinePlot,  ListPlot[N[
Transpose@{Prime[Range[PrimePi[2] PrimePi[range]]],
Table[(2 + (2*(E - Prime[x]))/range - 2/Prime[x])/range, {x, 1,
Length@qq}]}], PlotStyle -> {Thick, Red, PointSize[Medium]},
DataRange -> {1, range}]]


# Added

Bit of fun with game simulation:

games = 100; range = 30;
table = Prime[Range[PrimePi[range]]];
choice = Nearest[table, Round[Sqrt[range]]][[1]];
y = RandomChoice[Range[2, range], games];  z = Table[
Table[FactorInteger[y[[m]]][[n, 1]], {n, 1, PrimeNu[y[[m]]]}], {m, 1, games}];
Count[Table[If[Count[z, choice] == 0 && y[[m]] < choice \[Or]
Count[z, choice] > 0 && y[[m]] < choice, "lose", "win"],
{m, 1, games}], "win"]


& simulated wins against computer over variety of ranges

with

Clear[range]
highestRange = 1000;
ListLinePlot[Table[games = 100;
table = Prime[Range[PrimePi[range]]];
choice = Nearest[table, Round[Sqrt[range]]][[1]];
y = RandomChoice[Range[2, range], games];
z = Table[Table[FactorInteger[y[[m]]][[n, 1]], {n, 1, PrimeNu[y[[m]]]}], {m,
1, games}];
Count[Table[ If[Count[z, choice] == 0 && y[[m]] < choice \[Or]
Count[z, choice] > 0 && y[[m]] < choice, "lose", "win"], {m, 1,
games}], "win"], {range,2, highestRange}], Filling -> Axis, PlotRange-> All]


# Added 2

Plot of mean "best $n$" up to range$=1000$ with tentative conjectured error bound of $\pm\dfrac{\sqrt{\text{range}}}{\log(\text{range})}$ for range$>30$.

I could well be wrong here though. - In fact, on reflection, I think I am (related).

• Could you please answer the original question: what is the best play for $n=1000$? I can't find it. – Ross Millikan Oct 27 '14 at 2:24
• @RossMillikan apologies - now added – martin Oct 27 '14 at 2:55
• In the first plot, the huge spike at the highly composite number 30 is interesting. Why are there no more spikes after it? – user139000 Oct 27 '14 at 4:42
• Another very interesting result from the first plot is that the "best $n$" function is not monotonous, so a simple description such as "next smaller prime to $\sqrt{range}$" is not sufficient. – user139000 Oct 27 '14 at 4:45
• A "primes vs winchance" graph would probably be enlightening. – DanielV Oct 27 '14 at 7:39

Scheme:

(define (range a b) (if (> a b) '() (cons a (range (+ a 1) b))))

;; Probability of Winning if you Choose N out of NRange
(define (ProbOfWinning N NRange)
;; how many numbers k for which
;;   GCD(k, N)=1 and N < k
;; or
;;   GCD(k, N)>1 and N > k
;; N beats K
(define (Wins k)
(if (= (gcd k N) 1) (< N k) (> N k)))

(/ (length (filter Wins NRange)) (length NRange)))

(define GameRange (range 2 1000))

(define WinChances (map (λ (N) (list N (ProbOfWinning N GameRange))) GameRange))

(sort WinChances (λ (a b) (< (second a) (second b))) )


The best chances of winning are $31$ and $29$, with probability $938/999$, and close third $37$ at $937/999$

• I love it that we have two computer-generated "proofs" with different results :) – user139000 Oct 26 '14 at 20:24
• Haha yeah I notice, double checking now. This btw isn't a proof, just a computation. – DanielV Oct 26 '14 at 20:25
• If the computation is correct it can be turned into a proof by just printing out the probability for every number in the range. As a side note, Lisp never ceases to surprise me in being simultaneously elegant and incomprehensible... – user139000 Oct 26 '14 at 20:30
• Let's let the largest number that can be picked be $n$ (so in the problem $n = 1000$.) If the mathematician picks a prime $p$, then they lose if a multiple of $p$ or a number less than $p$ is picked. There are $n/p + p$ such numbers (approximately), and this is minimized when $p = \sqrt{n}$. – Michael Lugo Oct 26 '14 at 20:57
• @MichaelLugo: If you now add a proof that that prime cannot be beaten by a large composite number, you have an answer that I'll gladly accept. – user139000 Oct 27 '14 at 4:48

I used Javascript and brute-forced the solution. It wasn't terribly fast, but it worked. 29 or 31 are the best numbers to pick.

function find_primes(max) {
// returns all primes less than or equal to max
for (var primes = [], i = 2; i <= max; ++ i) {
for (var j = 0; j < primes.length; ++ j)
if (i % primes[j] == 0) break;  // stop if it's divisible by a prime
if (j >= primes.length) primes.push(i);
}
return primes;
}

function winner(a, b, primes) {
// primes must contain all primes less than or equal to the largest a or b
for (var j = 0; j < primes.length; ++ j)
if (!(a % primes[j] + b % primes[j])) return Math.max(a, b);
return Math.min(a, b);
}

var min = 2, max = 1000, primes = find_primes(max), plays = [];
for (var i = min; i <= max; ++ i) {
plays[i] = {play: i, opponents: 0, wins: 0};
for (var j = min; j <= max; ++ j) {
++ plays[i].opponents;
if (i != j && winner(i, j, primes) == i) ++ plays[i].wins;
}
}

// sort, highest number of wins first
plays.sort(function (a, b) { return a.wins < b.wins; });
// display the results
plays.map(function (e) {
return e.play + ": " + e.wins + " wins (" + (e.wins / e.opponents * 100).toFixed(2) * 1 + "%)";
});


(edit) you also get an interesting pattern if you print the entire list with prime numbers highlighted, like this:

• Your code seems to run to $2000$ while the problem specifies an upper bound of $1000$ for the range. – user139000 Oct 27 '14 at 15:56
• blink you're absolutely right... I coulda sworn that it said the upper bound was 2000. Ok, I changed my max variable to 1000 and the best choices are 29 and 31, same as everyone else got. – guest Oct 27 '14 at 16:04

I feel like all of the other responses are somewhat more complicated and hard to follow, so I hope to provide something more straight-forward. I'm always impressed with the beautiful formatting on MSE, though.

First, any number $n$ you choose will have some set of prime factors $(p1, p2, ...)$. For each prime factor $p_i$, there are $x_i =\lfloor \frac{1000}{p_i} \rfloor$ numbers to consider which will have that common factor. If $n$ is chosen to be a prime itself, there are $\lfloor \frac{1000}{n} \rfloor - 1$ numbers greater than $n$ which will beat it, and it loses to any numbers less than itself, for a total of $$f(n)=n + \lfloor \frac{1000}{n} \rfloor - 2$$

Now if $n$ is composite, it loses to all of those numbers greater than it which share common factors

$$g(n)=\sum{\lfloor{\frac{1000}{p_i}}\rfloor} - \lfloor\frac{n}{p_i}\rfloor$$

as well as $\phi(n)$ numbers less than n which don't share common factors, with $\phi$ being Euler's totient function.

Clearly large primes fail, since $\phi(n) =n-1$ when n is prime. Likewise, small composite numbers fail, since $g(n)$ dominates there.

If we examine the formula for primes, we can see that it is minimized roughly when $n^2=1000$ by taking a derivative and solving. This curve monotonically decreases as it is dominated by $\lfloor\frac{1000}{n}\rfloor$ and monotonically increases as it is dominated by $n$. Thus only a few primes need to be tried to discover that $n=29$ and $n=31$ give the best prime solutions.

If any composite solution is going to out-perform a prime one, it must have a totient smaller than $f(31)=61$, at the very least, and as a result it cannot exceed $210$. Thus it must also have less than $61$ numbers between $210$ and $1000$ which it shares a factor with. $\lfloor(1000-210) / 61 \rfloor= 12$, so none of the factors can be less than $12$, but if there is more than one such factor, we would have a totient too large. Therefore there can be no composite solution better than the prime one.

Intuitively, you are just comparing the quadratic growth of the small prime choice to the $n/\log{n}$ growth of the large composite choice

Let's assume that you stick with picking prime numbers first. In this case, to win you should choose $i$ to maximize the probability $$p_{\text{win}}(i)\approx p_{\text{common}}(i)p_{\text{larger}}(i)+p_{\text{no common}}(i)p_{\text{smaller}}(i)\approx\frac{1}{i}\frac{i}{N}+\frac{i-1}{i}\frac{N-i}{N}$$ where $N$ is the largest available number for picking. The chance to win with the larger number is $1/N$, independent of your choice of $i$, so you can focus on maximizing the probability for a getting a smaller number with no common prime factors. The derivative of the second term with respect to $i$ $$\frac{\text{d}}{\text{d}i}\left(\frac{i-1}{i}\frac{N-i}{N}\right)=\frac{N}{i^2}-1,$$ has a zero at $i=\sqrt{N}$, which turns out to be a maximum of the winning probability (31 is the closest prime to $\sqrt{1000}$, which is also, among others, the result from the simulations in other answers).

The question remains why you should go for prime numbers in the first place and not for large numbers with a lot of prime factors, increasing the chance of having one in common with the computer's number while also being bigger. For that, I do not have a definitive answer, however, according to this: http://mathworld.wolfram.com/DistinctPrimeFactors.html, the number of distinct prime factors in a natural number increases with $ln(ln(N))$, much slower than the number of prime numbers up to that same number ($\pi(N)=N/ln(N)$). This seems to suggest that for larger $N$, picking large numbers with many prime factors makes less and less sense.

• The first $=$ should be $\approx$ too: "common" and "larger" are not independent. – user14972 Oct 27 '14 at 15:23
• So what is the last $N$ for which the best choice is highly composite? (This may need an estimate about the prime gap neer $\sqrt N$) – Hagen von Eitzen Oct 27 '14 at 15:26
• @HagenvonEitzen: Based on the chart in martin's answer, I would guess it's 30 (for which, interestingly, the best choice is 30 as well). – user139000 Oct 27 '14 at 15:37
• @Hurkyl: you are right, changed it – Roberto Oct 27 '14 at 16:16
• @pew: Isn't 30 the best choice for 31 as well? It ties with 5 and 7 for 32, but for 33 and above, the best choice seems to be prime. – Peter Shor Oct 27 '14 at 19:22

Most of these answers take an empirical approach; I would rather take an analytical approach, although I'm not the most qualified to do so.

There are two kinds of numbers that could be good choices:

• Large numbers that share common factors with many smaller numbers, and
• Small numbers that share common factors with few larger numbers

Large Numbers

The best large numbers will be the product of as many primes as possible, but to be in range all of them should be less than $\sqrt{max}$. The more primes you use the smaller they must be, and smaller primes have more multiples.

The obvious first-guess here is the largest Primorial within range, but for the range 2 through 1000 that guess is only 210(=2×3×5×7), which is too small. The largest multiple in range is 840.

A number chosen in this way will win against half of lower numbers by the common factor 2, against $\frac{1}{3}$ of the other half by the common factor 3, against $\frac{1}{5}$ of the remaining $\frac{1}{3}$ by the common factor 5, and so on. It will lose against the same proportion of higher numbers.

840 will win against about $\frac{27}{35}$ of the 838 lower numbers and about $\frac{8}{35}$ of the 159 higher numbers, for a win rate of about $\frac{683}{998}\approx 68.3\%$.

I'm not sure the 7 is helping, so I'm also going to look at the next lower primorial, 30(=2×3×5); the largest multiple in range is 990, which will win against $\frac{11}{15}$ of the 989 lower numbers and $\frac{4}{15}$ of the 10 higher numbers, for a win rate of about $\frac{730}{998}\approx 73.1\%$

So it looks like being higher is more valuable than being more composite - but 1000, a multiple of 2 and 5, only wins against $\frac{3}{5}$ of the 998 lower numbers (for a win rate of about 59.9%), so there must be a limit to that.

The largest multiple of 6(=2×3) in range is 996, with a win rate of about 66.4%, so it looks like the limit on the value of compositeness in this case is between 5# and 7#.

Small Numbers

In order to have few common factors with larger numbers, the best small numbers will be (relatively) small primes, near $\sqrt{max}$. Any prime will lose against all smaller numbers, but will only lose against larger numbers which are multiples of itself. If $\exists p_n = \sqrt{max}$, that prime will lose against $p_n-2$ lower numbers and $p_n-1$ higher numbers.

Moving down from $\sqrt{max}$ loses fewer below and more above, and moving up loses move below and fewer above. Exactly how much these change depends on the gaps between primes, which I don't know how to consider in general, but the difference in wins and losses can be expressed in terms of the gap.

Choose some $p_n = \sqrt{max}-gap$; $p_n$ will lose to $p_n-2$ lower numbers and $\lfloor\frac{max}{p_n}\rfloor-1$ higher numbers. Losses are minimized where $p_n+\lfloor\frac{max}{p_n}\rfloor$ is minimized. If $p_n$ could be chosen continuously, that would be at $p_n=\sqrt{max}$; since it cannot be chose continuously, it is the $p_n$ nearest to $\sqrt{max}$ and any neighboring $p_n$ which lose equally.

In this case $\sqrt{1000}\approx 31.6$, so the nearest primes are 31, 29, 37, and 23, in that order.

• 31($\approx\sqrt{max}$-0.6) will lose to 29 lower numbers and $\lfloor \frac{1000}{31} \rfloor-1 = 31$ higher numbers, for a win rate of $\frac{938}{998}\approx 94.0\%$
• 29($\approx\sqrt{max}$-2.6) will lose to 27 lower numbers and $\lfloor\frac{1000}{29}\rfloor-1 = 33$ higher numbers, for a win rate of $\frac{938}{998}\approx 94.0\%$
• 37($\approx\sqrt{max}$+5.4) will lose to 35 lower numbers and $\lfloor\frac{1000}{37}\rfloor-1 = 26$ higher numbers, for a win rate of $\frac{937}{998}\approx 93.9\%$
• 23($\approx\sqrt{max}$-8.6) will lose to 21 lower numbers and $\lfloor\frac{1000}{23}\rfloor-1 = 42$ higher numbers, for a win rate of $\frac{935}{998}\approx 93.7\%$

The two primes closest to $\sqrt{1000}$ are a tie for the winningest choice.

• Looks like I defined the win rate differently than in other answers; I used $\frac{wins}{wins+losses}$, it looks like other answers use $\frac{wins}{wins+losses+ties}$ – ShadSterling Oct 29 '14 at 3:07
• ... I also missed that 990=30*11*3, so it also wins against lower multiples of 11, which my calculation did not account for - it should be the best of the high numbers but a slightly higher margin. 996=6*83*2 also has that problem, but it's probably negligible. 840=210*4 isn't affected. – ShadSterling Oct 31 '14 at 21:08

29 and 31 are the best options.

Python code (brute force) :

import operator
def hcf(x,y):
if y>x:
return hcf(y,x)
if x%y == 0:
return y
else:
return hcf(y,x%y)

def check(m,c):
if m==c:
return 0
if hcf(m,c)>1:
return 1 if m>c else -1
else:
return -1 if m>c else 1

h = {}
a = range(2,1001)

for i in a:
h[i] = sum(check(i,j) for j in a)

sorted_h = sorted(h.items(), key=operator.itemgetter(1),reverse=True)
print "\n".join(map(str,sorted_h[0:10]))


Result :

(29, 878)
(31, 878)
(37, 876)
(41, 874)
(23, 872)
(43, 872)
(47, 868)
(19, 862)
(53, 862)
(17, 854)

• No, the second value 877 is the expected total score for the mathematician for all possible choices by the computer assuming win = 1, loss = -1 and draw = 0 [which actually took me some time to figure out :) ] – FacePalm Oct 30 '14 at 3:03

I wrote some Matlab code to solve the problem as well.

Optimized optimized code time spent 0.0085 s

clear all

tic;

N = 1000;                           % max range
pmax = primes(N);                   % primes in max range

x = zeros(N, 1);                    % pre-allocating x
y = zeros(N, 1);                    % pre-allocating y

for p = pmax
y(p, 1) = y(p, 1) + ceil(N - p - (N - p)/p);
end

for i = 1:N
x(i, 1) = sum(y(i, :));
end

% optional to change y into probability vector
xx = x./(N - 1);         % turns the entries in to probabilities

% Sort probabilities and find the top five
[sortedValues, sortIndex] = sort(xx(:), 'descend');
maxIndex = sortIndex(1:5)
maxValues = sortedValues(1:5)

Time = toc;


Optimized code time spent is 0.0385s

clear all

tic;

N = 1000;                           % max range
pmax = primes(N);                   % primes in max range

x = zeros(N, N);                    % pre-allocating x
y = zeros(N, 1);                    % pre-allocating y

for p = pmax
% p + 1 so we skip all ties
for i = (p + 1:1:N)
if mod(i, p) == 0
% mod is 0 loss
x(p, i) = x(p, i) - 1;
else
% mod is not 0 win
x(p, i) = x(p, i) + 1;
end
end
% map -1 to 0 for counting purposes
for i = 2:N
if x(p, i) < 0
x(p, i) = 0;
end
end
end

% count number of wins and store in vector
for m = pmax
y(m, 1) = y(m, 1) + sum(x(m, :));
end

% optional to change y into probability vector
yy = y./(N - 1);         % turns the entries in to probabilities

% Sort probabilities and find the top five
[sortedValues, sortIndex] = sort(yy(:), 'descend');
maxIndex = sortIndex(1:5)
maxValues = sortedValues(1:5)

Time = toc;


time 30 secish

clear all
close all

Max = 1000;                   % max number range

x = zeros(Max - 1, 1);        % pre-allocating x
k = zeros(Max, 1);            % pre-allocating k

% determine the number of times each number in range 2:Max wins
for i = 2:Max
for j = i + 1:Max         % skip values on diagonal and below i + 1:
n = gcd(i, j);
if n == 1
k(min(i, j), 1) = k(min(i, j), 1) + 1;
else
k(max(i ,j), 1) = k(max(i, j), 1) + 1;
end
end
x(i, 1) = sum(k(i));
end

% optional to change x into probability vector
x = x./(Max - 1);         % turns the entries in to probabilities

% Sort probabilities and find the top five
[sortedValues, sortIndex] = sort(x(:), 'descend');
maxIndex = sortIndex(1:5)
maxValues = sortedValues(1:5)


Print out of sortIndex and sortedValues:

maxIndex =

29
31
37
41
23

maxValues =

0.9389
0.9389
0.9379
0.9369
0.9359


Since we know the solution is a prime number, the code can be re-written to take advantage of this knowledge. Time 14 secish

Max = 1000;                            % max number range
% since we know the solution is a prime this should speed up compilation
pMax = primes(1000);                   % primes under 1000

x = zeros(length(pMax) - 1, 1);        % pre-allocating x
k = zeros(Max, 1);                     % pre-allocating k

% determine the number of times each number in range pMax wins
for i = pMax
for j = 2:Max
n = gcd(i, j);
if n == 1
k(min(i, j), 1) = k(min(i, j), 1) + 1;
else
k(max(i ,j), 1) = k(max(i, j), 1) + 1;
end
end
x(i, 1) = sum(k(i));
end

% optional to change x into probability vector
x = x./(Max - 1);         % turns the entries in to probabilities

% Sort probabilities and find the top five
[sortedValues, sortIndex] = sort(x(:), 'descend');
maxIndex = sortIndex(1:5)
maxValues = sortedValues(1:5)


Here, if we choose the number 2, we see all multiples of 2, if chosen by the computer will make us loose the game. Thus, $1000/2$, thus 500-1 = 449 possibilities to loose.

However, if choose the largest prime, that is 997, all smaller numbers would make us loose the game.

So this is basically an optimization problem, where we need to minimize multiples(of its factors), as well as have less numbers smaller than it.

So we check for prime numbers, which would have the least number of numbers larger than it, that can win from it, and we need a small enough prime, that the numbers smaller than it are not many.

We see $\sqrt{1000} = 31$ (approx)

So, we see our answer would be close to 31, so we check the primes near it.

29->61 31->61 37->62

So our answer would be 29 and/or 31.

The mathematician should choose 990. It has the most prime factors of 2, 3 ,5, 11 at the high end of numbers to 1000. Since the computer is going to pick a random number, there is a high probability (large set) 1. That the number picked will be less than 990 2. That it will get an easily factored number containing 2,3,5 or 11 (there are many).

The mathematician does not want to be in this large set, and should win most games.

• It has already been proven in several answers that the optimal choices are 31 and 29, though 990 is a good initial guess. – user139000 Oct 29 '14 at 5:27