diff --git a/common/mathlib.cc b/common/mathlib.cc
index 55e330b7..89d5efd7 100644
--- a/common/mathlib.cc
+++ b/common/mathlib.cc
@@ -613,11 +613,10 @@ qvec3f GLM_PolyCentroid(const std::vector<qvec3f> &points)
     return poly_centroid;
 }
 
-qvec3f GLM_PolyRandomPoint(const std::vector<qvec3f> &points)
+poly_random_point_state_t GLM_PolyRandomPoint_Setup(const std::vector<qvec3f> &points)
 {
     Q_assert(points.size() >= 3);
-    
-    // FIXME: Precompute this
+
     float poly_area = 0;
     std::vector<float> triareas;
     
@@ -632,18 +631,28 @@ qvec3f GLM_PolyRandomPoint(const std::vector<qvec3f> &points)
         triareas.push_back(triarea);
         poly_area += triarea;
     }
-    
-    // Pick a random triangle, with probability proportional to triangle area
-    const float uniformRandom = Random();
     const std::vector<float> cdf = MakeCDF(triareas);
-    const int whichTri = SampleCDF(cdf, uniformRandom);
     
-    Q_assert(whichTri >= 0 && whichTri < triareas.size());
+    poly_random_point_state_t result;
+    result.points = points;
+    result.triareas = triareas;
+    result.triareas_cdf = cdf;
+    return result;
+}
 
-    const tri_t tri { points.at(0), points.at(1 + whichTri), points.at(2 + whichTri) };
+// r1, r2, r3 must be in [0, 1]
+qvec3f GLM_PolyRandomPoint(const poly_random_point_state_t &state, float r1, float r2, float r3)
+{
+    // Pick a random triangle, with probability proportional to triangle area
+    const float uniformRandom = r1;
+    const int whichTri = SampleCDF(state.triareas_cdf, uniformRandom);
+    
+    Q_assert(whichTri >= 0 && whichTri < state.triareas.size());
+
+    const tri_t tri { state.points.at(0), state.points.at(1 + whichTri), state.points.at(2 + whichTri) };
     
     // Pick random barycentric coords.
-    const qvec3f bary = Barycentric_Random(Random(), Random());
+    const qvec3f bary = Barycentric_Random(r2, r3);
     const qvec3f point = Barycentric_ToPoint(bary, tri);
     
     return point;
diff --git a/include/common/mathlib.hh b/include/common/mathlib.hh
index b55a0e12..ededf084 100644
--- a/include/common/mathlib.hh
+++ b/include/common/mathlib.hh
@@ -336,7 +336,16 @@ std::pair<bool, qvec3f> GLM_InterpolateNormal(const std::vector<qvec3f> &points,
 std::vector<qvec3f> GLM_ShrinkPoly(const std::vector<qvec3f> &poly, const float amount);
 /// Returns (front part, back part)
 std::pair<std::vector<qvec3f>,std::vector<qvec3f>> GLM_ClipPoly(const std::vector<qvec3f> &poly, const qvec4f &plane);
-qvec3f GLM_PolyRandomPoint(const std::vector<qvec3f> &points);
+
+class poly_random_point_state_t {
+public:
+    std::vector<qvec3f> points;
+    std::vector<float> triareas;
+    std::vector<float> triareas_cdf;
+};
+
+poly_random_point_state_t GLM_PolyRandomPoint_Setup(const std::vector<qvec3f> &points);
+qvec3f GLM_PolyRandomPoint(const poly_random_point_state_t &state, float r1, float r2, float r3);
 
 /// projects p onto the vw line.
 /// returns 0 for p==v, 1 for p==w
diff --git a/light/test_light.cc b/light/test_light.cc
index 150928de..a56c42b4 100644
--- a/light/test_light.cc
+++ b/light/test_light.cc
@@ -514,9 +514,11 @@ TEST(mathlib, RandomPointInPoly) {
     qvec3f max(-FLT_MAX);
     qvec3f avg(0);
     
+    const auto randomstate = GLM_PolyRandomPoint_Setup(poly);
+    
     const int N=100;
     for (int i=0; i<N; i++) {
-        const qvec3f point = GLM_PolyRandomPoint(poly);
+        const qvec3f point = GLM_PolyRandomPoint(randomstate, Random(), Random(), Random());
         ASSERT_TRUE(GLM_EdgePlanes_PointInside(edgeplanes, point));
         
         //std::cout << "point: " << qv::to_string(point) << std::endl;