90 views (last 30 days)

%shuffle a clean deck of 1e6 tails and heads

p=ones(1,500000); m=-ones(1,500000);C=[p ; m];a=C(:)';

%1st shuffle

A=a(1:500000);B=a(500001:1e6);C=[A ; B];C1=C(:)';

% next 1e3 shuffles

A=C1(1:5e5);B=C1(5e5+1:1e6);C=[A;B];C2=C(:)';

A=C2(1:5e5);B=C2(5e5+1:1e6);C=[A;B];C3=C(:)';

A=C3(1:5e5);B=C3(5e5+1:1e6);C=[A;B];C4=C(:)';

...

Brainlock :

How to use output of a command as input in following repeated commands

John D'Errico
on 12 Jul 2020 at 22:30

Edited: John D'Errico
on 12 Jul 2020 at 22:39

I put together this simple script. In it you will find perhaps the most efficient scheme for a riffle shuffle that I know of, since it can be done as a simple element permutation, just a re-indexing of the elements.

With N = 5e5, thus 1e6 total elements, you should see how fast it is.

timeit(@() B(riffleperm))

ans =

0.005221

So it would appear I can do nearly 200 such shuffles in one second.

Anyway, read down to the end, and then think about what you see.

%% Will a sequence of riffle shuffles ever produce a truly random sequence? Well, yes. And no.

% Consider a sequence of random bits, N zeros, followed by N ones. Now

% perform a perfect riffle shuffle. Repeat as often as you wish. The final

% sequence of bits after S such perfect riffle shuffles is perfectly

% deterministic, so not at all random. It may appear to be pseudo-random.

% It may approach the appearance of randomness at some level. We can test

% for things like that.

N = 5e5;

B0 = floor((0:2*N-1)/N);

% riffleperm has the property that B(riffleperm) will be a perfect riffle

% shuffle.

riffleperm = zeros(1,2*N);

riffleperm(1:2:end) = 1:N;

riffleperm(2:2:end) = N+1:2*N;

%% Now, perform S such riffleshuffles. For example, with S = 10, we might do this:

S = 10;

B = B0;

for ns = 1:S

B = B(riffleperm);

end

%% How random is B, after 10 such shuffles?

% One simple test might be to ask how many 1s are in the first half of the

% sequence. While a perfectly random sequence might seem to have exactly

% 1/2 of the 1's that would be wrong. In fact, the distribution of the

% number of 1's in either half of a truly random sequence would follow a

% binomial distribution.

sum(B(1:N))

ans =

250000

% In fact, that is a highly unlikely event. It is not at all random, that

% exactly 50% of the elements in the first half are 1, that is, IF the

% sequence was truly random. Worse, try the same thing for ANY number of

% perfect riffle shuffles.

S = 20;

B = B0;

for ns = 1:S

B = B(riffleperm);

disp([ns,sum(B(1:N))])

end

1 250000

2 250000

3 250000

4 250000

5 250000

6 249984

7 250000

8 250000

9 250000

10 250000

11 250000

12 249996

13 250000

14 250016

15 250000

16 250002

17 250000

18 249998

19 249992

20 249996

% Surprisingly, a fair number of those shuffles end up with exactly 50% of

% the elements 1 in each half o the sequence. Perhaps mor interesting yet,

% is to look at the number of elements in either half after EXACTLY 180

% such perfect riffle shuffles.

S = 180;

B = B0;

for ns = 1:S

B = B(riffleperm);

end

sum(B(1:N))

ans =

0

%% Yes, after 180 perfect riffle shuffles, that sequence of 1e6 bits

% really has returned to its original order. Had you started with

% B = 1:1e6;

% you would indeed find the final sequence in B is perfectly 1:1e6 again.

% 180 perfect riffle shuffles are sufficient there to completely and

% exactly restore the sequence of 1:1e6.

B = 1:2*N;

S = 180;

for ns = 1:S

B = B(riffleperm);

end

isequal(B,1:2*N)

ans =

logical

1

