Connect The Dots

Computer Science Level 4

Consider a $100 \times 100$ grid. Suppose four points are randomly chosen on the grid. What is the probability that the convex hull of these four points is a quadrilateral?

Give your answer to 3 decimal places.

For reference on Monte-Carlo, here is an excellent wiki.

Other approaches are definitely encouraged too.

The answer is 0.693.

6 solutions

Abdelhamid Saadi
Aug 4, 2015

Relevant wiki: Monte Carlo Simulation

For four points to satisfy the condition of the problem, every point must be outside the triangle formed by the others points.

This is a Monte Carlo solution written in python 3.4 :

from random import randint

def det(p1, p2, p3):
    "compute determinant of three points"
    return detv((p2[0] - p1[0], p2[1] - p1[1]), (p3[0] - p1[0], p3[1] - p1[1]))

def detv (v1, v2):
    "compute determinant of two vectors"
    return v1[0]*v2[1] - v2[0]*v1[1]

def inside(point, triangle):
    "test if point is inside triangel"
    d1 = det(triangle[0], triangle[1], triangle[2])
    if d1 == 0:
        return True
    else:
        a = det(triangle[0], triangle[1], point)/d1
        b = det(triangle[2], triangle[0], point)/d1
        return a >= 0 and b >= 0 and a + b <= 1

def isconvex(poly):
    "test if the convex hull is quadrilateral"
    for i in range(4):
        if inside(poly[i], poly[0:i] + poly[i+1:]): return False
    return True

def randpoint():
    "genere a random point"
    return (randint(0,100), randint(0,100))

def randquad():
    "genere a randdom quadrilateral"
    return [randpoint() for i in range(4)]

def solve(k):
    "Connect the dots"
    return sum([isconvex(randquad()) for i in range(k)])/k

for i in range(0, 11):
    print(solve(10**i))

Michael Mendrin
May 7, 2016

Determining whether or not $4$ points chosen at random makes for a convex quadrilateral is always a problematic one, but here's a different way to do it:

1) Taking groups of $3$ points out of $4$ , we determine the areas formed by the $4$ triangles, and find the largest area

2) Taking the $3$ different cyclic orders of $4$ points, we determine the $3$ quadrilateral areas formed by the $4$ points, and find the largest area

3) If the largest area in 2) is > than the largest area in 1), then count it as a convex quadrilateral

4) Repeat $100$ M times for a good Monte Carlo run with at least $3$ decimal digits accuracy

Note: Use the shoelace formula for finding areas

Piotr Idzik
Feb 14, 2020

Parallel C++ code doing a suitable simulation:

#include <iostream>
#include <chrono>
#include <thread>
#include <vector>
#include <numeric>

template<typename T>
bool IsInsideTriangle(const T& pX, const T& pY, const T& p0X, const T& p0Y, const T& p1X, const T& p1Y, const T& p2X, const T& p2Y)
{
    //checks if the point (pX, pY) is inside the triangle ((p0X, p0Y), (p1X, p1Y), (p2X, p2Y))
    const T s = p0Y*p2X-p0X*p2Y+(p2Y-p0Y)*pX+(p0X-p2X)*pY;
    const T t = p0X*p1Y-p0Y*p1X+(p0Y-p1Y)*pX+(p1X-p0X)*pY;

    if ((s < 0) != (t < 0))
    {
        return false;
    }

    const T A = -p1Y*p2X+p0Y*(p2X-p1X)+p0X*(p1Y-p2Y)+p1X*p2Y;

    return A < 0 ?
            (s <= 0 && s + t >= A) :
            (s >= 0 && s + t <= A);
}

bool SingleTest(const unsigned xSize, const unsigned ySize)
{
    //chooses randomly 4 points and checks if their convex hull forms a quadrilateral
    int xPosList[4];
    int yPosList[4];
    for(unsigned i = 0; i < 4; ++i)
    {
        xPosList[i] = rand()%xSize;
        yPosList[i] = rand()%ySize;
    }
    return !IsInsideTriangle(xPosList[0], yPosList[0], xPosList[1], yPosList[1], xPosList[2], yPosList[2], xPosList[3], yPosList[3]) &&
           !IsInsideTriangle(xPosList[1], yPosList[1], xPosList[0], yPosList[0], xPosList[2], yPosList[2], xPosList[3], yPosList[3]) &&
           !IsInsideTriangle(xPosList[2], yPosList[2], xPosList[1], yPosList[1], xPosList[0], yPosList[0], xPosList[3], yPosList[3]) &&
           !IsInsideTriangle(xPosList[3], yPosList[3], xPosList[1], yPosList[1], xPosList[2], yPosList[2], xPosList[0], yPosList[0]);
}

long double SingleRun(const unsigned xSize, const unsigned ySize, const unsigned totalIter)
{
    unsigned numberOfQuadrilateral = 0;
    for(unsigned i = 0; i<totalIter; ++i)
    {
        if(SingleTest(xSize, ySize))
        {
            ++numberOfQuadrilateral;
        }
    }
    return static_cast<long double>(numberOfQuadrilateral)/static_cast<long double>(totalIter);
}

class Worker
{
    private:
        unsigned xSize, ySize, localIter;
        long double resultValue;
    public:
        Worker(const unsigned _xSize, const unsigned _ySize, const unsigned _localIter) :
        xSize(_xSize),
        ySize(_ySize),
        localIter(_localIter),
        resultValue(-1.0)
        {}
        void operator()()
        {
            this->resultValue = SingleRun(this->xSize, this->ySize, this->localIter);
        }
        long double GetResultValue() const
        {
            return this->resultValue;
        }
};

