diff --git a/src/lib.rs b/src/lib.rs index 52894e1..e5c147d 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -59,6 +59,9 @@ const O3DERRBOUND_C: f64 = (26.0 + 288.0 * EPSILON) * EPSILON * EPSILON; const ICCERRBOUND_A: f64 = (10.0 + 96.0 * EPSILON) * EPSILON; const ICCERRBOUND_B: f64 = (4.0 + 48.0 * EPSILON) * EPSILON; const ICCERRBOUND_C: f64 = (44.0 + 576.0 * EPSILON) * EPSILON * EPSILON; +const ISPERRBOUND_A: f64 = (16.0 + 224.0 * EPSILON) * EPSILON; +const ISPERRBOUND_B: f64 = (5.0 + 72.0 * EPSILON) * EPSILON; +const ISPERRBOUND_C: f64 = (71.0 + 1408.0 * EPSILON) * EPSILON * EPSILON; /// Returns a positive value if the coordinates `pa`, `pb`, and `pc` occur in counterclockwise order /// (pc lies to the **left** of the directed line defined by coordinates pa and pb). @@ -1322,6 +1325,554 @@ fn incircleadapt( fin1[finlength - 1] } +/// Return a positive value if the point pe lies inside the sphere passing through `pa`, `pb`, `pc`, and `pd` +/// Returns a negative value if it lies outside +/// Returns `0` if the five points are **cospherical** +/// The points pa, pb, pc, and pd must be ordered so that they have a positive orientation +pub fn insphere>( + pa: Coord3D, + pb: Coord3D, + pc: Coord3D, + pd: Coord3D, + pe: Coord3D, +) -> f64 { + let pa = Coord3D { + x: pa.x.into(), + y: pa.y.into(), + z: pa.z.into(), + }; + let pb = Coord3D { + x: pb.x.into(), + y: pb.y.into(), + z: pb.z.into(), + }; + let pc = Coord3D { + x: pc.x.into(), + y: pc.y.into(), + z: pc.z.into(), + }; + let pd = Coord3D { + x: pd.x.into(), + y: pd.y.into(), + z: pd.z.into(), + }; + let pe = Coord3D { + x: pe.x.into(), + y: pe.y.into(), + z: pe.z.into(), + }; + + let aex = pa.x - pe.x; + let bex = pb.x - pe.x; + let cex = pc.x - pe.x; + let dex = pd.x - pe.x; + let aey = pa.y - pe.y; + let bey = pb.y - pe.y; + let cey = pc.y - pe.y; + let dey = pd.y - pe.y; + let aez = pa.z - pe.z; + let bez = pb.z - pe.z; + let cez = pc.z - pe.z; + let dez = pd.z - pe.z; + + let aexbey = aex * bey; + let bexaey = bex * aey; + let ab = aexbey - bexaey; + let bexcey = bex * cey; + let cexbey = cex * bey; + let bc = bexcey - cexbey; + let cexdey = cex * dey; + let dexcey = dex * cey; + let cd = cexdey - dexcey; + let dexaey = dex * aey; + let aexdey = aex * dey; + let da = dexaey - aexdey; + + let aexcey = aex * cey; + let cexaey = cex * aey; + let ac = aexcey - cexaey; + let bexdey = bex * dey; + let dexbey = dex * bey; + let bd = bexdey - dexbey; + + let abc = aez * bc - bez * ac + cez * ab; + let bcd = bez * cd - cez * bd + dez * bc; + let cda = cez * da + dez * ac + aez * cd; + let dab = dez * ab + aez * bd + bez * da; + + let alift = aex * aex + aey * aey + aez * aez; + let blift = bex * bex + bey * bey + bez * bez; + let clift = cex * cex + cey * cey + cez * cez; + let dlift = dex * dex + dey * dey + dez * dez; + + let det = (dlift * abc - clift * dab) + (blift * cda - alift * bcd); + + let aezplus = abs(aez); + let bezplus = abs(bez); + let cezplus = abs(cez); + let dezplus = abs(dez); + let aexbeyplus = abs(aexbey); + let bexaeyplus = abs(bexaey); + let bexceyplus = abs(bexcey); + let cexbeyplus = abs(cexbey); + let cexdeyplus = abs(cexdey); + let dexceyplus = abs(dexcey); + let dexaeyplus = abs(dexaey); + let aexdeyplus = abs(aexdey); + let aexceyplus = abs(aexcey); + let cexaeyplus = abs(cexaey); + let bexdeyplus = abs(bexdey); + let dexbeyplus = abs(dexbey); + let permanent = ((cexdeyplus + dexceyplus) * bezplus + + (dexbeyplus + bexdeyplus) * cezplus + + (bexceyplus + cexbeyplus) * dezplus) + * alift + + ((dexaeyplus + aexdeyplus) * cezplus + + (aexceyplus + cexaeyplus) * dezplus + + (cexdeyplus + dexceyplus) * aezplus) + * blift + + ((aexbeyplus + bexaeyplus) * dezplus + + (bexdeyplus + dexbeyplus) * aezplus + + (dexaeyplus + aexdeyplus) * bezplus) + * clift + + ((bexceyplus + cexbeyplus) * aezplus + + (cexaeyplus + aexceyplus) * bezplus + + (aexbeyplus + bexaeyplus) * cezplus) + * dlift; + let errbound = ISPERRBOUND_A * permanent; + if det > errbound || -det > errbound { + return det; + } + + insphereadapt(pa, pb, pc, pd, pe, permanent) +} + +fn insphereadapt>( + pa: Coord3D, + pb: Coord3D, + pc: Coord3D, + pd: Coord3D, + pe: Coord3D, + permanent: f64, +) -> f64 { + let aex = pa.x - pe.x; + let bex = pb.x - pe.x; + let cex = pc.x - pe.x; + let dex = pd.x - pe.x; + let aey = pa.y - pe.y; + let bey = pb.y - pe.y; + let cey = pc.y - pe.y; + let dey = pd.y - pe.y; + let aez = pa.z - pe.z; + let bez = pb.z - pe.z; + let cez = pc.z - pe.z; + let dez = pd.z - pe.z; + + let (aexbey1, aexbey0) = two_product(aex, bey); + let (bexaey1, bexaey0) = two_product(bex, aey); + let (ab3, ab2, ab1, ab0) = two_two_diff(aexbey1, aexbey0, bexaey1, bexaey0); + let ab = [ab0, ab1, ab2, ab3]; + + let (bexcey1, bexcey0) = two_product(bex, cey); + let (cexbey1, cexbey0) = two_product(cex, bey); + let (bc3, bc2, bc1, bc0) = two_two_diff(bexcey1, bexcey0, cexbey1, cexbey0); + let bc = [bc0, bc1, bc2, bc3]; + + let (cexdey1, cexdey0) = two_product(cex, dey); + let (dexcey1, dexcey0) = two_product(dex, cey); + let (cd3, cd2, cd1, cd0) = two_two_diff(cexdey1, cexdey0, dexcey1, dexcey0); + let cd = [cd0, cd1, cd2, cd3]; + + let (dexaey1, dexaey0) = two_product(dex, aey); + let (aexdey1, aexdey0) = two_product(aex, dey); + let (da3, da2, da1, da0) = two_two_diff(dexaey1, dexaey0, aexdey1, aexdey0); + let da = [da0, da1, da2, da3]; + + let (aexcey1, aexcey0) = two_product(aex, cey); + let (cexaey1, cexaey0) = two_product(cex, aey); + let (ac3, ac2, ac1, ac0) = two_two_diff(aexcey1, aexcey0, cexaey1, cexaey0); + let ac = [ac0, ac1, ac2, ac3]; + + let (bexdey1, bexdey0) = two_product(bex, dey); + let (dexbey1, dexbey0) = two_product(dex, bey); + let (bd3, bd2, bd1, bd0) = two_two_diff(bexdey1, bexdey0, dexbey1, dexbey0); + let bd = [bd0, bd1, bd2, bd3]; + + + let mut temp8a = [0f64; 8]; + let mut temp8b = [0f64; 8]; + let mut temp8c = [0f64; 8]; + let mut temp16 = [0f64; 16]; + let mut temp24 = [0f64; 24]; + let mut temp48 = [0f64; 48]; + let mut xdet = [0f64; 96]; + let mut ydet = [0f64; 96]; + let mut zdet = [0f64; 96]; + let mut xydet = [0f64; 192]; + let mut adet = [0f64; 288]; + let mut bdet = [0f64; 288]; + let mut cdet = [0f64; 288]; + let mut ddet = [0f64; 288]; + let mut abdet = [0f64; 576]; + let mut cddet = [0f64; 576]; + let mut fin1 = [0f64; 1152]; + + let mut temp8alen = scale_expansion_zeroelim( &cd[..4], bez, &mut temp8a); + let mut temp8blen = scale_expansion_zeroelim(&bd[..4], -cez, &mut temp8b); + let mut temp8clen = scale_expansion_zeroelim( &bc[..4], dez, &mut temp8c); + let mut temp16len = fast_expansion_sum_zeroelim(&temp8a[..temp8alen], &temp8b[..temp8blen], &mut temp16); + let mut temp24len = fast_expansion_sum_zeroelim(&temp8c[..temp8clen], &temp16[..temp16len], &mut temp24); + let mut temp48len = scale_expansion_zeroelim(&temp24[..temp24len], aex, &mut temp48); + let mut xlen = scale_expansion_zeroelim(&temp48[..temp48len], -aex, &mut xdet); + let mut temp48len = scale_expansion_zeroelim(&temp24[..temp24len], aey, &mut temp48); + let mut ylen = scale_expansion_zeroelim(&temp48[..temp48len], -aey, &mut ydet); + let mut temp48len = scale_expansion_zeroelim(&temp24[..temp24len], aez, &mut temp48); + let mut zlen = scale_expansion_zeroelim(&temp48[..temp48len], -aez, &mut zdet); + let mut xylen = fast_expansion_sum_zeroelim(&xdet[..xlen], &ydet[..ylen], &mut xydet); + let mut alen = fast_expansion_sum_zeroelim(&xydet[..xylen], &zdet[..zlen], &mut adet); + + temp8alen = scale_expansion_zeroelim(&da[..4], cez, &mut temp8a); + temp8blen = scale_expansion_zeroelim(&ac[..4], dez, &mut temp8b); + temp8clen = scale_expansion_zeroelim(&cd[..4], aez, &mut temp8c); + temp16len = fast_expansion_sum_zeroelim(&temp8a[..temp8alen], &temp8b[..temp8blen], temp16); + temp24len = fast_expansion_sum_zeroelim(&temp8c[..temp8clen], &temp16[..temp16len], temp24); + temp48len = scale_expansion_zeroelim(&temp24[..temp24len], bex, &mut temp48); + xlen = scale_expansion_zeroelim(&temp48[..temp48len], bex, &mut xdet); + temp48len = scale_expansion_zeroelim(&temp24[..temp24len], bey, &mut temp48); + ylen = scale_expansion_zeroelim(&temp48[..temp48len], bey, &mut ydet); + temp48len = scale_expansion_zeroelim(&temp24[..temp24len], bez, &mut temp48); + zlen = scale_expansion_zeroelim(&temp48[..temp48len], bez, &mut zdet); + xylen = fast_expansion_sum_zeroelim(&xdet[..xlen], &ydet[..ylen], &mut xydet); + let blen = fast_expansion_sum_zeroelim(&xydet[..xylen], &zdet[..zlen], &mut bdet); + + temp8alen = scale_expansion_zeroelim(&ab[..4], dez, &mut temp8a); + temp8blen = scale_expansion_zeroelim(&bd[..4], aez, &mut temp8b); + temp8clen = scale_expansion_zeroelim(&da[..4], bez, &mut temp8c); + temp16len = fast_expansion_sum_zeroelim(&temp8a[..temp8alen], &temp8b[..temp8blen], temp16); + temp24len = fast_expansion_sum_zeroelim(&temp8c[..temp8clen], &temp16[..temp16len], temp24); + temp48len = scale_expansion_zeroelim(&temp24[..temp24len], cex, &mut temp48); + xlen = scale_expansion_zeroelim(&temp48[..temp48len], -cex, &mut xdet); + temp48len = scale_expansion_zeroelim(&temp24[..temp24len], cey, &mut temp48); + ylen = scale_expansion_zeroelim(&temp48[..temp48len], -cey, &mut ydet); + temp48len = scale_expansion_zeroelim(&temp24[..temp24len], cez, &mut temp48); + zlen = scale_expansion_zeroelim(&temp48[..temp48len], -cez, &mut zdet); + xylen = fast_expansion_sum_zeroelim(&xdet[..xlen], &ydet[..ylen], xydet); + let clen = fast_expansion_sum_zeroelim(&xydet[..xylen], &zdet[..zlen], cdet); + + temp8alen = scale_expansion_zeroelim(&bc[..4], aez, &mut temp8a); + temp8blen = scale_expansion_zeroelim(&ac[..4], -bez, &mut temp8b); + temp8clen = scale_expansion_zeroelim(&ab[..4], cez, &mut temp8c); + temp16len = fast_expansion_sum_zeroelim(&temp8a[..temp8alen], &temp8b[..temp8blen], &mut temp16); + temp24len = fast_expansion_sum_zeroelim(&temp8c[..temp8clen], &temp16[..temp16len], &mut temp24); + temp48len = scale_expansion_zeroelim(&temp24[..temp24len], dex, &mut temp48); + xlen = scale_expansion_zeroelim(&temp48[..temp48len], dex, &mut xdet); + temp48len = scale_expansion_zeroelim(&temp24[..temp24len], dey, &mut temp48); + ylen = scale_expansion_zeroelim(&temp48[..temp48len], dey, &mut ydet); + temp48len = scale_expansion_zeroelim(&temp24[..temp24len], dez, &mut temp48); + zlen = scale_expansion_zeroelim(&temp48[..temp48len], dez, &mut zdet); + xylen = fast_expansion_sum_zeroelim(&xdet[..xlen], &ydet[..ylen], &mut xydet); + let dlen = fast_expansion_sum_zeroelim(&xydet[..xylen], &zdet[..zlen], &mut ddet); + + let ablen = fast_expansion_sum_zeroelim(&adet[..alen], &bdet[..blen], &mut abdet); + let cdlen = fast_expansion_sum_zeroelim(&cdet[..clen], &ddet[..dlen], &mut cddet); + let finlength = fast_expansion_sum_zeroelim(&abdet[..ablen], &cddet[..cdlen], &mut fin1); + + let det = estimate(&fin1[..finlength]); + let mut errbound = ISPERRBOUND_B * permanent; + if det >= errbound || -det >= errbound { + return det; + } + + let aextail = two_diff_tail(pa.x, pe.x, aex); + let aeytail = two_diff_tail(pa.y, pe.y, aey); + let aeztail = two_diff_tail(pa.z, pe.z, aez); + let bextail = two_diff_tail(pb.x, pe.x, bex); + let beytail = two_diff_tail(pb.y, pe.y, bey); + let beztail = two_diff_tail(pb.z, pe.z, bez); + let cextail = two_diff_tail(pc.x, pe.x, cex); + let ceytail = two_diff_tail(pc.y, pe.y, cey); + let ceztail = two_diff_tail(pc.z, pe.z, cez); + let dextail = two_diff_tail(pd.x, pe.x, dex); + let deytail = two_diff_tail(pd.y, pe.y, dey); + let deztail = two_diff_tail(pd.z, pe.z, dez); + if (aextail == 0.0) + && (aeytail == 0.0) + && (aeztail == 0.0) + && (bextail == 0.0) + && (beytail == 0.0) + && (beztail == 0.0) + && (cextail == 0.0) + && (ceytail == 0.0) + && (ceztail == 0.0) + && (dextail == 0.0) + && (deytail == 0.0) + && (deztail == 0.0) + { + return det; + } + + errbound = ISPERRBOUND_C * permanent + RESULTERRBOUND * abs(det); + let abeps = (aex * beytail + bey * aextail) - (aey * bextail + bex * aeytail); + let bceps = (bex * ceytail + cey * bextail) - (bey * cextail + cex * beytail); + let cdeps = (cex * deytail + dey * cextail) - (cey * dextail + dex * ceytail); + let daeps = (dex * aeytail + aey * dextail) - (dey * aextail + aex * deytail); + let aceps = (aex * ceytail + cey * aextail) - (aey * cextail + cex * aeytail); + let bdeps = (bex * deytail + dey * bextail) - (bey * dextail + dex * beytail); + det += (((bex * bex + bey * bey + bez * bez) + * ((cez * daeps + dez * aceps + aez * cdeps) + + (ceztail * da3 + deztail * ac3 + aeztail * cd3)) + + (dex * dex + dey * dey + dez * dez) + * ((aez * bceps - bez * aceps + cez * abeps) + + (aeztail * bc3 - beztail * ac3 + ceztail * ab3))) + - ((aex * aex + aey * aey + aez * aez) + * ((bez * cdeps - cez * bdeps + dez * bceps) + + (beztail * cd3 - ceztail * bd3 + deztail * bc3)) + + (cex * cex + cey * cey + cez * cez) + * ((dez * abeps + aez * bdeps + bez * daeps) + + (deztail * ab3 + aeztail * bd3 + beztail * da3)))) + + 2.0 + * (((bex * bextail + bey * beytail + bez * beztail) + * (cez * da3 + dez * ac3 + aez * cd3) + + (dex * dextail + dey * deytail + dez * deztail) + * (aez * bc3 - bez * ac3 + cez * ab3)) + - ((aex * aextail + aey * aeytail + aez * aeztail) + * (bez * cd3 - cez * bd3 + dez * bc3) + + (cex * cextail + cey * ceytail + cez * ceztail) + * (dez * ab3 + aez * bd3 + bez * da3))); + if ((det >= errbound) || (-det >= errbound)) { + return det; + } + + insphereexact(pa, pb, pc, pd, pe); +} + +fn insphereexact>( + pa: Coord3D, + pb: Coord3D, + pc: Coord3D, + pd: Coord3D, + pe: Coord3D, +) -> f64 { + let (axby1, axby0) = two_product(pa.x, pb.y); + let (bxay1, bxay0) = two_product(pb.x, pa.y); + let (ab3, ab2, ab1, ab0) = two_two_diff(axby1, axby0, bxay1, bxay0); + let ab = [ab0, ab1, ab2, ab3]; + + let (bxcy1, bxcy0) = two_product(pb.x, pc.y,); + let (cxby1, cxby0) = two_product(pc.x, pb.y); + let (bc3, bc2, bc1, bc0) = two_two_diff(bxcy1, bxcy0, cxby1, cxby0); + let bc = [bc0, bc1, bc2, bc3]; + + let (cxdy1, cxdy0) = two_product(pc.x, pd.y); + let (dxcy1, dxcy0) = two_product(pd.x, pc.y); + let (cd3, cd2, cd1, cd0) = two_two_diff(cxdy1, cxdy0, dxcy1, dxcy0); + let cd = [cd0, cd1, cd2, cd3]; + + let (dxey1, dxey0) = two_product(pd.x, pe.y); + let (exdy1, exdy0) = two_product(pe.x, pd.y); + let (de3, de2, de1, de0) = two_two_diff(dxey1, dxey0, exdy1, exdy0); + let de = [de0, de1, de2, de3]; + + let (exay1, exay0) = two_product(pe.x, pa.y); + let (axey1, axey0) = two_product(pa.x, pe.y); + let (ea3, ea2, ea1, ea0) = two_two_diff(exay1, exay0, axey1, axey0); + let ea = [ea0, ea1, ea2, ea3]; + + let (axcy1, axcy0) = two_product(pa.x, pc.y); + let (cxay1, cxay0) = two_product(pc.x, pa.y); + let (ac3, ac2, ac1, ac0) = two_two_diff(axcy1, axcy0, cxay1, cxay0); + let ac = [ac0, ac1, ac2, ac3]; + + let (bxdy1, bxdy0) = two_product(pb.x, pd.y); + let (dxby1, dxby0) = two_product(pd.x, pb.y); + let (bd3, bd2, bd1, bd0) = two_two_diff(bxdy1, bxdy0, dxby1, dxby0); + let bd = [bd0, bd1, bd2, bd3]; + + let (cxey1, cxey0) = two_product(pc.x, pe.y); + let (excy1, excy0) = two_product(pe.x, pc.y); + let (ce3, ce2, ce1, ce0) = two_two_diff(cxey1, cxey0, excy1, excy0); + let ce = [ce0, ce1, ce2, ce3]; + + let (dxay1, dxay0) = two_product(pd.x, pa.y); + let (axdy1, axdy0) = two_product(pa.x, pd.y); + let (da3, da2, da1, da0) = two_two_diff(dxay1, dxay0, axdy1, axdy0); + let da = [da0, da1, da2, da3]; + + let (exby1, exby0) = two_product(pe.x, pb.y); + let (bxey1, bxey0) = two_product(pb.x, pe.y); + let (eb3, eb2, eb1, eb0) = two_two_diff(exby1, exby0, bxey1, bxey0); + let eb = [eb0, eb1, eb2, eb3]; + + let mut temp8alen = scale_expansion_zeroelim(4, bc, pa.z, temp8a); + let mut temp8blen = scale_expansion_zeroelim(4, ac, -pb.z, temp8b); + let mut temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp8blen, temp8b, + temp16); + temp8alen = scale_expansion_zeroelim(4, ab, pc.z, temp8a); + let abclen = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp16len, temp16, + abc); + + temp8alen = scale_expansion_zeroelim(4, cd, pb.z, temp8a); + temp8blen = scale_expansion_zeroelim(4, bd, -pc.z, temp8b); + temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp8blen, temp8b, + temp16); + temp8alen = scale_expansion_zeroelim(4, bc, pd.z, temp8a); + let bcdlen = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp16len, temp16, + bcd); + + temp8alen = scale_expansion_zeroelim(4, de, pc.z, temp8a); + temp8blen = scale_expansion_zeroelim(4, ce, -pd.z, temp8b); + temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp8blen, temp8b, + temp16); + temp8alen = scale_expansion_zeroelim(4, cd, pe.z, temp8a); + cdelen = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp16len, temp16, + cde); + + temp8alen = scale_expansion_zeroelim(4, ea, pd.z, temp8a); + temp8blen = scale_expansion_zeroelim(4, da, -pe.z, temp8b); + temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp8blen, temp8b, + temp16); + temp8alen = scale_expansion_zeroelim(4, de, pa.z, temp8a); + dealen = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp16len, temp16, + dea); + + temp8alen = scale_expansion_zeroelim(4, ab, pe.z, temp8a); + temp8blen = scale_expansion_zeroelim(4, eb, -pa.z, temp8b); + temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp8blen, temp8b, + temp16); + temp8alen = scale_expansion_zeroelim(4, ea, pb.z, temp8a); + eablen = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp16len, temp16, + eab); + + temp8alen = scale_expansion_zeroelim(4, bd, pa.z, temp8a); + temp8blen = scale_expansion_zeroelim(4, da, pb.z, temp8b); + temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp8blen, temp8b, + temp16); + temp8alen = scale_expansion_zeroelim(4, ab, pd.z, temp8a); + abdlen = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp16len, temp16, + abd); + + temp8alen = scale_expansion_zeroelim(4, ce, pb.z, temp8a); + temp8blen = scale_expansion_zeroelim(4, eb, pc.z, temp8b); + temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp8blen, temp8b, + temp16); + temp8alen = scale_expansion_zeroelim(4, bc, pe.z, temp8a); + bcelen = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp16len, temp16, + bce); + + temp8alen = scale_expansion_zeroelim(4, da, pc.z, temp8a); + temp8blen = scale_expansion_zeroelim(4, ac, pd.z, temp8b); + temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp8blen, temp8b, + temp16); + temp8alen = scale_expansion_zeroelim(4, cd, pa.z, temp8a); + cdalen = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp16len, temp16, + cda); + + temp8alen = scale_expansion_zeroelim(4, eb, pd.z, temp8a); + temp8blen = scale_expansion_zeroelim(4, bd, pe.z, temp8b); + temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp8blen, temp8b, + temp16); + temp8alen = scale_expansion_zeroelim(4, de, pb.z, temp8a); + deblen = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp16len, temp16, + deb); + + temp8alen = scale_expansion_zeroelim(4, ac, pe.z, temp8a); + temp8blen = scale_expansion_zeroelim(4, ce, pa.z, temp8b); + temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp8blen, temp8b, + temp16); + let temp8alen = scale_expansion_zeroelim(4, ea, pc.z, temp8a); + eaclen = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp16len, temp16, + eac); + + let temp48alen = fast_expansion_sum_zeroelim(cdelen, cde, bcelen, bce, temp48a); + let temp48blen = fast_expansion_sum_zeroelim(deblen, deb, bcdlen, bcd, temp48b); + for i in 0..temp48blen { + temp48b[i] = -temp48b[i]; + + } + + bcdelen = fast_expansion_sum_zeroelim(&temp48a[..temp48alen], + &temp48b[..temp48blen], &mut bcde); + xlen = scale_expansion_zeroelim(&bcde[..bcdelen], pa.x, &mut temp192); + xlen = scale_expansion_zeroelim(&temp192[..xlen], pa.x, &mut det384x); + ylen = scale_expansion_zeroelim(&bcde[..bcdelen], pa.y, &mut temp192); + ylen = scale_expansion_zeroelim(&temp192[..ylen], pa.y, &mut det384y); + zlen = scale_expansion_zeroelim(&bcde[..bcdelen], pa.z, &mut temp192); + zlen = scale_expansion_zeroelim(&temp192[..zlen], pa.z, &mut det384z); + xylen = fast_expansion_sum_zeroelim(&det384x[..xlen], &det384y[..ylen], &mut detxy); + alen = fast_expansion_sum_zeroelim(&detxy[..xylen], &det384z[..zlen], &mut adet); + + temp48alen = fast_expansion_sum_zeroelim(&dea[..dealen], &cda[..cdalen], &mut temp48a); + temp48blen = fast_expansion_sum_zeroelim(&eac[..eaclen], cde[..cdelen], &mut temp48b); + for (let i = 0; i < temp48blen; i++) { + temp48b[i] = -temp48b[i]; + } + cdealen = fast_expansion_sum_zeroelim(temp48alen, temp48a, + temp48blen, temp48b, &mut cdea); + xlen = scale_expansion_zeroelim(cdealen, cdea, pb.x, &mut temp192); + xlen = scale_expansion_zeroelim(xlen, temp192, pb.x, &mut det384x); + ylen = scale_expansion_zeroelim(cdealen, cdea, pb.y, &mut temp192); + ylen = scale_expansion_zeroelim(ylen, temp192, pb.y, &mut det384y); + zlen = scale_expansion_zeroelim(cdealen, cdea, pb.z, &mut temp192); + zlen = scale_expansion_zeroelim(zlen, temp192, pb.z, &mut det384z); + xylen = fast_expansion_sum_zeroelim(xlen, det384x, ylen, det384y, &mut detxy); + blen = fast_expansion_sum_zeroelim(xylen, detxy, zlen, det384z, &mut bdet); + + temp48alen = fast_expansion_sum_zeroelim(eablen, eab, deblen, deb, &mut temp48a); + temp48blen = fast_expansion_sum_zeroelim(abdlen, abd, dealen, dea, &mut temp48b); + for (i = 0; i < temp48blen; i++) { + temp48b[i] = -temp48b[i]; + } + deablen = fast_expansion_sum_zeroelim(temp48alen, temp48a, + temp48blen, temp48b, &mut deab); + xlen = scale_expansion_zeroelim(deablen, deab, pc.x, &mut temp192); + xlen = scale_expansion_zeroelim(xlen, temp192, pc.x, &mut det384x); + ylen = scale_expansion_zeroelim(deablen, deab, pc.y, &mut temp192); + ylen = scale_expansion_zeroelim(ylen, temp192, pc.y, &mut det384y); + zlen = scale_expansion_zeroelim(deablen, deab, pc.z, &mut temp192); + zlen = scale_expansion_zeroelim(zlen, temp192, pc.z, &mut det384z); + xylen = fast_expansion_sum_zeroelim(xlen, det384x, ylen, det384y, &mut detxy); + clen = fast_expansion_sum_zeroelim(xylen, detxy, zlen, det384z, &mut cdet); + + temp48alen = fast_expansion_sum_zeroelim(abclen, abc, eaclen, eac, &mut temp48a); + temp48blen = fast_expansion_sum_zeroelim(bcelen, bce, eablen, eab, &mut temp48b); + for (i = 0; i < temp48blen; i++) { + temp48b[i] = -temp48b[i]; + } + eabclen = fast_expansion_sum_zeroelim(temp48alen, temp48a, + temp48blen, temp48b, &mut eabc); + xlen = scale_expansion_zeroelim(eabclen, eabc, pd.x, &mut temp192); + xlen = scale_expansion_zeroelim(xlen, temp192, pd.x, &mut det384x); + ylen = scale_expansion_zeroelim(eabclen, eabc, pd.y, &mut temp192); + ylen = scale_expansion_zeroelim(ylen, temp192, pd.y, &mut det384y); + zlen = scale_expansion_zeroelim(eabclen, eabc, pd.z, &mut temp192); + zlen = scale_expansion_zeroelim(zlen, temp192, pd.z, &mut det384z); + xylen = fast_expansion_sum_zeroelim(xlen, det384x, ylen, det384y, &mut detxy); + dlen = fast_expansion_sum_zeroelim(xylen, detxy, zlen, det384z, &mut ddet); + + temp48alen = fast_expansion_sum_zeroelim(bcdlen, bcd, abdlen, abd, &mut temp48a); + temp48blen = fast_expansion_sum_zeroelim(cdalen, cda, abclen, abc, &mut temp48b); + for (i = 0; i < temp48blen; i++) { + temp48b[i] = -temp48b[i]; + } + abcdlen = fast_expansion_sum_zeroelim(temp48alen, temp48a, + temp48blen, temp48b, &mut abcd); + xlen = scale_expansion_zeroelim(abcdlen, abcd, pe.x, &mut temp192); + xlen = scale_expansion_zeroelim(xlen, temp192, pe.x, &mut det384x); + ylen = scale_expansion_zeroelim(abcdlen, abcd, pe.y, &mut temp192); + ylen = scale_expansion_zeroelim(ylen, temp192, pe.y, &mut det384y); + zlen = scale_expansion_zeroelim(abcdlen, abcd, pe.z, &mut temp192); + zlen = scale_expansion_zeroelim(zlen, temp192, pe.z, &mut det384z); + xylen = fast_expansion_sum_zeroelim(xlen, det384x, ylen, det384y, &mut detxy); + elen = fast_expansion_sum_zeroelim(xylen, detxy, zlen, det384z, &mut edet); + + ablen = fast_expansion_sum_zeroelim(alen, adet, blen, bdet, &mut abdet); + cdlen = fast_expansion_sum_zeroelim(clen, cdet, dlen, ddet, &mut cddet); + cdelen = fast_expansion_sum_zeroelim(cdlen, cddet, elen, edet, &mut cdedet); + deterlen = fast_expansion_sum_zeroelim(ablen, abdet, cdelen, cdedet, &mut deter); + + return deter[deterlen - 1]; +} + fn scale_expansion_zeroelim(e: &[f64], b: f64, h: &mut [f64]) -> usize { let (bhi, blo) = split(b); let (mut Q, hh) = two_product_presplit(e[0], b, bhi, blo);