Dirichlet’s Approximation Theorem in Mathematica

The Dirichilet’s approximation Theorem is about using rational numbers to approximate any numbers. To state it formally:

If \alpha \in \mathbb{R}, n \in \mathbb{Z}^+, then there exist integer a and b with 1 \leq a \leq n such that |a \alpha - b | < 1/n.

I managed to write up a simple Mathematica code (when I say simple, it still took me about 1 hr to write up and satisfy with, I suck at programming)




Declaring all the necessary variables first.

\[Alpha] = E;
n = 1000;
seq = (Range[n + 1] - 1)*\[Alpha];
fracseq = FractionalPart[seq];
intervals = Partition[(Range[n + 1] - 1)/n, 2, 1];
found = 0;
k = 1;

Just like in the proof of the Dirichlet Therorem, we created a sequence, which are the fractional part of first n integer multiples of alpha. The intervals are created in a similar way.

A loop to find the very first instance that we have “two piegeons in the same piegeonhole”.

While[k < Length[intervals] && found == 0, If[ Total[Boole[IntervalMemberQ[Interval[intervals[[k]]], fracseq]]] > 1,
found = 1, found = 0

This is just checking that we have the piegeons in the right place.

k - 1
Total[Boole[IntervalMemberQ[Interval[intervals[[k - 1]]], fracseq]]]
intervals[[k - 1]]



{3/250, 13/1000}

We have a selection criterion here, to select the first two piegeons in the hole.

selectseqcriterion =
Boole[IntervalMemberQ[Interval[intervals[[k - 1]]], fracseq]];
Pick[seq, selectseqcriterion, 1];

i = Pick[Range[n + 1] – 1, selectseqcriterion, 1][[1]];
j = Pick[Range[n + 1] – 1, selectseqcriterion, 1][[2]];

‘a’ and ‘b’ are the denominator and numerator of the rational approximation to the real number alpha.

a = j - i
b = Floor[j*\[Alpha]] - Floor[i*\[Alpha]]



This is just a check that the function works.

Abs[a*\[Alpha] - b] < 1/n
Abs[\[Alpha] - b/a] <= 1/(a^2)
N[Abs[\[Alpha] - b/a]]





So this simple program works. The exercise and theorem came from Kenneth Rosen’s “Elementary Number Theory”, and it certainly a good exercise to do.

Tossing A Die With A Choice

This week, a guy from Jane Street asked the following question:

  1. Suppose in Game 1, a dealer toss a fair die. The player receive money equal to whatever shows up on the die. (So, if a ‘4’ shows up, you receive $4). But in order to play this game, you must pay an “entry fee” first. What is the maximum amount of “entry fee” should you pay?
  2. Suppose in Game 2, we have the same rule (assume the entry fee is $3.5, which is the answer for a fair Game 1). But after every toss, you can observe the value on the die and have the choice to either:

a)     Receive the amount equal to whatever showed up on the die, which is the same as Game 1.

b)    Have a re-roll and receive whatever the new amount shows up on the die.

The question: is this option to have a re-roll more favourable game to play for the player?

The answer I immediately gave was: yes. This was just a pure reflex, because possessing more information “usually” improve the decision making process. Of course, this should be confirmed using probability, but I am reasonably comfortable with this intuition.

The argument using probability is very simple: If $3.5 is the entry fee, then anything less than $3.5 will be a “loss”. The probability of a “loss” is the sum of probabilities of getting a ‘1’, a ‘2’ and a ‘3’. P(win)=P(loss)=0.5, which is a fair game. Now, if we choose to have a re-roll whenever we encounter a ‘1’ or a ‘2’ or a ‘3’, then intuitively, we should maximise our probability of winning. Hence,

P(win)=P(win at first roll) + P(loss at the first roll)*P(win at second roll)

= 0.5  + 0.5 * 0.5

= 0.75

Which is what we should expect if we run a simulation in Mathematica.

