Numerical precision of implementation for convex quadrilateral area

Question

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 < 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 <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?

Answer 1

Doubles have 52 bits of fraction. In base ten, that is about 15.7 digits. Your numbers have eight digits on the left side of the dot, so you can only expect to have seven digits correct on the right side of the dot - and that is just due to the representation, without even considering any cumulative errors produced by the calculations.

So, the answer is: No, there is no way. You have to use a format with more precision.