Now you should recognize something interesting. That is, after 1 perfect riffle shuffle, we will have identically the same sequence as 181 perfect riffle shuffles, since 180 shuffles returned us to the starting position. In fact, every 180 such perfect riffle shuffles are sufficient to return the sequence to the exact start point.

The most interesting question in all of this might be to show that for N = 5e5, thus a total sequence of 1e6 such elements, that 180 is the magic number of perfect riffle shuffles to restore the original sequence, when N = 5e5.

That probably would get more deeply into group theory then we want to go here.

John D'Errico
on 13 Jul 2020 at 13:59

There are also subtly alternative ways to shuffle a deck. But all of them can be turned into a shuffle as I did, by resequencing an index vector, and then using that index vector to perform the shuffle.

At a glance, the riffle shuffle you are using looks like an out shuffle. In those pages, you will see the number of perfect shuffles of that sort to return a deck with an even number of cards to its original order. In fact, they are even tablulated in the OEIS. Thus, for the outshuffle, apprarently it is sequence A002326 in the OEIS.

And again, once you have reached the number of shuffles to return the deck to order and the deck is in its original sequence, then futher shuffles act as if you started all over again at the very beginning.

An important thing to recognize is that for the problem at hand, since only 180 shuffles are used to restore the deck to perfect order, that there are only 180 distinct possible permutations of this deck, when using out shuffles. Of course, there are a vast number of possible permutations. For example, the number of ways we could shuffle those 1e6 bits is literally immense. With N = 5e5, we would have

factorial(2*N)/factorial(N)^2

We cannot compute that in MATLAB using double precision. Even using symbolic computations, it would be a number with a VAST number of digits. We can use the gamma function to compute the necessary power of 10 as:

>> (gammaln(2*N + 1) - 2*gammaln(N+1))/log(10)

ans =

300966.691648236

So a number on the order of 10^300966. there are roughly 301 THOUSAND decimal digits in the number of ways we could permute those bits. That is one BIG number.

However, we can never realize more than 180 possible such permutations using the riffle shuffles you have indicated. And that is effectively one SMALL number in context. So I don't think you would ever really call any permutation you will get here a truly random sequence. :) The only question is, how non-random is it, or how random does it look.

John D'Errico
on 17 Jul 2020 at 2:01

Your belief structure is not that relevant.

First, you need to accept that every perfect shuffle involves exactly the same resequencing of the cards. That is why the indexing solution I showed really is a valid shuffling sche,e. It is created by performing the shuffle on the numbers 1:2*N, Subsequent shuffles are the exact same permutation.

In fact, we can view a shuffle as a matrix multiply. This is why the corresponding matrix is called a permutation matrix. (A standard idea in mathematics.)

Consider a deck of size 8, so as you had it, N = 4.

N = 4;

% riffleperm has the property that B(riffleperm) will be a perfect riffle

% shuffle.

riffleperm = zeros(1,2*N);

riffleperm(1:2:end) = 1:N;

riffleperm(2:2:end) = N+1:2*N;

If you look at the sequence created, we see:

riffleperm

riffleperm =

1 5 2 6 3 7 4 8

In fact, that is exactly the out-shuffle you were using, but because I use a vector index, this is a faster way of doing it.

Now consider how it applies to the vector 1:8. According to the wikipedia page, on an 8 card deck, exactly 3 out-shuffles will restore the deck. We can find that same informatino from the OEIS - the online encyclopedia of integer sequences. Does it really work?

B = 1:8;

% after one shuffle:

B = B(riffleperm)

B =

1 5 2 6 3 7 4 8

% The second shuffle:

B = B(riffleperm)

B =

1 3 5 7 2 4 6 8

% and the third shuffle restores the sequence

B = B(riffleperm)

B =

1 2 3 4 5 6 7 8

So Wikipedia and the OEIS got it right. 3 perfect shuffles restore the sequence.

By the way, we can also view this as a matrix multiply. Again, they call it a permutation matrix for a good reason. In fact though, this would probably be a really inefficient way to do the permutation, as a matrix multiply, even if we used sparse matrices.

P = eye(8);

>> P = P(:,riffleperm)

P =

1 0 0 0 0 0 0 0

