# Monte-Carlo method to check triangle inequality

This is code for finding if inequality $$\abc > 2(\text{area of triangle})\$$ where $$\a\$$, $$\b\$$ and $$\c\$$ are the sides of any triangle.

The below, and repl, code work as intended. However they take a long time. Can the performance be improved without reducing max_tries?

import random
import math
hn=100000000
num=0
max_tries=10000000
tries = 0
while tries != max_tries:
a= random.randint(1,hn)
b= random.randint(1,hn)
c= random.randint(1,hn)
if a+b>c and c+b>a and a+c>b :
s=(a+b+c)/2
tries= tries +1
k =math.sqrt(s*(s-a)*(s-b)*(s-c))
if a*b*c> 2*k:
num= num+1

print(num/max_tries)

• Can you flesh out the description on what this does. I've read your post but you've mostly just complained about speed without explaining what you're tying to do. Mar 27, 2020 at 12:30
• Also, out of curiosity, did you ever find an example where it is not true or do you expect to? Mar 27, 2020 at 12:34
• My expectation what triangle inequality is about is #1 (the if condition above) - abc does show up in at least one of the inequalities for triangle area. Mar 27, 2020 at 14:17
• It would be interesting to see the impact of replacing just one of the lengths in each iteration. Mar 27, 2020 at 14:20
• Why can the lengths only be integers? If you want integers to prevent problems with floating point rounding, you might try harder to keep s and k integers too by slightly rewriting the inequality Mar 27, 2020 at 15:34

The problem of your code is the way you generate triangles. You just create three independent random numbers. In theory, you could always end up rolling numbers that don't end up in a valid triangle and thus your code would not even terminate. That is, because triangle side lengths are dependent.

To always create a valid triangle, you can instead roll two angles and one side and then use the law of sines to calculate the remaining sides.

The sum of angles in a triangle is 180 degrees (or pi radians). So the first angle alpha is randomly selected between (0, pi), the second angle beta is then between (0, pi-alpha) and the third angle gamma is gained by simple arithmetics: gamma = pi - alpha - beta

Or in code:

alpha = random.uniform(0, math.pi)
beta  = random.uniform(0, math.pi - alpha)
gamma = math.pi - alpha - beta


If we know two sides and the angle in between, we can easily calculate the area. Since your inequality uses twice the area, we can drop the factor of 0.5.

double_area = side_b * side_c * sin_alpha


Some micro-optimisations: math.sin is an expensive function (time-wise), so we save the result in a variable. Same goes for floating point division. Multiplication is faster, so once again, we save the result to use multiplication later on.

Lastly: Since you already know how often you want to execute the loop, you can use a for-loop instead of a while-loop although frankly, I don't know how well this goes in terms of memory consumption.

import math
import random
import time

max_side_length = 100000000
n_max_tries     = 10000000
n_successes     = 0

time_start      = time.time()

for k in range(n_max_tries):
alpha                   = random.uniform(0, math.pi)
beta                    = random.uniform(0, math.pi - alpha)
gamma                   = math.pi - alpha - beta

side_a                  = random.randint(1, max_side_length)
sin_alpha               = math.sin(alpha)
side_a_over_sin_alpha   = side_a / sin_alpha

side_b                  = side_a_over_sin_alpha * math.sin(beta)
side_c                  = side_a_over_sin_alpha * math.sin(gamma)

double_area             = side_b * side_c * sin_alpha

if side_a * side_b * side_c > double_area:
n_successes += 1

time_end = time.time()

print("successes: ", n_successes, " / ", n_max_tries, " (", n_successes/n_max_tries, ")")
print("time: ", time_end-time_start)


I didn't do excessive timings because I couldn't be bothered, but this code runs on my machine in about 37.8s, while three runs of your code averaged on 109.8s, resulting in an optimistic speed-up factor of around 3.

Disclaimer: I can somewhat write in python but I'm probably violating this famous pep-8. Yet, I strongly encourage you to give reasonable names to your variables as it makes the code more easy to understand. Not only for others, but also for yourself, when you look at it 2 weeks later.

• Thanks it is working much faster Mar 27, 2020 at 13:48
• Wouldn't it be easier to say double_area = side_a**2 * math.sin(alpha) * math.sin(beta) * math.sin(gamma)? Mar 27, 2020 at 14:34
• Alternative approach to generating a valid length triple from every pair: draw c from [abs(a-b), a+b]. Apr 11, 2020 at 5:39

One way to speed this up is to use a package designed for numerical evaluation, numpy. It is implemented in C and can take advantage of multiple cores if necessary. The only limitation is that you can only create arrays that fit into memory. At least on my machine (16GB RAM) your current values fit easily, though (if each integer takes 64 bytes, then I could store over 250 million of them).

Another limitation is that you are limited to long long integers (64 bytes on my machine), so you will have to adapt your upper ceiling for the numbers to that. The largest value being produced by this code is s * (s - a) * (s - b) * (s - c), so this needs to fit even if all of a, b and c are the maximum value. There is a programmatic way to get the size, so we can figure out what it is. Unfortunately, this formula has basically a $$\n^4\$$ in it, so the maximum value is quite small. On my machine I can generate values of up to just 83748 this way. Although practically this is probably not relevant, since the likelihood that you generate a larger number like that is very small (case in point, I had previsously assumed a * b * c /2 would be the largest value, but I still found 100% matching the inequality). The only problem with this is that if you do get a larger value, you get an overflow, the next number after the largest integer is the smallest integer, so you get negative numbers under the square root. But numpy will tell you about that with a warning, if it happens.

Here is a straight-forward implementation of your code using numpy. I did not make any effort to speed it up in addition, the only real difference is the value of limit and that I moved the factor 2 in the inequality to the other side (so I can make the limit a factor two larger).

import numpy as np

#limit = 100000000
max_int = np.iinfo(np.longlong).max
limit = int((16./3 * max_int) ** (1./4)) + 1   # 83748 with np.int64
max_tries = 10000000

a, b, c = np.random.randint(1, limit, size=(3, max_tries), dtype=np.longlong)
mask = (a + b > c) & (a + c > b) & (b + c > a)
s = (a + b + c) / 2
k = np.sqrt(s * (s - a) * (s - b) * (s - c))
inequality_true = a / 2 * b * c > k

print(inequality_true.sum() / len(inequality_true))


This script runs in about 1.1s on my machine. Your code (with the same limit for the random numbers) takes about 83s.

In addition to speeding up your code, you should also try to make it more readable. Python has an official style-guide, PEP8. It recommends using spaces around the = when doing assignments and around binary operators (like + and *). In addition, try to use clear variable names (hn tells me nothing about what it is for) and put your main calling code under a if __name__ == "__main__": guard if you want to be able to import from a script without everything running.

• I'd be curious about my solution in a numpy approach. Mar 27, 2020 at 14:21
• @infinitezero: It takes about 1.6s. Without doing the math for the maximum limit like I did for my script, which would probably speed it up more since the numbers are smaller. Although your script allows larger numbers, since a * b * c is the largest number there. Mar 27, 2020 at 14:38
• In numpy >= 1.18, randint() can take an array as the lower or upper bound. So I think you can do a, b = np.random.randint(1, limit, size=(2, max_tries), dtype=np.longlong) and then c = np.random.randint(abs(a-b+1), a+b, size=max_tries, dtype=np.longlong). That way each triple is guaranteed to be the sides of a triangle. Mar 29, 2020 at 1:51