Beyond Varsity Math: 
The red-and-blue-balls puzzle

An exploration of Diophantine equations

({{Robert K. Moniot}, {Fordham University}})

This work is licensed under the Creative Commons Attribution-ShareAlike 4.0 International License. To view a copy of this license, visit http://creativecommons.org/licenses/by-sa/4.0/.

Juggler image modified from https://commons.wikimedia.org/wiki/File:Juggling_man.jpg by Hhamilton / Creative Commons Attribution-Share Alike.

Introduction

“From a bag containing red and blue balls, two are removed at random. The chances are 50-50 that they will differ in color.
What were the possible numbers of balls initially in the bag?”

This is the National Museum of Mathematics’ Varsity Math puzzle 117.

We will solve that problem and then explore the solutions for other values of the probability.

The results of this exploration are in an article accepted for publication in the American Mathematical Monthly, “Solution of an Odds Inversion Problem,” by Robert K. Moniot, 2020 (to appear).

Solution of Varsity Math problem

Let r,b be the numbers of red and blue balls, respectively.  Probability of picking red, then blue:

P(red,blue)==(r/(b+r))·(b/(b+r-1))

Probability of picking blue, then red:

P(blue,red)==(b/(b+r))·(r/(b+r-1))

These are equal.  So, probability of picking balls of different colors are:

P(red,blue)+P(blue,red)==((2 b r)/((b+r)(b+r-1)))

Equate this to 50% or 1/2.

((2b  r)/((b+r)(b+r-1)))==(1/2)

Cross multiply and expand:

4b r==(b+r)(b+r-1)==(b^2)+2b r+(r^2)-(b+r)
(b^2)-2b r+(r^2)==b+r
(((b-r))^2)==b+r

This implies total number of balls b+r must be a square.  Call it (n^2).  Then

b+r==(n^2), b-r==±n

Taking the positive branch so b≥r, and adding and subtracting these, we obtain

2b==(n^2)+n=> b==(((n^2)+n)/2)       2r==(n^2)-n=> r==(((n^2)-n)/2)

So r and b are successive triangular numbers, whose sum is square.  That’s pretty cool!

First 10 solutions

TableForm[Table[{(((n^2)-n)/2),(((n^2)+n)/2),(n^2)},{n,1,10}],
TableHeadings->{None,{"r","b","r+b"}}]

The first row formally satisfies the equation, but is not admissible since there need to be at least 2 balls.

Generalizing the problem

That was fun.  Now, what can we say about other odds than 50-50?

We seek to answer the questions:

For a given odds, is there a solution to the problem?  That is, does there exist a pair of numbers of red and blue balls that produce a given probability of drawing balls of different colors?

If there is a solution, is the number of solutions finite or infinite?

For instances with a finite number of solutions, is there a way to list them all?

For instances with an infinite number of solutions, can we obtain a formula or recurrence that will generate as many solutions as desired?

Are the methods used able to find all solutions?  Or do they miss some?

All of these questions can be answered very satisfactorily.

Some examples

For the following probability ratios there is no solution.

({{(4/9), (5/6), (3/4), (3/8), (98058/176501)}})

For these probabilities there are only these solutions:

({{(8/15), ->, ({{2, 4}, {4, 6}, {7, 8}, {8, 8}}),         , (12/25), ->({{9, 16}})}})

For the following probabilities there are an infinite number of solutions:

({{Odds, Smallest Solutions}, {2/7, ({{1, 6}, {6, 30}, {30, 145}, {145, 696}, {696, 3336}})}, {6/13, ({{4, 9}, {9, 18}, {18, 34}, {70, 126}, {126, 225}})}, {4/51, ({{281063265339959257, 6593004860639355312}, {6593004860639355312, 154654550959684890576}, {154654550959684890576, 3627788942691955573225}})}})

Notice that a number from one solution often appears in the next.

Diophantine equation

Let the given probability be p/q, in lowest terms.  We are solving

((2b r)/((b+r)(b+r-1)))==(p/q)

Cross multiplying and rearranging, we obtain

p (b^2)-2(q-p)b r+p (r^2)-p b-p r==0

This is a Diophantine equation, i.e. an equation in integers.

We will call any pair of integers (r,b) satisfying Equation () a formal solution.

Not all formal solutions are admissible as solutions to the original problem.    For a solution to be admissible it must satisfy

r≥0, b≥0

(we cannot have a negative number of balls) and

r+b≥2

(we need at least two balls in the bag in order to be able to draw two out).

Trivial solutions

Equation () has three trivial solutions, i.e. formal solutions that hold regardless of the values of p and q:

(r,b)==(0,1), (0,0), (1,0)

These are not admissible solutions, but they will prove useful in finding admissible solutions.  
When plugged into the original odds equation,

((2b r)/((b+r)(b+r-1)))==(p/q)

they yield p/q==0/0.  This is undefined, which is how they are able to satisfy the equation for any values of p and q.

Symmetry

Equation () is symmetric in r and b:

If (r,b) is a solution then (b,r) is also.

To list just distinct solutions, we impose the condition r≤b.  (We already did that in solving the Varsity Math problem by only using the positive sign on n.)

Reverse search

To get a feel for the problem, I wrote a Python program to generate the odds for all values of red, blue less than 1000.

Pseudocode:

for b in {1..max}
  for r in {1..b}
    odds=(2 r b)/((r+b)(r+b-1))
    append {odds,{r,b}} to result
sort result by odds

Python has a Fraction class and unlimited-size integers, which make this easy.

The results showed interesting patterns, such as the reappearance of a number from one solution in the next noted earlier.