When I started to program in Mathematica, I made several generalisations:

  1. Instead of using the strategy “re-roll if the die shows less than 3.5”, I produced a list, when the “cut-off value for a re-roll ranges from 1 to 6”.
  2. Instead of “entry fee is 3.5”, I used Manipulate function to display the probability of making a profit in the game, if the entry fee ranges from 1 to 6.

I am still unfamiliar with a lot of functions in Mathematica, for instance, I am sure there is a more efficient way of counting than using Length[Select[….]].

nexp = 10000;
 rolls = RandomInteger[{1, 6}, nexp];
 nnoroll = Length[Select[rolls, # > 3.5 &]]
 f[i_, entry_] := (switch = RandomInteger[{1, 6}, nexp];
     switch[[k]] < i,
     switch[[k]] = RandomInteger[{1, 6}];
    , {k, 1, Length[switch]}
   Length[Select[switch, # > entry &]]/nexp)

 Manipulate[{ListPlot[Table[f[i, entry], {i, 1, 6}], Filling -> Axis, 
    DataRange -> {1, 6}, PlotRange -> {0, 1}], 
   Table[f[i, entry], {i, 1, 6}] // N}, {entry, 1, 6, 1}]

The Birthday Problem (2)

So, this is an improved version of my previous Birthday Problem. 

The major improvement was that I used “DeleteDuplicates” instead of “Union”. 

The difference in computation time is quite remarkable. DeleteDuplicates took 29 seconds, while Union took about 72 seconds. 

f[k_] := (
  count = 0;
  repeat = 10000;
   a = RandomInteger[{1, 365}, k];
   b = DeleteDuplicates[a];
   If[Length[a] != Length[b], count++
   , {i, 1, repeat}
  prob = N[count/repeat];
list = Table[f[i], {i, 1, 365}];

So, yes! Sleep is help to gain new perspective about how to do things.

And yes, there is always room for improvement. Maybe the test using Length isn’t the most efficient. I shall explore more about this next time. 

Birthday Problem (1)

I am in the process of learning how to do some basic programming in Mathematica, and so, with my very limited abilities, I attempted to write a simulation in Mathematica today. And I thought the birthday problem seems like a good thing to model.

The problem itself is often phrased something along the lines of:

In a party of N people, how many people we should select, so that there is 50% chance that at least two people shares the same birthday?

We are assuming that there isn’t leap year (no one is born on 29 Feb) and each birthday is equally probable etc.

Few weeks back, I wrote a MATLAB program, and it was basic enough so that when I input a particular k value, it would give me a probability out of 100,000 experiments.


for i=1:times
    x = round(365*rand(k,1));
        if length(x) ~= length(unique(x))


Just as probability theory predicted, we expect selecting 23 people would be enough to give 50%.


So, that was MATLAB, but I desire a plot of points, which I could vary the value of k, something like this:


But I ran into some problems, because I am not particular familiar with how to manipulate vectors and functions in Mathematica.


 f[k_] := (
   count = 0;
   repeat = 10000;
    a = RandomInteger[{1, 365}, k];
    b = Union[a];
    If[Length[a] != Length[b], count++
    , {i, 1, repeat}
   prob = N[count/repeat];
 list = Table[f[i], {i, 1, 365}];

Birthday Problem

  1. The whole thing took about 75 seconds to complete. Which is much longer than I expected to be. Starting from 100, the values are so close to probability of 1, it is taking quite a long time to run the simulation. With some playing around, I found most of the time were spent on the “Union” function.
  2. My original intent was to create a function that list out each probability for each k value. I planned to do something like (as if I am using MATLAB)
    list=(list, prob);

    HoweverI never figure out how to do this concatenation. (Probably it is for the better, this repeated concatenation might further increase computation time)

  3. I managed to skip problem 2 with the self defined function f[k_]. However, it took me quite a while to figure out how to use Return[], which gave me exactly what I desire from the function f[k_].

So, the next thing I should do is to consider how to improve the efficiency of the program by modify the “Union” function into something faster.