Lazy Larry's Random Algorithm

There are $4!=24$ permutations of the list, and each is possible with Larry's algorithm. This is a lot to analyze, but fortunately, we can organize these permutations into a set of "states":

• State P: Exactly two of the elements in the list are out of order. Example: bandicoot, aardvark, capybara, dingo

• State Q: All four of the elements in the list are out of order, but they can be placed in order by switching two pairs. This is the state the list is in before Larry runs the algorithm: capybara, dingo, aardvark, bandicoot

• State R: Exactly three of the elements in the list are out of order. Example: capybara, aardvark, bandicoot, dingo

• State S: Exactly four of the elements are out of order, and the order cannot be restored by switching two pairs. Example: bandicoot, capybara, dingo, aardvark

• State T: All four of the elements are in order: aardvark, bandicoot, capybara, dingo

Now we can analyze how we can move between states with each iteration of the algorithm. There are 6 possible switches with each iteration. Suppose we are in state P, like the example given above: bandicoot, aardvark, capybara, dingo

There is 1 switch that will bring us to state Q:

bandicoot, aardvark, capybara, dingo $\Rightarrow$ bandicoot, aardvark, dingo, capybara

There are 4 switches that will bring us to state R:

bandicoot, aardvark, capybara, dingo $\Rightarrow$ capybara, aardvark, bandicoot, dingo

bandicoot, aardvark, capybara, dingo $\Rightarrow$ dingo, aardvark, capybara, bandicoot

bandicoot, aardvark, capybara, dingo $\Rightarrow$ bandicoot, capybara, aardvark, dingo

bandicoot, aardvark, capybara, dingo $\Rightarrow$ bandicoot, dingo, capybara, aardvark

There is 1 switch that will bring us to state T:

bandicoot, aardvark, capybara, dingo $\Rightarrow$ aardvark, bandicoot, capybara, dingo

This gives us probabilities of $\frac{1}{6},$ $\frac{2}{3},$ and $\frac{1}{6}$ to go from state P to states Q, R, and T, respectively. The other states can be likewise analyzed. Below is a graph showing their relationships:

Let $X$ be the number of iterations until we reach state T. Then we have the expected values:

$\begin{aligned} p = \text{E}[X \mid P] \\ q = \text{E}[X \mid Q] \\ r = \text{E}[X \mid R] \\ s = \text{E}[X \mid S] \\ \end{aligned}$

We can write a system of equations using linearity of expectation .

$\begin{aligned} p &= 1+\frac{1}{6}q+\frac{2}{3}r \\ q &= 1+\frac{1}{3}p+\frac{2}{3}s \\ r &= 1+\frac{1}{2}p+\frac{1}{2}s \\ s &= 1+\frac{1}{3}q+\frac{2}{3}r \\ \end{aligned}$

Solving this system gives $q=27,$ and so the expected value of the number of iterations from the initial list is $\boxed{27}.$

Below you will find my code doing a brute-force Monte-Carlo simulation. The code uses multi-threading.

  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
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147

//Monte-Carlo simulation for https://brilliant.org/problems/lazy-larrys-random-algorithm/

#include <iostream>
#include <string>
#include <random>
#include <vector>
#include <chrono>
#include <thread>
#include <ctime>
#include <algorithm>

template <typename DataType, typename RandomEngineType>
std::uintmax_t randomSort(const DataType& inData, RandomEngineType& randomEngine)
{
    DataType curData = inData;
    std::uintmax_t curStep = 0;
    const std::size_t dataLen = inData.size();
    std::uniform_int_distribution<> posGenerator(0, dataLen-1);
    while(!std::is_sorted(curData.begin(), curData.end()))
    {
        const std::size_t posA = posGenerator(randomEngine);
        std::size_t posB = posGenerator(randomEngine);
        while(posA == posB)
        {
            posB = posGenerator(randomEngine);
        }
        const auto tmpVal = curData[posA];
        curData[posA] = curData[posB];
        curData[posB] = tmpVal;
        ++curStep;
    }
    return curStep;
}