0 0 1 0 0 0 0 0

0 0 0 0 1 0 0 0

0 0 0 0 0 0 1 0

0 1 0 0 0 0 0 0

0 0 0 1 0 0 0 0

0 0 0 0 0 1 0 0

0 0 0 0 0 0 0 1

>> B*P

ans =

1 5 2 6 3 7 4 8

>> B*P^2

ans =

1 3 5 7 2 4 6 8

>> B*P^3

ans =

1 2 3 4 5 6 7 8

Yes, there it is again. 3 applications of P to the original sequence restored the sequence. Mathematics is a pretty thing some days. A nice way to visualise this is to see that for out-shuffle permutation, we have P^3 equals the 8x8 identiy matrix.

P^3

ans =

1 0 0 0 0 0 0 0

0 1 0 0 0 0 0 0

0 0 1 0 0 0 0 0

0 0 0 1 0 0 0 0

0 0 0 0 1 0 0 0

0 0 0 0 0 1 0 0

0 0 0 0 0 0 1 0

0 0 0 0 0 0 0 1

Note that not ALL perfect shuffles of an 8 card deck restore the vector in exactly 3 shuffles. For example, I'll generate a random permutation.

riffleperm = randperm(8)

riffleperm =

3 7 1 5 8 4 2 6

B

B =

1 2 3 4 5 6 7 8

P = eye(8);

P = P(:,riffleperm);

B*P

ans =

3 7 1 5 8 4 2 6

B*P^2

ans =

1 2 3 8 6 5 7 4

B*P^3

ans =

3 7 1 6 4 8 2 5

B*P^4

ans =

1 2 3 4 5 6 7 8

So with this other specific shuffle, 4 perfect shuffles of that form are sufficient to restore the 8 card deck. A 5th shuffle is equivalent to what you get after exactly one shuffle. We can show this by understanding that P^4 power is the identity matrix. Therefore, we should see that P^5=P^4*P, so a 5th such shuffle of this form will be equivalent to the first shuffle.

P^4

ans =

1 0 0 0 0 0 0 0

0 1 0 0 0 0 0 0

0 0 1 0 0 0 0 0

0 0 0 1 0 0 0 0

0 0 0 0 1 0 0 0

0 0 0 0 0 1 0 0

0 0 0 0 0 0 1 0

0 0 0 0 0 0 0 1

I hope you understand how these things work, that is permutations and how to view a perfect shuffle as a permutation, and finally how to incorporate that into a permutation matrix. The nice thing about the permutation matrix form is it is easy to understand how a permutation can be viewed as a matrix multily, and how powers of the permutation matrix then result in a repeated shuffle. It all holds together rather neatly. As I said before, we could even get more deeply into the group theory concepts of such things.

The same ideas apply to larger "decks", even a deck with 1e6 total elements in it. As it turns out, with 1e6 total elements, that specific shuffle does indeed repeat after 180 permutations. Of course, I don't want to raise a 1e6 square matrix to the 180'th power. but the vector index form is sufficient, and I imagine it would be faster. :)

N = 5e5;

riffleperm = zeros(1,2*N);

riffleperm(1:2:end) = 1:N;

riffleperm(2:2:end) = N+1:2*N;

B0 = 1:2*N;

B = B0;

for i = 1:180

B = B(riffleperm);

end

isequal(B,B0)

ans =

logical

1

Again MATLAB sees the result after exactly 180 shuffles are sufficient to recover the original sequence. The 200th shuffle would then be the same sequence as the 20th shuffle.

B20 = B0;

for i = 1:20

B20 = B20(riffleperm);

end

B200 = B0;

for i = 1:200

B200 = B200(riffleperm);

end

isequal(B20,B200)

ans =

logical

1

Again, this is how it works. A belief structure is not necessary, because you can always prove it to be fact. That is the beauty of mathematics.

Now, how can we show that 180 is the recovery "time" for a shuffle of this form on 1e6 total elements? Again, the wiki page suggests this is based on the multiplicative order of 2, modulo 2*N-1. 2*N-1 is used here because your shuffle involved 2*N cards. Effectively, this asks us to compute the smallest value of X, such that