The recycling recurrence

Often solutions for a particular p/q occur in series in which the same number appears in two solutions, once as the smaller value and once as the larger, i.e.

(r,b), (b,x)

where r<b<x.  (The Varsity Math case shows this pattern, repeating indefinitely.)

By equating the probability ratios for these two solutions and imposing the condition r≠x one obtains this formula:

Solve[{((2 r b)/((r+b)(r+b-1)))==((2b x)/((b+x)(b+x-1))),x≠r},x]

I call this the recycling recurrence since the new solution re-uses one of the values in the old solution.

If r<b then the formula produces a larger solution, whereas if r≥b it goes the other direction.

The result x given by the formula is not necessarily integer.  One can show that if p==1 or p==2, the recurrence produces integers ad infinitum.  For p>2, it  goes fractional after a few steps, with only 2 or 3 integer solutions in any series.

Recycling recurrence examples

Here is an example of a series of solutions for the probability of drawing balls of different color equal to 2/7.

TableForm[RecurrenceTable[{r[i+1]==b[i],b[i+1]==b[i](b[i]-1)/r[i],r[1]==1,b[1]==6},{r[i],b[i]},{i,6}],
TableHeadings->{None,{"r","b"}}]

The recycling pattern continues indefinitely.

Here is an example of a series of solutions for which only three are integer pairs.  The probability ratio is 10/21:

TableForm[RecurrenceTable[{r[i+1]==b[i],b[i+1]==b[i](b[i]-1)/r[i],r[1]==2/5,b[1]==2},{r[i],b[i]},{i,5}],
TableHeadings->{None,{"r","b"}}]

Change of variables

Equation () simplifies considerably if we change variables.  Let

t==b+r, v==b-r

Plug in, and after some algebra, Equation () becomes

(q-2p)(t^2)+2p t-q (v^2)==0

Much simpler.  But we can do even better.

Complete the square on t.

To make the next steps clearer, set a==(q-2p):

a (t^2)+2p t-q (v^2)==0

First step is to make the coefficient of (t^2) square by multiplying by a:

(a^2)(t^2)+2p a t-q a (v^2)==0

Add (p^2) to both sides:

(a^2)(t^2)+2p a t+(p^2)-q a (v^2)==(p^2)
(((a t+p))^2)-q a (v^2)==(p^2)

Put back  a==(q-2p) and let

u==a t+p==(q-2p)t+p

Then the equation becomes

(u^2)-q(q-2p)(v^2)==(p^2)

We put Equations (), () and () into formulas for convenient use later, with p and q as parameters.

rbequation[r_,b_]:=p (b^2)-2(q-p)b r+p (r^2)-p b-p r==0

tvequation[t_,v_]:=(q-2p)(t^2)+2p t-q (v^2)==0

uvequation[u_,v_]:=(u^2)-q(q-2p)(v^2)==(p^2)

The discriminant D

Set D==q(q-2p), and f==(p^2).  Then Equation () is

(u^2)-D (v^2)==f

Call D the discriminant because it determines the character of the equation.

D<0 => ellipse

D>0 => hyperbola

If D==0 Equation () breaks down.  Recall that in deriving Equation ()  from Equation () (the t,v equation), we multiplied by a==(q-2p) to complete the square, so () is not valid when that is 0.  Equation ()  is a parabola for that case.  It was already solved as the Varsity Math problem, so no problem.

Character of equation vs p/q

D==q(q-2p)

Elliptic: D<0 => q<2p => (p/q)>(1/2)

Parabolic: D==0 => q==2p => (p/q)==(1/2)

Hyperbolic: D>0 => q>2p =>(p/q)<(1/2)

Elliptical case: D<0

For the elliptic case I have not found any methods other than direct search, i.e. testing all possibilities.  Since the ellipse is a finite curve, direct search can in principle find all solutions that exist.

It is helpful to graph an example.  Here is one, p/q==8/15, that has several integer solutions.  First we graph it in r-b space.  The grid lines show where solutions occur, for instance (4,6).

Plot[b/.Solve[rbequation[r,b]/.{p->8,q->15}],{r,-1,9},
PlotRange->{{-1,9},{-1,9}},AspectRatio->1,
GridLines->{{2,4,6,8},{2,4,6,8}},AxesLabel->{"r","b"}]

Elliptical case, example plotted in t,v-space

Now plot the p/q==8/15 example in t-v space.  This ellipse is aligned with the axes and is tangent to the horizontal axis.

Plot[t/.Solve[tvequation[t,v]/.{p->8,q->15}],{v,-5/2,5/2},
PlotRange->{{-10,10},{-2,18}},AspectRatio->1,
GridLines->{{-2,-1,1,2},{6,10,15,16}},AxesLabel->{"v","t"}]

Elliptical case, example plotted in u-v space

Now plot the p/q==8/15 example in u-v space.  This ellipse is aligned with the axes and centered on the origin.

Plot[u/.Solve[uvequation[u,v]/.{p->8,q->15}],{v,-5/2,5/2},
PlotRange->{{-10,10},{-10,10}},AspectRatio->1,GridLines->{{-2,-1,1,2},{-8,-7,-2,2,7,8}},AxesLabel->{"v","u"}]

Elliptical cases vs. p/q

Plot of ellipse in r-b space for a series of odds ratios of form p/(2p-1):

listforplot=Table[{b/.Solve[rbequation[r,b]/.q->2p-1,b]},{p,5}]; 
Plot[listforplot,
{r,-1/2,11/2},
PlotRange->{{-1/2,11/2},{-1/2,11/2}},AspectRatio->1,
PlotLegends->Table[p/(2p-1),{p,5}],AxesLabel->{"r","b"}]

