Eulers Egg, Random Numbers and Integration

The yellow line ( $y=x+2$ ) intersects the circle $(x+2)^2+y^2=16$ at ( $2√2-2, 2√2$ ). Therefore the radius of the smallest red circle, $r$ , is given by $r^2=2(2√2-2)^2=8(3-2√2)$ . Area of the region within the egg is twice the sum of three areas : $A_1=\dfrac{πr^2}{8}=π(3-2√2)$ , $A_2=\dfrac{π\times {4^2}}{8}-2=2π-2$ , and $A_3=\dfrac{π\times {2^2}}{4}=π$ . The area is thus $π(12-4√2)-4=\boxed {15.927580090444...}$

Below is my parallel C++ code. I have found the radius of the smallest red circle in a way similar to the one described by @Alak Bhattacharya . I use the symmetry: the experiment takes place only for the positive $x$ coordinated.

  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
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183

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

template<std::size_t Dimension, typename ValueType, typename RandomEngineType, typename IsInsideType>
class McWorker
{
    private:
        const IsInsideType isInside;
        const std::uintmax_t localIter;
        std::array<std::uniform_real_distribution<ValueType>, Dimension> pointGenerator;
        RandomEngineType randomEngine;
        ValueType amountInside;

    public:
        explicit McWorker(const std::array<std::pair<ValueType, ValueType>, Dimension>& _limitBox,
                          const IsInsideType& _isInside,
                          const std::uintmax_t _localIter,
                          const int _seed) :
        isInside(_isInside),
        localIter(_localIter),
        pointGenerator(),
        randomEngine(_seed),
        amountInside()
        {
            for(std::size_t curDim = 0; curDim < Dimension; ++curDim)
            {
                const auto curLimit = _limitBox[curDim];
                this->pointGenerator[curDim] = std::uniform_real_distribution<ValueType>(curLimit.first,
                                                                                         curLimit.second);
            }
        }
        void operator()()
        {
            for(std::uintmax_t curIter = 0; curIter < this->localIter; ++curIter)
            {
                const auto curPos = this->getRandomPoint();
                if(this->isInside(curPos))
                {
                    ++this->amountInside;
                }
            }
        }
        ValueType GetAmountInside() const
        {
            return this->amountInside;
        }
    private:
        std::array<ValueType, Dimension> getRandomPoint()
        {
            std::array<ValueType, Dimension> res;
            for(std::size_t curDim = 0; curDim < Dimension; ++curDim)
            {
                res[curDim] = this->pointGenerator[curDim](this->randomEngine);
            }
            return res;
        }
};

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 <std::size_t Dimension, typename ValueType>
ValueType getTotalMeasure(const std::array<std::pair<ValueType, ValueType>, Dimension>& limitBox)
{
    auto prodFun = [](const ValueType tmpProd, const std::pair<ValueType, ValueType>& curDim)->ValueType
    {
        return tmpProd*(curDim.second-curDim.first);
    };
    return std::accumulate(limitBox.begin(), limitBox.end(), 1.0L, prodFun);
}

template <std::size_t Dimension, typename ValueType, typename RandomEngineType, typename IsInsideType>
ValueType MeasureMcPar(const std::array<std::pair<ValueType, ValueType>, Dimension>& limitBox,
                       const IsInsideType& isInside,
                       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;
    typedef McWorker<Dimension, ValueType, RandomEngineType, IsInsideType> WorkerType;
    threadVector.reserve(numberOfThreads);
    std::vector<WorkerType> workerVector;
    workerVector.reserve(numberOfThreads);
    const auto mainSeed = time(NULL);
    for(unsigned i = 0; i < numberOfThreads; ++i)
    {
        workerVector.push_back(WorkerType(limitBox, isInside, localIter, i+mainSeed));
        threadVector.push_back(std::thread(std::ref(workerVector[i])));
    }
    for(auto & th : threadVector)
    {
        th.join();
    }

    auto sumFun = [](const std::uintmax_t tmpSum, const WorkerType& inWorker)->std::uintmax_t
    {
        return tmpSum+inWorker.GetAmountInside();
    };
    const std::uintmax_t totalSum = std::accumulate(workerVector.begin(), workerVector.end(), 0ULL, sumFun);
    const auto div = gcd(totalSum, totalIter);
    const auto limitingMeasure = getTotalMeasure<Dimension, ValueType>(limitBox);
    return limitingMeasure*ValueType(totalSum/div)/ValueType(totalIter/div);
}

template <typename ValueType>
class IsInsideEgg
{
    const ValueType r2 = 8.0L*(3.0L-2.0L*sqrt(2.0L));
    public:
        bool operator()(std::array<ValueType, 2> inPos) const
        {
            const auto x = std::abs(inPos[0]);
            const auto y = inPos[1];
            bool res = false;
            if(y <= 0 && x*x+y*y <= 4.0L)
            {
                res = true;
            }
            else if(y >= 0 && y <= x+2 && (x+2)*(x+2)+y*y <= 16.0L)
            {
                res = true;
            }
            else if(y >= x+2 && x*x+(y-2)*(y-2) <= this->r2)
            {
                res = true;
            }
            return res;
        }
};


int main()
{
    typedef long double ValueType;
    typedef IsInsideEgg<ValueType> IsInsideType;
    const IsInsideType isInEgg = IsInsideType();
    const std::array<std::pair<ValueType, ValueType>, 2> boxEgg = {std::make_pair(0.0L, 2.0L),
                                                                   std::make_pair(-2.0L, 4.0L)};

    const unsigned numberOfThreads = std::thread::hardware_concurrency();
    const std::uintmax_t totalIter = std::uintmax_t(numberOfThreads)*50000000ULL;

    const auto startTime = std::chrono::steady_clock::now();
    const auto res = 2*MeasureMcPar<2, ValueType, std::mt19937_64, IsInsideType>(boxEgg,
                                                                                 isInEgg,
                                                                                 totalIter,
                                                                                 numberOfThreads);
    const auto endTime = std::chrono::steady_clock::now();
    std::chrono::duration<double> runTime = endTime-startTime;
    std::cout<<"Runtime: "<<runTime.count()<<" [s]"<<std::endl;
    std::cout<<"result: "<<res<<std::endl;
    return 0;
}