mod(2^X,1e6-1) == 1

Be careful here, as large powers of 2 will rapidly exceed the dynamic range of double precision arithmetic. 2^180 is far beyond the limit, which is 2^53-1 to represent an integer exactly. However MATLAB offers the powermod function, which allows us to compute the necessary remainder as:

powermod(2,180,1e6-1)

ans =

1

Indeed, we see the remainder is 1. 180 is the smallest such power of 2 that does what we need.

Going back to the 8 card deck, all we need to know is:

mod(2^3,8 - 1)

ans =

1

So for that size deck, exactly 3 out-shuffles really should recover the original sequence.

It takes only a wee bit of work to show, for example, that given a deck of size

D = 12345678;

powermod(2,10136,D-1)

ans =

1

Therefore, given that size deck, exactly 10136 perfect out-shuffles will resoire the original deck to perfect order.

Again, belief is not necessary, just mathematics. And with the correct mathematics, you can prove what you need.

Steven Lord
on 12 Jul 2020 at 17:26

function out = shuffle(in)

out = in; % Make out the correct size and type first

out(1:2:end) = in(1:end/2);

out(2:2:end) = in((end/2)+1:end);

end

This is for an out shuffle on a deck with an even number of cards. Modifying it to handle an in shuffle is an exercise for the reader, as it modifying this to work for a deck with an odd number of cards (if required by the terms of your assignment, like if you wanted to shuffle the thirteen cards of an individual suit rather than all fifty-two cards of a standard poker deck.)

x = 1:10;

for numShuffles = 1:6

x = shuffle(x);

fprintf('The deck after shuffle %d is %s.\n', numShuffles, mat2str(x))

end

george simpson
on 12 Jul 2020 at 22:18

Steven Lord
on 16 Jul 2020 at 20:54

I suspect your question is along the line of “who is george simpson “, “what is he up to’, and how to get him to solve his own problem” ,’does he have a license to drive his Malab”,….

Do try accepting...

No. I want to know what problem you're trying to solve to understand if there's a more effective tool already in MATLAB that we can suggest to solve the problem you're trying to solve.

If you ask the question "How do I dig a hole?" the answer will differ depending on how big a hole you want. If you're planting flower bulbs, the answer is "Use a garden trowel." If you're putting in a tree, you probably want a shovel. If you're trying to lay foundation for a new house, the right tool for the job is probably an excavator.

But to answer the final question you asked, don't define new variables C2, C3, C4, etc. each time.

deck = [ones(1, 10), -ones(1, 10)];

numShuffles = 16;

history = repmat(deck, numShuffles+1, 1);

for whichDeckToProcess = 1:numShuffles

currentDeck = history(whichDeckToProcess, :);

A = currentDeck(1, 1:10);

B = currentDeck(1, 11:20);

C = reshape([A; B], 1, 20);

history(whichDeckToProcess+1, :) = C;

end

Opportunities for recent engineering grads.

Apply Today
## 3 Comments

## Direct link to this comment

https://www.mathworks.com/matlabcentral/answers/563711-when-does-randomness-enter-a-shuffled-deck-of-coin-tosses#comment_933710

⋮## Direct link to this comment

https://www.mathworks.com/matlabcentral/answers/563711-when-does-randomness-enter-a-shuffled-deck-of-coin-tosses#comment_933710

## Direct link to this comment

https://www.mathworks.com/matlabcentral/answers/563711-when-does-randomness-enter-a-shuffled-deck-of-coin-tosses#comment_933803

⋮## Direct link to this comment

https://www.mathworks.com/matlabcentral/answers/563711-when-does-randomness-enter-a-shuffled-deck-of-coin-tosses#comment_933803

## Direct link to this comment

https://www.mathworks.com/matlabcentral/answers/563711-when-does-randomness-enter-a-shuffled-deck-of-coin-tosses#comment_933980

⋮## Direct link to this comment

https://www.mathworks.com/matlabcentral/answers/563711-when-does-randomness-enter-a-shuffled-deck-of-coin-tosses#comment_933980

Sign in to comment.