long double mcQuadrilateralPar(const unsigned xSize, const unsigned ySize, const unsigned totalIter, const unsigned numberOfThreads)
{
    const unsigned localIter = totalIter/numberOfThreads;
    std::vector<std::thread> threadVector;
    threadVector.reserve(numberOfThreads);
    std::vector<Worker> workerVector;
    workerVector.reserve(numberOfThreads);
    for(unsigned i = 0; i < numberOfThreads; ++i)
    {
        workerVector.push_back(Worker(xSize, ySize, localIter));
        threadVector.push_back(std::thread(std::ref(workerVector[i])));
    }
    for(auto & th : threadVector)
    {
        th.join();
    }

    auto sumFun = [](const long double tmpSum, const Worker& inWorker)
    {
        return tmpSum+inWorker.GetResultValue();
    };
    long double sum = std::accumulate(workerVector.begin(), workerVector.end(), 0.0, sumFun);
    return sum/static_cast<long double>(numberOfThreads);
}

int main()
{
    const unsigned xSize = 100;
    const unsigned ySize = xSize;
    const unsigned totalIter = 20000000000;
    const unsigned threadNum = std::thread::hardware_concurrency()/2 > 0 ? std::thread::hardware_concurrency()/2 : 1;
    auto startTime = std::chrono::steady_clock::now();
    std::cout<<"Result: "<<mcQuadrilateralPar(xSize, ySize, totalIter, threadNum)<<std::endl;
    auto endTime = std::chrono::steady_clock::now();
    std::chrono::duration<double> runTime = endTime-startTime;
    std::cout<<"Runtime: "<<runTime.count()<<" [s]"<<std::endl;

    return 0;
}

Sourajit Debnath
Aug 26, 2016

Compute the average area of a triangle when any three of the points are i.i.d. chosen from a square. The fourth point must lie within the boundary of this triangle to not be part of the convex hull. This point can be chosen in four ways. The complement of this probability is the solution.

This is the solution in MATLAB.

    side = 100;
    n = 10000000;
    data = side * rand(n, 6);
    op = data(:, 1) .* data(:, 4) - data(:, 2) .* data(:, 3) + ...
         data(:, 3) .* data(:, 6) - data(:, 4) .* data(:, 5) + ...
         data(:, 5) .* data(:, 2) - data(:, 6) .* data(:, 1);
    op = 0.5 * abs(op) / (side ^ 2);
    disp(1 - 4 * mean(op))

Christopher Boo
May 23, 2016

Relevant wiki: Monte-Carlo Simulation

I always liked Monte-Carlo simulation. Everything is so mysterious until you run it and see how the probability converges.

For this problem, if the convex hull formed by the 4 points are not a quadrilateral, this implies that one of the points must be within the triangle formed by the remaining three points. Fortunately there are only 4 points and thus we don't have to code the convex hull algorithm. Instead, we can check if one point is within the triangle. To do that, we first compute the areas of the 4 possible triangles (with the shoelace formula ) and see if any one of them is the sum of the remaining three. If this holds, the convex hull is a triangle and is not a quadrilateral.

from __future__ import division
from random import randint

def area(x1,y1,x2,y2,x3,y3):
    # did not divide by 2 to avoid precision errors, but the equation still holds
    return abs(x1 * y2 + x2 * y3 + x3 * y1 - x3 * y2 - x1 * y3 - x2 * y1)

def in_triangle(x1,y1,x2,y2,x3,y3,x4,y4):
    A = area(x1,y1,x2,y2,x3,y3)
    B = area(x1,y1,x2,y2,x4,y4)
    C = area(x1,y1,x3,y3,x4,y4)
    D = area(x2,y2,x3,y3,x4,y4)

    # if a point is within the bigger triangle
    T = max(A,B,C,D)

    return A + B + C + D - T == T

cnt = 0    # cases where the convex hull is a quadrilateral
tst = 0    # total number of cases
while True:
    tst += 1

    x1 = randint(0,100)
    x2 = randint(0,100)
    x3 = randint(0,100)
    x4 = randint(0,100)
    y1 = randint(0,100)
    y2 = randint(0,100)
    y3 = randint(0,100)
    y4 = randint(0,100)

    if not in_triangle(x1,y1,x2,y2,x3,y3,x4,y4):
        cnt += 1

    print tst, cnt / tst

With this code you can see how the probability converges to 0.6925 .

Bill Bell
Aug 9, 2015

Graham's scan is one of several ways of computing the convex hull of a set of points in the plane. This codes uses a version due to Mr Tom Switzer.

# Graham Scan - Tom Switzer <thomas.switzer@gmail.com>

TURN_LEFT, TURN_RIGHT, TURN_NONE = (1, -1, 0)

def turn(p, q, r):
    return cmp((q[0] - p[0])*(r[1] - p[1]) - (r[0] - p[0])*(q[1] - p[1]), 0)

def _keep_left(hull, r):
    while len(hull) > 1 and turn(hull[-2], hull[-1], r) != TURN_LEFT:
            hull.pop()
    if not len(hull) or hull[-1] != r:
        hull.append(r)
    return hull

def convex_hull(points):
    points = sorted(points)
    l = reduce(_keep_left, points, [])
    u = reduce(_keep_left, reversed(points), [])
    return l.extend(u[i] for i in xrange(1, len(u) - 1)) or l

from random import randint

N=100000
count=0
for n in range(N):
    points=[(randint(1,100),randint(1,100)) for i in range(4)]
    count+=len(convex_hull(points))==4
print 1.*count/N

Connect The Dots

Other approaches are definitely encouraged too.

The answer is 0.693.

6 solutions

0 pending reports