python Decimal() with extreme precision acting funky

Question

I am trying to evaluate a rather long, complicated mathematical expression, and am getting some very weird results out of the python decimal package. (I already tried using regular floats and they don't have sufficient precision). Here is the code (it's a bit long, sorry):

from decimal import *

def docalc (xi0,yi0,xi1,yi1,xi2,yi2,xi3,yi3,xo0,yo0,xo1,yo1,xo2,yo2,xo3,yo3):
    res = (((xi2 * (-yi0 + yi1) + xi1 * (yi0 - yi2) + xi0 * (-yi1 + yi2)) *
          ((xi2 * (-xo0 + xo1) * yi0 * yi1 + xi1 * (xo0 - xo2) * yi0 * yi2 +
            xi0 * (-xo1 + xo2) * yi1 * yi2) *
           (xi2 * xi3 * (xo2 - xo3) * yi0 + xi0 * xi3 * (-xo0 + xo3) * yi2 +
            xi0 * xi2 * (xo0 - xo2) * yi3) -
           (xi1 * xi2 * (-xo1 + xo2) * yi0 + xi0 * xi2 * (xo0 - xo2) * yi1 +
            xi0 * xi1 * (-xo0 + xo1) * yi2) *
           (xi3 * (xo0 - xo2) * yi0 * yi2 + xi2 * (-xo0 + xo3) * yi0 * yi3 +
            xi0 * (xo2 - xo3) * yi2 * yi3)) *
          (xi1 * xi2 * (yi0 - yi3) * (yo1 - yo2) +
           xi0 * (xi1 * (yi2 - yi3) * (yo0 - yo1) - xi2 * (yi1 - yi3) * (yo0 - yo2) +
                  xi3 * (yi1 - yi2) * (yo0 - yo3)) - xi1 * xi3 * (yi0 - yi2) * (yo1 - yo3) +
           xi2 * xi3 * (yi0 - yi1) * (yo2 - yo3)) +
          (-(xi2 * yi0) + xi0 * yi2) * (xi2 * xi3 * (xo2 - xo3) * (yi0 - yi1) +
                                        xi0 * xi3 * (xo0 - xo3) * (yi1 - yi2) +
                                        xi1 * (-(xi3 * (xo1 - xo3) * (yi0 - yi2)) + xi2 * (xo1 - xo2) * (yi0 - yi3) +
                                               xi0 * (xo0 - xo1) * (yi2 - yi3)) - xi0 * xi2 * (xo0 - xo2) * (
                                        yi1 - yi3)) *
          ((-(xi2 * (-xo0 + xo1) * yi0 * yi1) - xi1 * (xo0 - xo2) * yi0 * yi2 -
            xi0 * (-xo1 + xo2) * yi1 * yi2) *
           (xi1 * xi2 * (yi0 - yi3) * (yo1 - yo2) +
            xi0 * (xi1 * (yi2 - yi3) * (yo0 - yo1) - xi2 * (yi1 - yi3) * (yo0 - yo2) +
                   xi3 * (yi1 - yi2) * (yo0 - yo3)) -
            xi1 * xi3 * (yi0 - yi2) * (yo1 - yo3) + xi2 * xi3 * (yi0 - yi1) * (yo2 - yo3))
           - (xi1 * xi2 * (-xo1 + xo2) * yi0 + xi0 * xi2 * (xo0 - xo2) * yi1 +
              xi0 * xi1 * (-xo0 + xo1) * yi2) *
           (xi1 * yi0 * yi2 * yo0 - xi1 * yi0 * yi3 * yo0 + xi3 * yi0 * yi1 * (yo0 - yo1) -
            xi0 * yi1 * yi2 * yo1 + xi0 * yi1 * yi3 * yo1 + xi3 * yi1 * yi2 * (yo1 - yo2) -
            xi1 * yi0 * yi2 * yo2 + xi0 * yi1 * yi2 * yo2 - xi0 * yi2 * yi3 * yo2 +
            xi1 * yi2 * yi3 * yo2 + xi3 * yi0 * yi2 * (-yo0 + yo2) +
            xi2 * yi0 * (yi1 * (-yo0 + yo1) + yi3 * (yo0 - yo3)) + xi1 * yi0 * yi3 * yo3 -
            xi0 * yi1 * yi3 * yo3 + xi0 * yi2 * yi3 * yo3 - xi1 * yi2 * yi3 * yo3 +
            xi2 * yi1 * yi3 * (-yo1 + yo3)))))

    return (res);