For p/q==1, it is a circle.  As p/q->1/2, the ellipses elongate, approaching a parabola at p/q==1/2.

Elliptical cases in u-v space

Same set of cases plotted in u-v space.

listforplot=Table[u/.Solve[uvequation[u,v]/.q->2p-1,u],{p,5}];
Plot[listforplot,
{v,-2,2},
PlotRange->{{-11/2,11/2},{-11/2,11/2}},AspectRatio->1,
PlotLegends->Table[p/(2p-1),{p,5}],AxesLabel->{"v","u"}]

The widths grow slowly. This means a search on v is more efficient than a search on u or on r or b.

Symmetry allows the testing to be limited to positive u,v values (first quadrant arc).  Need to use both signs on u to get all solutions for t.  Can use only positive v to keep r≤b.

Elliptical case: range of r,b

p (b^2)-2(q-p)b r+p (r^2)-p b-p r==0               (, repeated)

Solve for the far endpoint where b==r:

p (r^2)-2(q-p)(r^2)+p (r^2)-p r-p r==0
4p (r^2)-2q (r^2)-2p r==0
r((2p-q)r-p)==0
(Subscript[r, max])==(p/(2p-q))==((p/q)/(2p/q-1))

Solve for p/q value for a given (Subscript[r, max]):

(p/q)==((Subscript[r, max])/(2(Subscript[r, max])-1))

For (Subscript[r, max])==5,  p/q==5/9.  For larger p/q ratios, solutions must have r,b<5.  We can exhaustively list them.

Note: if q==2p-1, then (Subscript[r, max])==p.  This implies there are solutions for all probabilities of form p/(2p-1).  The recycling recurrence from (p,p) gives neighboring solutions (p-1,p) and (p,p-1).  These three solutions at the far vertex of the ellipse are the symmetrical counterparts of the three trivial solutions around the origin.

Elliptical case: exhaustive enumeration

Table enumerating all distinct solutions for r,b≤5:

TableForm[Table[{{r,b,2r b/((r+b)(r+b-1))}},{b,1,5},{r,1,b}]]

Some of these ratios are below the threshold p/q==5/9 so the search is not exhaustive for them.

For ratios in the range 5/9 to 1 the list is exhaustive.  These ratios are:

(5/9), (4/7), (3/5), (2/3), (1/1)

All of their solutions appear in this table.  No other ratios in this range have any solutions.

For example: 3/4, 4/5, 5/6, 5/7, 6/7, 5/8, 7/8 all have no solution.

This covers nearly half the range of possible probabilities!

Elliptical case: direct search in u-v space

Maximum v is for u==0:

Solve[uvequation[u,v]/.u->0,v]

We choose the positive root,

(Subscript[v, max])==(p/(Sqrt[q(2p-q)]))==((p/q)/(Sqrt[2p/q-1]))

For an ellipse that extends to a million red balls,

(p/q)==((Subscript[r, max])/(2(Subscript[r, max])-1)),  (Subscript[r, max])==(10^6) =>   (p/q)==((10^6)/(2 (10^6)-1))

(This is very close to 1/2.)  Plugging this into the formula for (Subscript[v, max]) yields

Floor[(p/(Sqrt[q(2p-q)]))/.{p->(10^6),q->2 (10^6)-1}]

Searching this range for solutions takes 0.04 seconds on my laptop.  The only solutions found are for the far endpoint (b,b) and adjacent (b-1,b):

({{r, b}, {999999, 1000000}, {1000000, 1000000}})

Even quite large solutions are within reach of direct search.

Hyperbolic case: D>0

The hyperbolic case can be divided into two categories:

D square: solution by factoring the equation

D non-square: solution by continued fractions

Hyperbolic case, D square

If D is square, then Equation () can be factored:

(u^2)-D (v^2)==f -> (u-(Sqrt[D])v)(u+(Sqrt[D])v)==f

Here (Sqrt[D]) is integer by assumption.  The two terms on the left must correspond to divisors of  f.

Let divisors of  f be {(Subscript[d, 1]),(Subscript[d, 2]),...,(Subscript[d, k])}.  Set

(u-(Sqrt[D])v)==(Subscript[d, i]), (u+(Sqrt[D])v)==(f/(Subscript[d, i]))

Solution:

u==((f+(Subscript[d, i]^2))/(2(Subscript[d, i]))), v==((f-(Subscript[d, i]^2))/(2(Sqrt[D])(Subscript[d, i])))

Reject any solutions that are not integer.

Exchanging (Subscript[d, i]) and f/(Subscript[d, i]) simply changes the sign of v, so only the divisors (Subscript[d, i])≤(Sqrt[f]) need to be tested.

Negative divisors simply change sign of u, v, so only positive divisors need to be tested.

At most one solution (u,v) for each divisor, so total number of solutions is finite.  Some ratios have no solutions.

Hyperbolic case, D square example: p/q==12/25

Here is an example that has an admissible solution.

D==q(q-2p)==25(25-24)==25==(5^2)
f==(p^2)==(12^2)==144

Divisors of f: {1,2,3,4,6,8,9,12,16,18,24,36,48,72,144}

Too many for a demo.  Use only the one that gives a solution.

d==2 => u==((f+(d^2))/(2d))==((144+4)/(2 2))==(148/4)==37,  v==((f-(d^2))/(2(Sqrt[D])d))==((144-4)/(2 5 2))==(140/20)==7
t==((u-p)/(q-2p))==((37-12)/(25-24))==25
r==((t-v)/2)==((25-7)/2)==(18/2)==9, b==((t+v)/2)==((25+7)/2)==(32/2)==16

