Numerical precision of implementation for convex quadrilateral area I have implemented a method to compute a convex quadrilateral area in R3. The method works fine, but I am having numerical precision problems in the 8th decimal place. Take a look on the method:

internal static double GetTriangleArea(double ax, double ay, double az,
            double bx, double by, double bz,
            double cx, double cy, double cz
            )
        {
            /**
             * AB = B-A = (ux, uy, uz)
             * AC = C-A = (wx, wy, wz)
             * 
             * S = 0.5*sqrt{(uy*wz - uz*wy)² + (uz*wx - ux*wz)² + (ux*wy - uy*wx)²}
             * */

            var ux = bx - ax;
            var uy = by - ay;
            var uz = bz - az;

            var wx = cx - ax;
            var wy = cy - ay;
            var wz = cz - az;

            var t1 = uy*wz - uz*wy;
            var t2 = uz*wx - ux*wz;
            var t3 = ux*wy - uy*wx;
            var s = 0.5*Math.Sqrt(t1*t1 + t2*t2 + t3*t3);

            return s;
        }

        internal static double GetConvexQuadrilateralArea(double ax, double ay, double az,
            double bx, double by, double bz,
            double cx, double cy, double cz,
            double dx, double dy, double dz)
        {
            var triangle1 = GetTriangleArea(ax, ay, az, bx, by, bz, cx, cy, cz);
            var triangle2 = GetTriangleArea(ax, ay, az, cx, cy, cz, dx,dy,dz);

            return triangle1 + triangle2;
        }


And this is the test:

[TestMethod]
        public void ParallelogramOfBaseBAndHeightHMustHaveAreaEqualToBTimesH()
        {
            var random = new Random(1);
            const double scale = 10000;
            for (var counter = 0; counter &lt; 1000; counter++)
            {
                double baseLength = random.NextDouble() * scale;
                double height = random.NextDouble() * scale;

                double dx = random.NextDouble()*scale;

                var a = new[] { 0, 0, 0 };
                var b = new[] { baseLength, 0, 0 };
                var c = new[] { baseLength+dx, height, 0 };
                var d = new[] { 0F+dx, height, 0 };

                double expected = baseLength * height;

                var result = MathUtils.GetConvexQuadrilateralArea(a[0], a[1], a[2], b[0], b[1], b[2], c[0], c[1], c[2],
                    d[0], d[1], d[2]);

                Assert.AreEqual(expected, result, Epsilon*scale,
                    string.Format("sideA: {0}, height: {1}, dx: {2}", baseLength, height, dx));
            }
        }


This test fails with the following message: Expected a difference no greater than &lt;1E-09> between expected value 74813926.2967871 and actual value 74813926.2967871. sideA: 8552.44307245707, height: 8747.66726454146, dx: 4721.64729829954.

My question is: is there a way to increase the numerical precision of my implementation while still using double precision numbers?
Debugging floating-point precision issues is, as @Bobby mentioned, off-topic for this site, so this is going to be just a code review :)

The first thing that strikes me, is the cryptic two-letter identifiers.

Let's start with the signature for GetTriangleArea. You're taking 9 parameters, this is bad. In reality you really need 3 groups of 3 numbers. Your code has failed to make this abstraction, and this is why you end up with a1-a2-a3 -style identifiers.

What if GetTriangleArea had a signature like this?

internal static double GetTriangleArea(Point3D pointA, Point3D pointB, Point3D pointC)


Given a Point3D struct with members X, Y and Z, the rest of the code would be much less cryptic.

The same applies to GetConvexQuadrilateralArea, which could look much less crowded like this:

internal static double GetConvexQuadrilateralArea(Point3D pointA, Point3D pointB, Point3D pointC, Point3D pointD)


Lastly, I'm not sure, but I think GetConvexQuadrilateralArea is lying - the method isn't actively verifying the convex-ness of the points it's given, which may or may not be problematic.