template <typename DataType, typename RandomEngineType>
std::uintmax_t randomSortMcSimulation(const DataType& testData,
                                      const std::uintmax_t localMaxIter,
                                      const int inSeed)
{
    RandomEngineType randomEngine(inSeed);
    std::uintmax_t totalSum = 0;
    for(std::uintmax_t curIter = 0; curIter < localMaxIter; ++curIter)
    {
        totalSum += randomSort(testData, randomEngine);
    }
    return totalSum;
}

template<typename DataType, typename RandomEngineType>
class Worker
{
    private:
        DataType testData;
        std::uintmax_t localIter;
        int seed;
        std::uintmax_t totalSum;
    public:
        explicit Worker(const DataType& _testData,
                        const std::uintmax_t _localIter,
                        const int _seed) :
        testData(_testData),
        localIter(_localIter),
        seed(_seed),
        totalSum()
        {}
        void operator()()
        {
            this->totalSum = randomSortMcSimulation<DataType, RandomEngineType>(this->testData,
                                                                                this->localIter,
                                                                                this->seed);
        }
        std::uintmax_t GetTotalSum() const
        {
            return this->totalSum;
        }
};

template<typename UnsignedIntType>
UnsignedIntType gcd(const UnsignedIntType _a, UnsignedIntType _b)
{
    UnsignedIntType a, b;
    if(_a > _b)
    {
        a = _a;
        b = _b;
    }
    else
    {
        a = _b;
        b = _a;
    }
    while(b > 0)
    {
        const UnsignedIntType div = a%b;
        a = b;
        b = div;
    }
    return a;
}

template <typename DataType, typename RandomEngineType>
long double randomSortMcSimulation(const DataType& testData,
                                   const std::uintmax_t totalIter,
                                   const unsigned numberOfThreads)
{
    if(totalIter%std::uintmax_t(numberOfThreads) != 0)
    {
        throw std::invalid_argument("totalIter must be a multiple of numberOfThreads.");
    }
    const std::uintmax_t localIter = totalIter/std::uintmax_t(numberOfThreads);
    std::vector<std::thread> threadVector;
    threadVector.reserve(numberOfThreads);
    std::vector<Worker<DataType, RandomEngineType> > workerVector;
    workerVector.reserve(numberOfThreads);
    for(unsigned i = 0; i < numberOfThreads; ++i)
    {
        workerVector.push_back(Worker<DataType, RandomEngineType>(testData, localIter, i+time(NULL)));
        threadVector.push_back(std::thread(std::ref(workerVector[i])));
    }
    for(auto & th : threadVector)
    {
        th.join();
    }

    auto sumFun = [](const std::uintmax_t tmpSum, const Worker<DataType, RandomEngineType>& inWorker)->std::uintmax_t
    {
        return tmpSum+inWorker.GetTotalSum();
    };
    const std::uintmax_t totalSum = std::accumulate(workerVector.begin(), workerVector.end(), 0ULL, sumFun);
    const std::uintmax_t div = gcd(totalSum, totalIter);
    return static_cast<long double>(totalSum/div)/static_cast<long double>(totalIter/div);
}

int main()
{
    const unsigned numberOfThreads = std::thread::hardware_concurrency();
    const std::uintmax_t totalIter = std::uintmax_t(numberOfThreads)*5000000ULL;
    std::vector<std::string> testData{"capybara", "dingo", "aardvark", "bandicoot"};

    auto startTime = std::chrono::steady_clock::now();
    const auto res = randomSortMcSimulation<std::vector<std::string>, std::mt19937_64>(testData, totalIter, numberOfThreads);
    auto endTime = std::chrono::steady_clock::now();
    std::chrono::duration<double> runTime = endTime-startTime;
    std::cout<<"Runtime: "<<runTime.count()<<" [s]"<<std::endl;
    std::cout<<"Expected number of steps: "<<res<<std::endl;
    return 0;
}