Hyperbolic case, D square example: p/q==4/9

D==q(q-2p)==9(9-8)==9==(3^2)
f==(p^2)==(4^2)==16

Divisors of f: {1,2,4,8,16}.

Trying them all, one finds none lead to an admissible solution.
This case has no solutions.

Hyperbolic case, D nonsquare

(u^2)-D (v^2)==f==(p^2)               (, repeated)

Divide through by (p^2):

((((u/p)))^2)-D((((v/p)))^2)==1, or (x^2)-D (y^2)==1

This is known as the Pell equation.  It has been well studied for a few centuries.

When D>0 is nonsquare, its solutions are found from x/y== a convergent of the continued fraction (Sqrt[D]).

These concepts are explained in the next few sections.

Continued fractions

Continued fraction of a real number r:

r==(Subscript[a, 0])+(1/((Subscript[a, 1])+(1/((Subscript[a, 2])+(1/(...))))))

where the (Subscript[a, i]) are integers.  Continued fractions are usually written in the compact form

r==[(Subscript[a, 0]),(Subscript[a, 1]),(Subscript[a, 2]),...]

There is a simple algorithm to compute continued fractions.

If r is rational, the fraction terminates.  Otherwise it continues indefinitely.

If r is the irrational root of a quadratic equation with rational coefficients, the continued fraction repeats.

Convergents

Terminating the continued fraction at any point yields a rational number.

(Subscript[a, 0])
(Subscript[a, 0])+(1/(Subscript[a, 1]))==(((Subscript[a, 0])(Subscript[a, 1])+1)/(Subscript[a, 1]))
(Subscript[a, 0])+(1/((Subscript[a, 1])+(1/(Subscript[a, 2]))))=(((Subscript[a, 2])((Subscript[a, 1])(Subscript[a, 0])+1)+(Subscript[a, 0]))/((Subscript[a, 2])(Subscript[a, 1])+1))
...

These are called convergents of the continued fraction for r.

There is a simple algorithm to generate the sequence of convergents.

The convergents are the best rational approximations to r using a denominator of that size.  Example: 22/7 for π.

Pell equation

(x^2)-D (y^2)==1

For D>0 nonsquare, the Pell equation always has a nontrivial solution.  (Trivial solution is x==±1, y==0.)

The smallest nontrivial solution is given by setting x/y equal to the convergent of (Sqrt[D]) at the end of the first repeat cycle if the repeat cycle is even length; otherwise the end of the second cycle.

Example: p/q==4/11,  D==11(11-2 4)==33

Convergents[(Sqrt[33])]

Mathematica indicates the start of the repeat cycle by repeating the argument.  Show that this solution works:

(23^2)-33 (4^2)

Solution:

x==(u/p)==23, y==(v/p)==4

Get u,v:

u==p x==4 23==92,  v==p y==4 4==16

Map u back to t==r+b.

t==((u-p)/(q-2p))==((92-4)/(11-2 4))==(88/3)

Fractional, not admissible.  But don’t lose hope.

Formula to generate more solutions

If ((Subscript[x, 1]),(Subscript[y, 1])) is the smallest solution of the Pell equation, other solutions can be generated by

(Subscript[x, n])+(Subscript[y, n])(Sqrt[D])==((((Subscript[x, 1])+(Subscript[y, 1])(Sqrt[D])))^n), n==1, 2, ...

When the expression on the RHS is expanded, after collecting terms it is always in the form of x+y(Sqrt[D]).

Example: p/q==4/11,  D==33, ((Subscript[x, 1]),(Subscript[y, 1]))==(23,4)

(((23+4(Sqrt[33])))^2)==(23^2)+2 4  23(Sqrt[33])+16 33==1057+184(Sqrt[33]) => (Subscript[x, 2])==1057, (Subscript[y, 2])==184

Verify that this is also a solution to the Pell equation.

(1057^2)-33 (184^2)

Convert back to r,b.

u==p x==4 1057==4228,  v==p y==4 184==736
t==((u-p)/(q-2p))==((4228-4)/(11-2 4))==(4224/3)==1408
r==((t-v)/2)==336, b==((t+v)/2)==1072

Infinite number of solutions

Example shows smallest solution of Pell equation does not necessarily give admissible (r,b) solution.

However, it can be proved that if ((Subscript[x, 1]),(Subscript[y, 1])) is the smallest nontrivial solution of the Pell equation, then the solution given by

((((Subscript[x, 1])+(Subscript[y, 1])(Sqrt[D])))^(2n)), n>0

always yields an admissible (r,b) solution.

Hence the number of solutions to our problem is always infinite for D>0 and non-square.

Recurrence starting from trivial solutions

The Pell recurrence can be unwound algebraically to work directly on r,b.

It can then be applied to the trivial solutions (0,0), (0,1), and (1,0) to generate solutions.

Starting from (0,0) gives the Pell solution series.  Starting from (0,1) or (1,0) gives neighbors related via the recycling recurrence.

We’re not done yet!

We have a method that always finds solutions to the hyperbolic, non-square D case.

However, it misses some solutions.  Example: p/q==20/41.  Solving the Pell equation yields

({{r, b}, {17280, 23680}})

Additional solutions found using the recurrence are only larger.  Yet the following are solutions:

({{r, b}, {16, 25}, {85, 120}, {120, 168}, {552, 760}})

Note that the first is not related to another via the recycling recurrence.  The next two are a doublet, not a triplet.
The Pell equation method can only yield solutions that are members of a triplet related by the recycling recurrence.
This is rooted in the fact that its solutions are the same as found by applying the Pell recurrence to the trivial solutions, which are three in number.