def test():

    xi0 = 1.1;      yi0 = 1.1;      xi1 = 2.1;      yi1 = 1.1;      xi2 = 2.1;      yi2 = 2.1;      xi3 = 1.1;      yi3 = 2.1
    xo0 = 91.23;    yo0 = 82.34;    xo1 = 73.45;    yo1 = 64.56;    xo2 = 55.67;    yo2 = 46.78;    xo3 = 37.89;    yo3 = 28.90
    c_float = docalc(xi0, yi0, xi1, yi1, xi2, yi2, xi3, yi3, xo0, yo0, xo1, yo1, xo2, yo2, xo3, yo3)

    getcontext().prec = 19
    xi0d = Decimal(xi0);   yi0d = Decimal(yi0);    xi1d = Decimal(xi1);   yi1d = Decimal(yi1);    xi2d = Decimal(xi2);   yi2d = Decimal(yi2);    xi3d = Decimal(xi3);   yi3d = Decimal(yi3);
    xo0d = Decimal(xo0);   yo0d = Decimal(yo0);    xo1d = Decimal(xo1);   yo1d = Decimal(yo1);    xo2d = Decimal(xo2);   yo2d = Decimal(yo2);    xo3d = Decimal(xo3);   yo3d = Decimal(yo3);
    c_decimal_19 = docalc(xi0d, yi0d, xi1d, yi1d, xi2d, yi2d, xi3d, yi3d, xo0d, yo0d, xo1d, yo1d, xo2d, yo2d, xo3d, yo3d)

    getcontext().prec = 500
    xi0d = Decimal(xi0);   yi0d = Decimal(yi0);    xi1d = Decimal(xi1);   yi1d = Decimal(yi1);    xi2d = Decimal(xi2);   yi2d = Decimal(yi2);    xi3d = Decimal(xi3);   yi3d = Decimal(yi3);
    xo0d = Decimal(xo0);   yo0d = Decimal(yo0);    xo1d = Decimal(xo1);   yo1d = Decimal(yo1);    xo2d = Decimal(xo2);   yo2d = Decimal(yo2);    xo3d = Decimal(xo3);   yo3d = Decimal(yo3);
    c_decimal_500 = docalc(xi0d, yi0d, xi1d, yi1d, xi2d, yi2d, xi3d, yi3d, xo0d, yo0d, xo1d, yo1d, xo2d, yo2d, xo3d, yo3d)

    print("using floats: ",c_float);
    print("using Decimals() w/ prec=19: ",c_decimal_19);
    print("using Decimals() w/ prec=500: ",c_decimal_500);

test();

Here is the output I get:

C:\dev\Python36-32\python.exe C:/dev/PythonProjects/untitled4/main.py
using floats:  -1.1703014024533338e-10
using Decimals() w/ prec=19:  -4.289200000000001210E-13
using Decimals() w/ prec=500:  0E-700

Process finished with exit code 0

For the record, Mathematica says the actual answer is 7.80201*10^-11. I'm not surprised that using floats with such a long complicated expression results in horrendous rounding errors, but I am confused that the Decimal() type, (even with precision set to 500), returns weird spurious results. In fact, increasing the precision from 19 to 500 seems to cause the result to get worse. Any ideas? Is this a bug in the python decimals package, or maybe a bug in python in general with extremely long, convoluted mathematical expressions? Thanks.