| 1708 | } |
| 1709 | |
| 1710 | REAL orient3dexact(pa, pb, pc, pd) |
| 1711 | REAL *pa; |
| 1712 | REAL *pb; |
| 1713 | REAL *pc; |
| 1714 | REAL *pd; |
| 1715 | { |
| 1716 | INEXACT REAL axby1, bxcy1, cxdy1, dxay1, axcy1, bxdy1; |
| 1717 | INEXACT REAL bxay1, cxby1, dxcy1, axdy1, cxay1, dxby1; |
| 1718 | REAL axby0, bxcy0, cxdy0, dxay0, axcy0, bxdy0; |
| 1719 | REAL bxay0, cxby0, dxcy0, axdy0, cxay0, dxby0; |
| 1720 | REAL ab[4], bc[4], cd[4], da[4], ac[4], bd[4]; |
| 1721 | REAL temp8[8]; |
| 1722 | int templen; |
| 1723 | REAL abc[12], bcd[12], cda[12], dab[12]; |
| 1724 | int abclen, bcdlen, cdalen, dablen; |
| 1725 | REAL adet[24], bdet[24], cdet[24], ddet[24]; |
| 1726 | int alen, blen, clen, dlen; |
| 1727 | REAL abdet[48], cddet[48]; |
| 1728 | int ablen, cdlen; |
| 1729 | REAL deter[96]; |
| 1730 | int deterlen; |
| 1731 | int i; |
| 1732 | |
| 1733 | INEXACT REAL bvirt; |
| 1734 | REAL avirt, bround, around; |
| 1735 | INEXACT REAL c; |
| 1736 | INEXACT REAL abig; |
| 1737 | REAL ahi, alo, bhi, blo; |
| 1738 | REAL err1, err2, err3; |
| 1739 | INEXACT REAL _i, _j; |
| 1740 | REAL _0; |
| 1741 | |
| 1742 | Two_Product(pa[0], pb[1], axby1, axby0); |
| 1743 | Two_Product(pb[0], pa[1], bxay1, bxay0); |
| 1744 | Two_Two_Diff(axby1, axby0, bxay1, bxay0, ab[3], ab[2], ab[1], ab[0]); |
| 1745 | |
| 1746 | Two_Product(pb[0], pc[1], bxcy1, bxcy0); |
| 1747 | Two_Product(pc[0], pb[1], cxby1, cxby0); |
| 1748 | Two_Two_Diff(bxcy1, bxcy0, cxby1, cxby0, bc[3], bc[2], bc[1], bc[0]); |
| 1749 | |
| 1750 | Two_Product(pc[0], pd[1], cxdy1, cxdy0); |
| 1751 | Two_Product(pd[0], pc[1], dxcy1, dxcy0); |
| 1752 | Two_Two_Diff(cxdy1, cxdy0, dxcy1, dxcy0, cd[3], cd[2], cd[1], cd[0]); |
| 1753 | |
| 1754 | Two_Product(pd[0], pa[1], dxay1, dxay0); |
| 1755 | Two_Product(pa[0], pd[1], axdy1, axdy0); |
| 1756 | Two_Two_Diff(dxay1, dxay0, axdy1, axdy0, da[3], da[2], da[1], da[0]); |
| 1757 | |
| 1758 | Two_Product(pa[0], pc[1], axcy1, axcy0); |
| 1759 | Two_Product(pc[0], pa[1], cxay1, cxay0); |
| 1760 | Two_Two_Diff(axcy1, axcy0, cxay1, cxay0, ac[3], ac[2], ac[1], ac[0]); |
| 1761 | |
| 1762 | Two_Product(pb[0], pd[1], bxdy1, bxdy0); |
| 1763 | Two_Product(pd[0], pb[1], dxby1, dxby0); |
| 1764 | Two_Two_Diff(bxdy1, bxdy0, dxby1, dxby0, bd[3], bd[2], bd[1], bd[0]); |
| 1765 | |
| 1766 | templen = fast_expansion_sum_zeroelim(4, cd, 4, da, temp8); |
| 1767 | cdalen = fast_expansion_sum_zeroelim(templen, temp8, 4, ac, cda); |