It can be shown that if p is a prime greater than 2, then the method based on the Pell equation is in fact complete.  But for composite p, there can be additional solutions.

Finding all solutions

Hua (1982) gives a method that finds all solutions to Equation ().

(u^2)-D (v^2)==f               (4, repeated)

If |f|<(Sqrt[D]), then if it has a solution, the solution will be found among the convergents of (Sqrt[D]), similarly to solving the Pell equation, but with the difference that the solution is not the last convergent in the cycle.  One has to search the convergents for solutions, and there may be more than one in the first cycle.

If this inequality is not satisfied, then we proceed as follows.

Hua, Luogeng (Loo Keng) (1982).  Introduction to Number Theory, translated from the Chinese by Peter Shiu. Springer-Verlag.  Chapter 11.

Method of solution

The method requires seeking integers l,h satisfying

(l^2)==D+f h

which can be solved by searching.  It suffices to search using the range -(Subscript[h, max])<h<(Subscript[h, max]) where

(Subscript[h, max])==Max(|(f/4)|,|(D/f)|)

Since f≥(Sqrt[D]) it is guaranteed that (Subscript[h, max])<\[LeftBracketingBar]f\[RightBracketingBar].  Thus the RHS is always reduced. Then solve

(u^2)-D (v^2)==h

If |h|<(Sqrt[D]) solve directly using continued fractions; otherwise repeat recursively.  Once one has a solution

(x^2)-D (y^2)==h

then solutions to (u^2)-D (v^2)==f are given by

u==((D y±l x)/h), v==((x±l y)/h)

One more important element is needed to ensure completeness.   The continued fraction method only yields solutions in which gcd(u,v)==1.  If gcd(u,v)==g>1, then (g^2) must divide f==(p^2), i.e. g must divide p.  To find those solutions, then, one can divide Equation () by each of the divisors of  p.  Thus, for each g that is a divisor of p, solve

((((u/g)))^2)-D((((v/g)))^2)==((((p/g)))^2)

Then multiply the solution by g to obtain (non coprime) u,v satisfying Equation ().

This method finds all solutions.

Appendix: functions to solve the equation

For the reader who would like to explore this problem with examples, I provide here functions that solve Equation () given any probability p/q in the range 0 to 1.  These functions draw on analysis that is not laid out in the foregoing summary of results.  A single function solvePuzzle[p/q,[n,[tableform]]] is provided that uses the appropriate function to solve for any ratio 0≤p/q≤1 and print up to n solutions (default 10), as a list or in table form (the default).

Function to solve the elliptical case by search

Argument z==p/q.  Include a test to proceed only for ratios in the elliptical range.  Optional argument n is the maximum number of solutions to print (default 10).  Only distinct solutions with r≤b are printed.

Outline of the function:

Compute range of v to search.

Compute positive u values from Equation () corresponding to each integer value of v in range.

Sift out integer values of u and pair with corresponding v values.  Note that v values are offset from Position values by 1 since v starts at 0 while Position starts at 1.  Augment list with -u values since they are also solutions and give distinct (r,b).  Use only non-negative v values to have distinct solutions.

Convert (u,v) solutions to (t,v).

Convert (t,v) solutions to (r,b).

Sift out positive integer values of r and b to yield only admissible solutions.  (The trivial solutions are not output.)

solveElliptical[z_,n_:10]:=Module[
{p,q,maxv,ufromv,uvaluesall,hits,uvalues,vvalues,tvalues,rvalues,bvalues,rbvalues},
p=Numerator[z];
q=Denominator[z];
If[1/2<z≤1, (* make sure this is an elliptical case *)
maxv=Floor[(( z)/(Sqrt[2 z-1]))];
ufromv=(Sqrt[(p^2)-q(2p-q)(v^2)]);
uvaluesall=Table[ufromv,{v,0,maxv}];
hits=Position[uvaluesall,_Integer,{1}];
uvalues=Join[Extract[uvaluesall,hits],-Extract[uvaluesall,hits]];
vvalues=Flatten[Join[hits-1,hits-1]];
tvalues=(u-p)/(q-2p)/.u->uvalues;
rvalues=Table[(t-v)/2/.{t->tvalues[[i]],v->vvalues[[i]]},{i,Length[vvalues]}];
bvalues=Table[(t+v)/2/.{t->tvalues[[i]],v->vvalues[[i]]},{i,Length[vvalues]}];
(* extract the admissible solutions *)
bvalues=Extract[bvalues,Position[rvalues,_Integer?Positive,{1}]];
rvalues=Extract[rvalues,Position[rvalues,_Integer?Positive,{1}]];
rbvalues=Sort[Table[{rvalues[[i]],bvalues[[i]]},{i,Length[rvalues]}]];
rbvalues[[;;Min[Length[rbvalues],n]] ](* print up to n solutions *)
, (* Else *) Print["argument out of range"]]
]

Examples

For 100% odds of drawing different colors, the only solution is 1 red and 1 blue.

solveElliptical[1]

Here is a ratio of form p/(2p-1) which hence has far-vertex solutions, around (p,p).  It also has other solutions.

TableForm[solveElliptical[200/399],TableHeadings->{None,{"r","b"}}]

This one has only the far-vertex solutions.

solveElliptical[3/5]

Here is one that has no admissible solution.

solveElliptical[4/5]

This ratio is not elliptical.

solveElliptical[1/9]

Function to solve the parabolic case

This has no argument z since it must be 1/2.  It has an optional argument n to specify how many solutions to print (default 10).  The solutions are pairs of successive triangular numbers.

solveParabolic[n_:10]:=Module[{},
Table[{((i(i-1))/2),((i(i+1))/2)},{i,2,n+1}]
]

Example

solveParabolic[]

Function to solve the hyperbolic square-D case

This function takes the probability ratio z==p/q as argument, with optional argument n to limit number of solutions printed (default 10).  Only distinct admissible solutions are printed.

Outline of function:

Make sure D>0 and square.  Handle special case p==0, D==1, for which all balls are one color.

Get list of positive divisors of f that are ≤f/2.

For each divisor, calculate u and v by the factoring method.  These are rational, not necessarily integer.

Sift for integer values of u and v.

Convert the list of (u,v) values to (t,v) values.  Sift again for integer values of t.

Convert list of (t,v) values to (r,b) values.  Sift again for integer values.

Clear[solveHyperbolicSquareD]

solveHyperbolicSquareD[z_,n_:10]:=Module[{p,q,D,f,dlist,
ucandidates,vcandidates,hits,uvalues,vvalues,tvaluesall,tvalues,rvalues,bvalues,rbvalues},
p=Numerator[z];q=Denominator[z];
D=q(q-2p);f=(p^2);
If[D>0&&Element[(Sqrt[D]),Integers],
If[D==1,
Table[{0,i},{i,2,n+1}],(* p==0 gives D==1, degenerate case *)

dlist=Divisors[f];
dlist=dlist[[1;;Ceiling[Length[dlist]/2]]];
ucandidates=Table[((f+(d^2))/(2d)),{d,dlist}];
vcandidates=Table[((f-(d^2))/(2 (Sqrt[D]) d)),{d,dlist}];
(* first pass: pick out positions where u is integer *)
hits=Position[ucandidates,_Integer];
uvalues=Extract[ucandidates,hits];
vvalues=Extract[vcandidates,hits];
(* second pass: pick out positions where v is integer.  It is always ≥0 *)
hits=Position[vvalues,_Integer];
uvalues=Extract[uvalues,hits];
vvalues=Extract[vvalues,hits];
tvaluesall=(u-p)/(q-2p)/.u->uvalues;
hits=Position[tvaluesall,t_/;t>1&&IntegerQ[t]]; (* require admissible *)
tvalues=Extract[tvaluesall,hits];
vvalues=Extract[vvalues,hits];
rvalues=(tvalues-vvalues)/2;
bvalues=(tvalues+vvalues)/2;
(* r,b may still be fractional if t, v not same parity *)
bvalues=Extract[bvalues,Position[rvalues,_Integer,{1}]];
rvalues=Extract[rvalues,Position[rvalues,_Integer,{1}]];
rbvalues=Sort[Table[{rvalues[[i]],bvalues[[i]]},{i,Length[rvalues]}]];
rbvalues[[;;Min[Length[rbvalues],n]] ](* print up to n solutions *)
]
,
(* else *) Print["D=",D, " not OK"]]
]

Exercise the discriminant test by giving it a hyperbolic ratio with nonsquare D.

solveHyperbolicSquareD[1/3]

Examples

This one has no solutions.  D==9.

solveHyperbolicSquareD[4/9]

This one has only one solution.  D==25.

solveHyperbolicSquareD[12/25]

Degenerate case, p/q==0/1, D==1, hyperbola coincides with axes.  Give first 5 admissible solutions.

solveHyperbolicSquareD[0,5]

Function to solve the Pell equation

Here is a function to find the primary, i.e., the smallest nontrivial, solution of the Pell equation (h^2)-D (k^2)==1 where D>0 is nonsquare.  Argument is D==q(q-2p).

Outline of function:

Check D>0 and nonsquare

Calculate first cycle of continued fraction using the Mathematica function ContinuedFraction.  The result is of form {(Subscript[a, 0]),(Subscript[a, 1]),...,(Subscript[a, k]),{(Subscript[b, 1]),(Subscript[b, 2]),...,(Subscript[b, r])}} where the (Subscript[b, i]) are the repeating cycle.

If the repeat length r is even, then this result, flattened to remove the structure, is passed to the Mathematica function Convergents to produce the list of convergents.  The r-th convergent gives the Pell solution.

If the repeat length r is odd, then the repeating cycle is appended to the result, which is flattened and passed to Convergents.  The 2r-th convergent gives the Pell solution.

Return numerator and denominator of the chosen convergent.  Since this is the base case, they are named (h,k).

solvePell[D_]:=Module[{cf,repeat,replen,hoverk} ,
If[D>0 &&(Sqrt[D])∉Rationals,
cf=ContinuedFraction[(Sqrt[D])];
repeat=cf[[-1]]; (* repeat cycle is a sub-list as last element *)
replen=Length[repeat];
If[Mod[replen,2]==0,
hoverk=Convergents[Flatten[cf]][[;;replen]], (* even replen *)
hoverk=Convergents[Join[Flatten[cf],repeat]][[;;2replen]] (* odd *)
];
{Numerator[hoverk[[-1]]],Denominator[hoverk[[-1]]]}  (* return last h,k as solution *)
,(* else *) Print["D=",D, " not OK"]
]
]

Example

solvePell[33]

Function to solve the hyperbolic non-square D case

This is the most complex case, and takes a team of functions to solve it.

First we define a pair of functions to find the solutions (u,v) using Hua’s method of recursive reduction.  Then another function uses those solutions and the Pell recurrence to get as many r,b values satisfying the original problem as desired.  A couple of other utility functions are needed.

Function to solve for (u,v)

Since the method of continued fractions only yields relatively prime solutions, it is necessary to also solve  Equation () divided by the square divisors of f, in order to find all solutions.

The recursive reduction routine includes a section that checks whether the RHS of the equation being solved is a perfect square, and if so includes the trivial solution ((Sqrt[f]),0).  This is not necessary in order to find all solution classes, but if it is not included, then fundamental solutions of some classes are not found.  (Solution classes are discussed in the next section.)

Either sign of u,v satisfies Equation ().  The method seeks only positive values.  For the hyperbolic case, negative u values yield negative r,b, which are not admissible.  Only positive values of v are returned.  The admissibility tests use Positive as the criterion to exclude the trivial solutions; replacing that by NonNegative would include them.

The following function implements the recursion to find (u,v) solutions.  It is meant to be called by the next function, a driver that collects all the solutions for Equation () and the sub-equations from dividing by square divisors of f.

Outline of the function:

If f<(Sqrt[D]), search the convergents of (Sqrt[D]) to find solutions of (u^2)-D (v^2)==f.  If f is a square, include the trivial solution u==(Sqrt[f]), v==0.

If f≥(Sqrt[D]), reduce the RHS using Hua’s method, and call self recursively to solve that equation.

solveRecursiveReduction[D_,f_]:=Module[
{l,h,lvalues,hvalues,contfrac,repeat,replen,convergents,rruvvalues,uvvalues,uvlist,r,s,unextgen,wnextgen},
If[(f^2)<D, (* solution by convergents is guaranteed *)
contfrac=ContinuedFraction[(Sqrt[D])]; (* List is of form {a0, a1, ... an,{b0, b1, ... bm}} *)
repeat=contfrac[[-1]]; (* extract the repeat cycle at end *)
replen=Length[repeat];
If[Mod[replen,2]==0,
convergents=Convergents[Flatten[contfrac][[;;replen]]], (* even replen *)
convergents=Convergents[Join[Flatten[contfrac],repeat][[;;2replen]]] (* odd *)
];
(* turn convergent fractions into {u,v} values *)
uvvalues=Transpose[{Numerator/@convergents,Denominator/@convergents}];
(* find {u,v} values satisfying the equation *)
uvvalues=Select[uvvalues,(u^2)-D (v^2)==f/.{u->#[[1]],v->#[[2]]}&];
(* If f is square, include the trivial solution u==(Sqrt[f]),v=0 *)
If[(Sqrt[f])∈Integers,
uvvalues=Append[uvvalues,{(Sqrt[f]),0}];
];
uvvalues (* return result *)

,(* else |f|≥(Sqrt[D]) *)
(* recursive reduction of RHS.  Seek l, h satisfying (l^2)==D+f h *)
lvalues=Flatten[Position[Table[Mod[(l^2)-D,f]==0,{l,Ceiling[f/2-1]}],True]];
hvalues=Table[((l^2)-D)/f,{l,lvalues}];
uvlist={};
Do[(* solve with reduced RHS for each l,h pair *)
l=lvalues[[j]]; h=hvalues[[j]];
rruvvalues=solveRecursiveReduction[D,h];
uvvalues=Flatten[Table[{(D v+δ l u)/h,(u+δ l v)/h}/.
{u->rruvvalues[[i]][[1]],v->rruvvalues[[i]][[2]]},
{δ,{-1,1}},{i,Length[rruvvalues]}],1];
(* Sign of u, v can be either, so take absolute value.  Remove duplicates to fend off combinatorial explosion. *)
uvlist=Sort[DeleteDuplicates[Abs[Join[uvvalues,uvlist]]]]
,{j,Length[hvalues]} (* Do control *)
];
uvlist (* return result *)
]
]

Next is the top-level driver of the recursive routine.  It calls it to find solutions of Equation () divided by the square divisors of f.

Outline of function:

Make sure D>0 and nonsquare.

For each square divisor d of f, call the previous routine to solve ((((u/d)))^2)-((D((v/d)))^2)==(f/(d^2)) and obtain solution to Equation () by multiplying the solution by d.

Delete duplicate solutions.

solveuvByRecursiveReduction[z_]:=Module[
{p,q,D,f,uvvalues,divisorlist},
p=Numerator[z];q=Denominator[z];
D=q(q-2p); f=(p^2);
If[ D>0 &&(Sqrt[D])∉Rationals ,
divisorlist=Divisors[p];
uvvalues={};
Do[
uvvalues=Join[uvvalues,
Cases[d(solveRecursiveReduction[D,f/(d^2)]),{_Integer,_Integer}]];
,{d,divisorlist}];
DeleteDuplicates[Sort[uvvalues]] (* return result *)
,(* else *) Print["D=",D, " not OK"]
]
]

Function to reduce set of (u,v) solutions to representatives of distinct classes

Solutions of Equation () can be put into classes based on whether they are related through the Pell recurrence.  See Nagell (1964), Section 58.

A fundamental solution of a class is the member of a class for which u≥0 is least, and conjugate solutions are obtained by changing the sign of v.

Outline of the function:

Calculate f==(u^2)-D (v^2) using the first solution in the given list.

Ensure u≥0 for all solutions and augment list with conjugates by negating v.

Sort into ascending order so that the fundamental solutions are first.  Delete duplicates such as any with v==0 or if the list already contained conjugates.

Initialize list uvclass of distinct class members to empty.

For each {u,v} in the list,

If it is not in the same class as any element of uvclass, append it to uvclass.

Return the list of distinct class members, sorted in ascending order.

getClasses[D_,uvlist_]:=Module[{f,uvall,uvclass},
f=((uvlist[[1]][[1]])^2)-D ((uvlist[[1]][[2]])^2);
uvall=Sort[DeleteDuplicates[Join[
Table[{Abs[uv[[1]]],uv[[2]]},{uv,uvlist}],
Table[{Abs[uv[[1]]],-uv[[2]]},{uv,uvlist}]]]];
uvclass={};
Do[
If[Position[Table[(* test each given in uvall against growing set of class representatives *)
((u U-D v V)/f∈Integers &&(v U-u V)/f∈Integers )/.{u->uv[[1]],v->uv[[2]],U->UV[[1]],V->UV[[2]]},{UV,uvclass}]
,True]=={},
uvclass=Append[uvclass,uv]
]
,{uv,uvall}];(* end Do *)
(* return *)  uvclass
]

Example: p/q==10/21.

getClasses[21,{{10,0},{11,1}}]

Trygve Nagell, Introduction to Number Theory, 2nd Ed., Chelsea Publishing Co., New York (1964).

Function to find (r,b) solutions

Now the function to generate (r,b) solutions.  Optional argument n is the maximum number of solutions to print (default 10).  It first gets the set of all fundamental class solutions.  It applies the Pell recurrence to these to generate more solutions, until the required number n of (r,b) solutions has been obtained.

The Pell recurrence for solutions of Equation () is obtained by equating rational and irrational parts of each side of the following equation:

(Subscript[u, n+1])+(Subscript[v, n+1])(Sqrt[D])==(h+k(Sqrt[D]))((Subscript[u, n])+(Subscript[v, n])(Sqrt[D]))

where (h,k) is the primary solution of the Pell equation.  Thus

(Subscript[u, n+1])==h (Subscript[u, n])+D k (Subscript[v, n])
(Subscript[v, n+1])==k (Subscript[u, n])+h (Subscript[v, n])

Outline of function:

Solve Pell equation (h^2)-D (k^2)==1.  (The result will be empty if D is not positive and nonsquare, in which case skip the rest.)

Solve Equation () using the above functions.  Result is a list of (u,v) pairs.

Reduce the list of (u,v) pairs to one representative of each class.  Conjugate solutions with v<0 are included in order to generate all solutions in succeeding generations.

Initialize list of (r,b) solutions to the empty set.

Loop until sufficient solutions have been found:

Calculate list of (r,b) values from list of (u,v) solutions.  Sift for integer positive values and append to growing list of solutions.

Replace list of (u,v) pairs by next generation obtained using the Pell recurrence.

solveHyperbolicNonsquareD[z_,n_:10]:=
Module[{p,q,D,hkpell,h,k,uvsolns,rbsolns,newrb},
p=Numerator[z];
q=Denominator[z];
D=q(q-2p);
hkpell=solvePell[D]; (* get the primary solution *)
If[Length[hkpell]==2,
{h,k}=hkpell;
uvsolns=getClasses[D,Evaluate[solveuvByRecursiveReduction[z]]];
rbsolns={};
While[Length[rbsolns]<n,newrb=Table[(1/2){(((p-u)/(2 p-q))-v),(((p-u)/(2 p-q))+v)}/.{u->uv[[1]],v->uv[[2]]},
{uv,uvsolns}]; (* calculate (r,b) values from (u,v) values *)
rbsolns=Join[rbsolns,DeleteDuplicates[Sort[Cases[Sort/@newrb,{_Integer?Positive,_Integer?Positive}]]]];
uvsolns=Table[{h u+D k v,k u+h v}/.{u->uv[[1]],v->uv[[2]]},
{uv,uvsolns}] (* replace list of (u,v) with next generation *)
];
rbsolns[[1;;n]]
]
]

Examples

TableForm[solveHyperbolicNonsquareD[4/11],TableHeadings->{None, {"r","b"}}]

Here is an example that has singlet and recycling doublet solutions missed by the Pell equation method.

TableForm[solveHyperbolicNonsquareD[20/41,12],TableHeadings->{None,{"r","b"}}]

Function to solve the problem for all cases

This function uses the appropriate method based on D to solve Equation () for any given probability from 0 to 1.  First argument is the probability p/q, optional second argument n is the limit of how many solutions to print (default 10).  Optional third argument tableform to format results as a table (default True).  Only distinct (r≤b) admissible solutions are printed.

solvePuzzle[z_,n_:10,tableform_:True]:=Module[{p,q,D,rbsolns},
If[z≥0&&z≤1,
p=Numerator[z];
q=Denominator[z];
D=q(q-2p);
If[D<0, rbsolns=solveElliptical[z,n]
,If[D==0, rbsolns=solveParabolic[n]
,If[D>0&&Element[(Sqrt[D]),Integers],rbsolns=solveHyperbolicSquareD[z,n]
,rbsolns=solveHyperbolicNonsquareD[z,n]
]
]
];
If[tableform,
TableForm[rbsolns,TableHeadings->{None,{"r","b"}}],
rbsolns
]
,(* else z<0 or z>1 *)
Print[z," is out of range"]
]
]

Examples

The Varsity Math case.

solvePuzzle[1/2]

An elliptical case.

solvePuzzle[4/7]

An elliptical case with no solution.

solvePuzzle[4/5]

Special case of unit probability.

solvePuzzle[1]

Special case of 0 probability.

solvePuzzle[0,5]

Hyperbolic square D case with no solutions.

solvePuzzle[4/9]

Hyperbolic square D case with one solution.  Format as list instead of table.

solvePuzzle[21/50,5,False]

Hyperbolic nonsquare D case with infinite number of solutions.  Produce the first 10.

solvePuzzle[5/13]