#![cfg_attr(feature = "no_std", no_std)]

#![doc(html_logo_url = "https://raw.githubusercontent.com/georust/meta/master/logo/logo.png")]

// Copyright 2017 The Spade Developers.

// Copyright 2020 The GeoRust Project Developers.

//

// Licensed under the Apache License, Version 2.0 <LICENSE-APACHE or

// http://www.apache.org/licenses/LICENSE-2.0> or the MIT license

// <LICENSE-MIT or http://opensource.org/licenses/MIT>, at your

// option. This file may not be copied, modified, or distributed

// except according to those terms.

#![allow(non_snake_case)]

//! # Adaptive Precision Floating-Point Arithmetic and Fast Robust Predicates for Computational Geometry

//! This is a direct transcript of the source code and algorithms provided by

//! Jonathan Richard Shewchuk ([https://www.cs.cmu.edu/~quake/robust.html](https://www.cs.cmu.edu/~quake/robust.html))

//! See the paper and the source code for more information.

//!

//! The module offers adaptive and precise calculations for orientation queries

//! – "on which side of a line (2d) or plane (3d) does a point lie?" – and in-circle / in-sphere queries

//! – "is a given point contained in the circumference of a triangle?".  

//! The "adaptive" nature will increase performance only if a simpler calculation

//! cannot be guaranteed to be accurate enough, yielding higher performance on

//! average.

//!

//! The public API will accept both `f32` and `f64` input points for predicate checking, with input being converted to

//! `f64` values for internal use.

//! This has no effect on precision, as the [IEEE-754 standard](https://drive.google.com/file/d/0B3O3Ys97VjtxYXBCY08wanNoZ1U/view) (section 5.3)

//! guarantees that conversion from `f32` to `f64` must be exact.

//!

//! # Features

//! - `no_std`: Build without the Rust standard library

#[cfg(test)]

mod tests;

/// A two dimensional coordinate.

#[derive(Copy, Clone, Debug, PartialEq)]

pub struct Coord<T: Into<f64>> {

    pub x: T,

    pub y: T,

}

/// A three dimensional coordinate.

#[derive(Copy, Clone, Debug, PartialEq)]

pub struct Coord3D<T: Into<f64>> {

    pub x: T,

    pub y: T,

    pub z: T,

}

// These values are precomputed from the "exactinit" method of the c-source code. They should? be

// the same in all IEEE-754 environments, including rust f64

const SPLITTER: f64 = 134_217_729f64;

const EPSILON: f64 = 0.000_000_000_000_000_111_022_302_462_515_65;

const RESULTERRBOUND: f64 = (3.0 + 8.0 * EPSILON) * EPSILON;

const CCWERRBOUND_B: f64 = (2.0 + 12.0 * EPSILON) * EPSILON;

const CCWERRBOUND_C: f64 = (9.0 + 64.0 * EPSILON) * EPSILON * EPSILON;

const O3DERRBOUND_A: f64 = (7.0 + 56.0 * EPSILON) * EPSILON;

const O3DERRBOUND_B: f64 = (3.0 + 28.0 * EPSILON) * EPSILON;

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;

#[inline(always)]

fn abs(x: f64) -> f64 {

    x.abs()

}

/// 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`).  

/// Returns a negative value if they occur in clockwise order (`pc` lies to the **right** of the directed line `pa, pb`).  

/// Returns `0` if they are **collinear**.  

pub fn orient2d<T: Into<f64>>(pa: Coord<T>, pb: Coord<T>, pc: Coord<T>) -> f64 {

    let pa = Coord {

        x: pa.x.into(),

        y: pa.y.into(),

    };

    let pb = Coord {

        x: pb.x.into(),

        y: pb.y.into(),

    };

    let pc = Coord {

        x: pc.x.into(),

        y: pc.y.into(),

    };

    let detleft = (pa.x - pc.x) * (pb.y - pc.y);

    let detright = (pa.y - pc.y) * (pb.x - pc.x);

    let det = detleft - detright;

    // The errbound calculation was changed to require only one branch on the likely execution path.

    // This improves performance on modern processors as discussed by Ozaki et al in

    // https://doi.org/10.1007/s10543-015-0574-9

    // The underflow guard "+ u_N" was omitted because orient2dadapt(...) would not guarantee

    // correct results in cases of underflow, the derivation of THETA is given in the reference.

    let detsum = abs(detleft + detright);

    const THETA: f64 = 3.3306690621773722e-16;

    let errbound = THETA * detsum;

    if det >= errbound || -det >= errbound {

        det

    } else {

        orient2dadapt(pa, pb, pc, detsum)

    }

}

Unexecuted instantiation: robust::orient2d::<f64>

Unexecuted instantiation: robust::orient2d::<_>

fn orient2dadapt(pa: Coord<f64>, pb: Coord<f64>, pc: Coord<f64>, detsum: f64) -> f64 {

    let acx = pa.x - pc.x;

    let bcx = pb.x - pc.x;

    let acy = pa.y - pc.y;

    let bcy = pb.y - pc.y;

    let (detleft, detlefttail) = two_product(acx, bcy);

    let (detright, detrighttail) = two_product(acy, bcx);

    let (B3, B2, B1, B0) = two_two_diff(detleft, detlefttail, detright, detrighttail);

    let B = [B0, B1, B2, B3];

    let mut det = estimate(&B);

    let errbound = CCWERRBOUND_B * detsum;

    if det >= errbound || (-det >= errbound) {

        return det;

    }

    let acxtail = two_diff_tail(pa.x, pc.x, acx);

    let bcxtail = two_diff_tail(pb.x, pc.x, bcx);

    let acytail = two_diff_tail(pa.y, pc.y, acy);

    let bcytail = two_diff_tail(pb.y, pc.y, bcy);

    if acxtail == 0.0 && acytail == 0.0 && bcxtail == 0.0 && bcytail == 0.0 {

        return det;

    }

    let errbound = CCWERRBOUND_C * detsum + RESULTERRBOUND * abs(det);

    det += (acx * bcytail + bcy * acxtail) - (acy * bcxtail + bcx * acytail);

    if det >= errbound || -det >= errbound {

        return det;

    }

    let (s1, s0) = two_product(acxtail, bcy);

    let (t1, t0) = two_product(acytail, bcx);

    let (u3, u2, u1, u0) = two_two_diff(s1, s0, t1, t0);

    let U = [u0, u1, u2, u3];

    let mut C1 = [0.0f64; 8];

    let c1length = fast_expansion_sum_zeroelim(&B, &U, &mut C1);

    let (s1, s0) = two_product(acx, bcytail);

    let (t1, t0) = two_product(acy, bcxtail);

    let (u3, u2, u1, u0) = two_two_diff(s1, s0, t1, t0);

    let U = [u0, u1, u2, u3];

    let mut C2 = [0.0f64; 12];

    let c2length = fast_expansion_sum_zeroelim(&C1[..c1length], &U, &mut C2);

    let (s1, s0) = two_product(acxtail, bcytail);

    let (t1, t0) = two_product(acytail, bcxtail);

    let (u3, u2, u1, u0) = two_two_diff(s1, s0, t1, t0);

    let U = [u0, u1, u2, u3];

    let mut D = [0.0f64; 16];

    let dlength = fast_expansion_sum_zeroelim(&C2[..c2length], &U, &mut D);

    D[dlength - 1]

}

/// Returns a positive value if the point `pd` lies below the plane passing through `pa`, `pb`, and `pc`

/// ("below" is defined so that `pa`, `pb`, and `pc` appear in counterclockwise order when viewed from above the plane).  

/// Returns a negative value if `pd` lies above the plane.  

/// Returns `0` if they are **coplanar**.  

pub fn orient3d<T: Into<f64>>(

    pa: Coord3D<T>,

    pb: Coord3D<T>,

    pc: Coord3D<T>,

    pd: Coord3D<T>,

) -> 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 adx = pa.x - pd.x;

    let bdx = pb.x - pd.x;

    let cdx = pc.x - pd.x;

    let ady = pa.y - pd.y;

    let bdy = pb.y - pd.y;

    let cdy = pc.y - pd.y;

    let adz = pa.z - pd.z;

    let bdz = pb.z - pd.z;

    let cdz = pc.z - pd.z;

    let bdxcdy = bdx * cdy;

    let cdxbdy = cdx * bdy;

    let cdxady = cdx * ady;

    let adxcdy = adx * cdy;

    let adxbdy = adx * bdy;

    let bdxady = bdx * ady;

    let det = adz * (bdxcdy - cdxbdy) + bdz * (cdxady - adxcdy) + cdz * (adxbdy - bdxady);

    let permanent = (abs(bdxcdy) + abs(cdxbdy)) * abs(adz)

        + (abs(cdxady) + abs(adxcdy)) * abs(bdz)

        + (abs(adxbdy) + abs(bdxady)) * abs(cdz);

    let errbound = O3DERRBOUND_A * permanent;

    if det > errbound || -det > errbound {

        return det;

    }

    orient3dadapt(pa, pb, pc, pd, permanent)

}

fn orient3dadapt(

    pa: Coord3D<f64>,

    pb: Coord3D<f64>,

    pc: Coord3D<f64>,

    pd: Coord3D<f64>,

    permanent: f64,

) -> f64 {

    let adx = pa.x - pd.x;

    let bdx = pb.x - pd.x;

    let cdx = pc.x - pd.x;

    let ady = pa.y - pd.y;

    let bdy = pb.y - pd.y;

    let cdy = pc.y - pd.y;

    let adz = pa.z - pd.z;

    let bdz = pb.z - pd.z;

    let cdz = pc.z - pd.z;

    let (bdxcdy1, bdxcdy0) = two_product(bdx, cdy);

    let (cdxbdy1, cdxbdy0) = two_product(cdx, bdy);

    let (bc3, bc2, bc1, bc0) = two_two_diff(bdxcdy1, bdxcdy0, cdxbdy1, cdxbdy0);

    let bc = [bc0, bc1, bc2, bc3];

    let mut adet = [0f64; 8];

    let alen = scale_expansion_zeroelim(&bc[..4], adz, &mut adet);

    let (cdxady1, cdxady0) = two_product(cdx, ady);

    let (adxcdy1, adxcdy0) = two_product(adx, cdy);

    let (ca3, ca2, ca1, ca0) = two_two_diff(cdxady1, cdxady0, adxcdy1, adxcdy0);

    let ca = [ca0, ca1, ca2, ca3];

    let mut bdet = [0f64; 8];

    let blen = scale_expansion_zeroelim(&ca[..4], bdz, &mut bdet);

    let (adxbdy1, adxbdy0) = two_product(adx, bdy);

    let (bdxady1, bdxady0) = two_product(bdx, ady);

    let (ab3, ab2, ab1, ab0) = two_two_diff(adxbdy1, adxbdy0, bdxady1, bdxady0);

    let ab = [ab0, ab1, ab2, ab3];

    let mut cdet = [0f64; 8];

    let clen = scale_expansion_zeroelim(&ab[..4], cdz, &mut cdet);

    let mut abdet = [0f64; 16];

    let ablen = fast_expansion_sum_zeroelim(&adet[..alen], &bdet[..blen], &mut abdet);

    let mut fin1 = [0f64; 192];

    let finlength = fast_expansion_sum_zeroelim(&abdet[..ablen], &cdet[..clen], &mut fin1);

    let mut det = estimate(&fin1[..finlength]);

    let mut errbound = O3DERRBOUND_B * permanent;

    if det >= errbound || -det >= errbound {

        return det;

    }

    let adxtail = two_diff_tail(pa.x, pd.x, adx);

    let bdxtail = two_diff_tail(pb.x, pd.x, bdx);

    let cdxtail = two_diff_tail(pc.x, pd.x, cdx);

    let adytail = two_diff_tail(pa.y, pd.y, ady);

    let bdytail = two_diff_tail(pb.y, pd.y, bdy);

    let cdytail = two_diff_tail(pc.y, pd.y, cdy);

    let adztail = two_diff_tail(pa.z, pd.z, adz);

    let bdztail = two_diff_tail(pb.z, pd.z, bdz);

    let cdztail = two_diff_tail(pc.z, pd.z, cdz);

    if (adxtail == 0.0)

        && (bdxtail == 0.0)

        && (cdxtail == 0.0)

        && (adytail == 0.0)

        && (bdytail == 0.0)

        && (cdytail == 0.0)

        && (adztail == 0.0)

        && (bdztail == 0.0)

        && (cdztail == 0.0)

    {

        return det;

    }

    errbound = O3DERRBOUND_C * permanent + RESULTERRBOUND * abs(det);

    det += (adz * ((bdx * cdytail + cdy * bdxtail) - (bdy * cdxtail + cdx * bdytail))

        + adztail * (bdx * cdy - bdy * cdx))

        + (bdz * ((cdx * adytail + ady * cdxtail) - (cdy * adxtail + adx * cdytail))

            + bdztail * (cdx * ady - cdy * adx))

        + (cdz * ((adx * bdytail + bdy * adxtail) - (ady * bdxtail + bdx * adytail))

            + cdztail * (adx * bdy - ady * bdx));

    if det >= errbound || -det >= errbound {

        return det;

    }

    let mut finnow = fin1;

    let mut finother = [0f64; 192];

    let mut at_b = [0f64; 4];

    let mut at_c = [0f64; 4];

    let mut bt_c = [0f64; 4];

    let mut bt_a = [0f64; 4];

    let mut ct_a = [0f64; 4];

    let mut ct_b = [0f64; 4];

    let at_blen: usize;

    let at_clen: usize;

    let bt_clen: usize;

    let bt_alen: usize;

    let ct_alen: usize;

    let ct_blen: usize;

    if adxtail == 0.0 {

        if adytail == 0.0 {

            at_b[0] = 0.0;

            at_blen = 1;

            at_c[0] = 0.0;

            at_clen = 1;

        } else {

            let negate = -adytail;

            (at_b[1], at_b[0]) = two_product(negate, bdx);

            at_blen = 2;

            (at_c[1], at_c[0]) = two_product(adytail, cdx);

            at_clen = 2;

        }

    } else if adytail == 0.0 {

        (at_b[1], at_b[0]) = two_product(adxtail, bdy);

        at_blen = 2;

        let negate = -adxtail;

        (at_c[1], at_c[0]) = two_product(negate, cdy);

        at_clen = 2;

    } else {

        let (adxt_bdy1, adxt_bdy0) = two_product(adxtail, bdy);

        let (adyt_bdx1, adyt_bdx0) = two_product(adytail, bdx);

        (at_b[3], at_b[2], at_b[1], at_b[0]) =

            two_two_diff(adxt_bdy1, adxt_bdy0, adyt_bdx1, adyt_bdx0);

        at_blen = 4;

        let (adyt_cdx1, adyt_cdx0) = two_product(adytail, cdx);

        let (adxt_cdy1, adxt_cdy0) = two_product(adxtail, cdy);

        (at_c[3], at_c[2], at_c[1], at_c[0]) =

            two_two_diff(adyt_cdx1, adyt_cdx0, adxt_cdy1, adxt_cdy0);

        at_clen = 4;

    }

    if bdxtail == 0.0 {

        if bdytail == 0.0 {

            bt_c[0] = 0.0;

            bt_clen = 1;

            bt_a[0] = 0.0;

            bt_alen = 1;

        } else {

            let negate = -bdytail;

            (bt_c[1], bt_c[0]) = two_product(negate, cdx);

            bt_clen = 2;

            (bt_a[1], bt_a[0]) = two_product(bdytail, adx);

            bt_alen = 2;

        }

    } else if bdytail == 0.0 {

        (bt_c[1], bt_c[0]) = two_product(bdxtail, cdy);

        bt_clen = 2;

        let negate = -bdxtail;

        (bt_a[1], bt_a[0]) = two_product(negate, ady);

        bt_alen = 2;

    } else {

        let (bdxt_cdy1, bdxt_cdy0) = two_product(bdxtail, cdy);

        let (bdyt_cdx1, bdyt_cdx0) = two_product(bdytail, cdx);

        (bt_c[3], bt_c[2], bt_c[1], bt_c[0]) =

            two_two_diff(bdxt_cdy1, bdxt_cdy0, bdyt_cdx1, bdyt_cdx0);

        bt_clen = 4;

        let (bdyt_adx1, bdyt_adx0) = two_product(bdytail, adx);

        let (bdxt_ady1, bdxt_ady0) = two_product(bdxtail, ady);

        (bt_a[3], bt_a[2], bt_a[1], bt_a[0]) =

            two_two_diff(bdyt_adx1, bdyt_adx0, bdxt_ady1, bdxt_ady0);

        bt_alen = 4;

    }

    if cdxtail == 0.0 {

        if cdytail == 0.0 {

            ct_a[0] = 0.0;

            ct_alen = 1;

            ct_b[0] = 0.0;

            ct_blen = 1;

        } else {

            let negate = -cdytail;

            (ct_a[1], ct_a[0]) = two_product(negate, adx);

            ct_alen = 2;

            (ct_b[1], ct_b[0]) = two_product(cdytail, bdx);

            ct_blen = 2;

        }

    } else if cdytail == 0.0 {

        (ct_a[1], ct_a[0]) = two_product(cdxtail, ady);

        ct_alen = 2;

        let negate = -cdxtail;

        (ct_b[1], ct_b[0]) = two_product(negate, bdy);

        ct_blen = 2;

    } else {

        let (cdxt_ady1, cdxt_ady0) = two_product(cdxtail, ady);

        let (cdyt_adx1, cdyt_adx0) = two_product(cdytail, adx);

        (ct_a[3], ct_a[2], ct_a[1], ct_a[0]) =

            two_two_diff(cdxt_ady1, cdxt_ady0, cdyt_adx1, cdyt_adx0);

        ct_alen = 4;

        let (cdyt_bdx1, cdyt_bdx0) = two_product(cdytail, bdx);

        let (cdxt_bdy1, cdxt_bdy0) = two_product(cdxtail, bdy);

        (ct_b[3], ct_b[2], ct_b[1], ct_b[0]) =

            two_two_diff(cdyt_bdx1, cdyt_bdx0, cdxt_bdy1, cdxt_bdy0);

        ct_blen = 4;

    }

    let mut bct = [0f64; 8];

    let mut cat = [0f64; 8];

    let mut abt = [0f64; 8];

    let mut u = [0f64; 4];

    let mut v = [0f64; 12];

    let mut w = [0f64; 16];

    let mut vlength: usize;

    let mut wlength: usize;

    let bctlen = fast_expansion_sum_zeroelim(&bt_c[..bt_clen], &ct_b[..ct_blen], &mut bct);

    wlength = scale_expansion_zeroelim(&bct[..bctlen], adz, &mut w);

    let mut finlength =

        fast_expansion_sum_zeroelim(&finnow[..finlength], &w[..wlength], &mut finother);

    ::core::mem::swap(&mut finnow, &mut finother);

    let catlen = fast_expansion_sum_zeroelim(&ct_a[..ct_alen], &at_c[..at_clen], &mut cat);

    wlength = scale_expansion_zeroelim(&cat[..catlen], bdz, &mut w);

    finlength = fast_expansion_sum_zeroelim(&finnow[..finlength], &w[..wlength], &mut finother);

    ::core::mem::swap(&mut finnow, &mut finother);

    let abtlen = fast_expansion_sum_zeroelim(&at_b[..at_blen], &bt_a[..bt_alen], &mut abt);

    wlength = scale_expansion_zeroelim(&abt[..abtlen], cdz, &mut w);

    finlength = fast_expansion_sum_zeroelim(&finnow[..finlength], &w[..wlength], &mut finother);

    ::core::mem::swap(&mut finnow, &mut finother);

    if adztail != 0.0 {

        vlength = scale_expansion_zeroelim(&bc[..4], adztail, &mut v);

        finlength = fast_expansion_sum_zeroelim(&finnow[..finlength], &v[..vlength], &mut finother);

        ::core::mem::swap(&mut finnow, &mut finother);

    }

    if bdztail != 0.0 {

        vlength = scale_expansion_zeroelim(&ca[..4], bdztail, &mut v);

        finlength = fast_expansion_sum_zeroelim(&finnow[..finlength], &v[..vlength], &mut finother);

        ::core::mem::swap(&mut finnow, &mut finother);

    }

    if cdztail != 0.0 {

        vlength = scale_expansion_zeroelim(&ab[..4], cdztail, &mut v);

        finlength = fast_expansion_sum_zeroelim(&finnow[..finlength], &v[..vlength], &mut finother);

        ::core::mem::swap(&mut finnow, &mut finother);

    }

    if adxtail != 0.0 {

        if bdytail != 0.0 {

            let (adxt_bdyt1, adxt_bdyt0) = two_product(adxtail, bdytail);

            (u[3], u[2], u[1], u[0]) = two_one_product(adxt_bdyt1, adxt_bdyt0, cdz);

            finlength = fast_expansion_sum_zeroelim(&finnow[..finlength], &u[..4], &mut finother);

            ::core::mem::swap(&mut finnow, &mut finother);

            if cdztail != 0.0 {

                (u[3], u[2], u[1], u[0]) = two_one_product(adxt_bdyt1, adxt_bdyt0, cdztail);

                finlength =

                    fast_expansion_sum_zeroelim(&finnow[..finlength], &u[..4], &mut finother);

                ::core::mem::swap(&mut finnow, &mut finother);

            }

        }

        if cdytail != 0.0 {

            let negate = -adxtail;

            let (adxt_cdyt1, adxt_cdyt0) = two_product(negate, cdytail);

            (u[3], u[2], u[1], u[0]) = two_one_product(adxt_cdyt1, adxt_cdyt0, bdz);

            finlength = fast_expansion_sum_zeroelim(&finnow[..finlength], &u[..4], &mut finother);

            ::core::mem::swap(&mut finnow, &mut finother);

            if bdztail != 0.0 {

                (u[3], u[2], u[1], u[0]) = two_one_product(adxt_cdyt1, adxt_cdyt0, bdztail);

                finlength =

                    fast_expansion_sum_zeroelim(&finnow[..finlength], &u[..4], &mut finother);

                ::core::mem::swap(&mut finnow, &mut finother);

            }

        }

    }

    if bdxtail != 0.0 {

        if cdytail != 0.0 {

            let (bdxt_cdyt1, bdxt_cdyt0) = two_product(bdxtail, cdytail);

            (u[3], u[2], u[1], u[0]) = two_one_product(bdxt_cdyt1, bdxt_cdyt0, adz);

            finlength = fast_expansion_sum_zeroelim(&finnow[..finlength], &u[..4], &mut finother);

            ::core::mem::swap(&mut finnow, &mut finother);

            if adztail != 0.0 {

                (u[3], u[2], u[1], u[0]) = two_one_product(bdxt_cdyt1, bdxt_cdyt0, adztail);

                finlength =

                    fast_expansion_sum_zeroelim(&finnow[..finlength], &u[..4], &mut finother);

                ::core::mem::swap(&mut finnow, &mut finother);

            }

        }

        if adytail != 0.0 {

            let negate = -bdxtail;

            let (bdxt_adyt1, bdxt_adyt0) = two_product(negate, adytail);

            (u[3], u[2], u[1], u[0]) = two_one_product(bdxt_adyt1, bdxt_adyt0, cdz);

            finlength = fast_expansion_sum_zeroelim(&finnow[..finlength], &u[..4], &mut finother);

            ::core::mem::swap(&mut finnow, &mut finother);

            if cdztail != 0.0 {

                (u[3], u[2], u[1], u[0]) = two_one_product(bdxt_adyt1, bdxt_adyt0, cdztail);

                finlength =

                    fast_expansion_sum_zeroelim(&finnow[..finlength], &u[..4], &mut finother);

                ::core::mem::swap(&mut finnow, &mut finother);

            }

        }

    }

    if cdxtail != 0.0 {

        if adytail != 0.0 {

            let (cdxt_adyt1, cdxt_adyt0) = two_product(cdxtail, adytail);

            (u[3], u[2], u[1], u[0]) = two_one_product(cdxt_adyt1, cdxt_adyt0, bdz);

            finlength = fast_expansion_sum_zeroelim(&finnow[..finlength], &u[..4], &mut finother);

            ::core::mem::swap(&mut finnow, &mut finother);

            if bdztail != 0.0 {

                (u[3], u[2], u[1], u[0]) = two_one_product(cdxt_adyt1, cdxt_adyt0, bdztail);

                finlength =

                    fast_expansion_sum_zeroelim(&finnow[..finlength], &u[..4], &mut finother);

                ::core::mem::swap(&mut finnow, &mut finother);

            }

        }

        if bdytail != 0.0 {

            let negate = -cdxtail;

            let (cdxt_bdyt1, cdxt_bdyt0) = two_product(negate, bdytail);

            (u[3], u[2], u[1], u[0]) = two_one_product(cdxt_bdyt1, cdxt_bdyt0, adz);

            finlength = fast_expansion_sum_zeroelim(&finnow[..finlength], &u[..4], &mut finother);

            ::core::mem::swap(&mut finnow, &mut finother);

            if adztail != 0.0 {

                (u[3], u[2], u[1], u[0]) = two_one_product(cdxt_bdyt1, cdxt_bdyt0, adztail);

                finlength =

                    fast_expansion_sum_zeroelim(&finnow[..finlength], &u[..4], &mut finother);

                ::core::mem::swap(&mut finnow, &mut finother);

            }

        }

    }

    if adztail != 0.0 {

        wlength = scale_expansion_zeroelim(&bct[..bctlen], adztail, &mut w);

        finlength = fast_expansion_sum_zeroelim(&finnow[..finlength], &w[..wlength], &mut finother);

        ::core::mem::swap(&mut finnow, &mut finother);

    }

    if bdztail != 0.0 {

        wlength = scale_expansion_zeroelim(&cat[..catlen], bdztail, &mut w);

        finlength = fast_expansion_sum_zeroelim(&finnow[..finlength], &w[..wlength], &mut finother);

        ::core::mem::swap(&mut finnow, &mut finother);

    }

    if cdztail != 0.0 {

        wlength = scale_expansion_zeroelim(&abt[..abtlen], cdztail, &mut w);

        finlength = fast_expansion_sum_zeroelim(&finnow[..finlength], &w[..wlength], &mut finother);

        ::core::mem::swap(&mut finnow, &mut finother);

    }

    finnow[finlength - 1]

}

/// Returns a positive value if the coordinate `pd` lies **inside** the circle passing through `pa`, `pb`, and `pc`.  

/// Returns a negative value if it lies **outside** the circle.  

/// Returns `0` if the four points are **cocircular**.  

/// **Note**: The points `pa`, `pb`, and `pc` must be in **counterclockwise order**, or the sign of the result will be reversed.

pub fn incircle<T: Into<f64>>(pa: Coord<T>, pb: Coord<T>, pc: Coord<T>, pd: Coord<T>) -> f64 {

    let pa = Coord {

        x: pa.x.into(),

        y: pa.y.into(),

    };

    let pb = Coord {

        x: pb.x.into(),

        y: pb.y.into(),

    };

    let pc = Coord {

        x: pc.x.into(),

        y: pc.y.into(),

    };

    let pd = Coord {

        x: pd.x.into(),

        y: pd.y.into(),

    };

    let adx = pa.x - pd.x;

    let bdx = pb.x - pd.x;

    let cdx = pc.x - pd.x;

    let ady = pa.y - pd.y;

    let bdy = pb.y - pd.y;

    let cdy = pc.y - pd.y;

    let bdxcdy = bdx * cdy;

    let cdxbdy = cdx * bdy;

    let alift = adx * adx + ady * ady;

    let cdxady = cdx * ady;

    let adxcdy = adx * cdy;

    let blift = bdx * bdx + bdy * bdy;

    let adxbdy = adx * bdy;

    let bdxady = bdx * ady;

    let clift = cdx * cdx + cdy * cdy;

    let det = alift * (bdxcdy - cdxbdy) + blift * (cdxady - adxcdy) + clift * (adxbdy - bdxady);

    let permanent = (abs(bdxcdy) + abs(cdxbdy)) * alift

        + (abs(cdxady) + abs(adxcdy)) * blift

        + (abs(adxbdy) + abs(bdxady)) * clift;

    let errbound = ICCERRBOUND_A * permanent;

    if det > errbound || -det > errbound {

        return det;

    }

    incircleadapt(pa, pb, pc, pd, permanent)

}

fn incircleadapt(

    pa: Coord<f64>,

    pb: Coord<f64>,

    pc: Coord<f64>,

    pd: Coord<f64>,

    permanent: f64,

) -> f64 {

    let mut temp8 = [0f64; 8];

    let mut temp16a = [0f64; 16];

    let mut temp16b = [0f64; 16];

    let mut temp16c = [0f64; 16];

    let mut temp32a = [0f64; 32];

    let mut temp32b = [0f64; 32];

    let mut temp48 = [0f64; 48];

    let mut temp64 = [0f64; 64];

    let adx = pa.x - pd.x;

    let bdx = pb.x - pd.x;

    let cdx = pc.x - pd.x;

    let ady = pa.y - pd.y;

    let bdy = pb.y - pd.y;

    let cdy = pc.y - pd.y;

    let (bdxcdy1, bdxcdy0) = two_product(bdx, cdy);

    let (cdxbdy1, cdxbdy0) = two_product(cdx, bdy);

    let (bc3, bc2, bc1, bc0) = two_two_diff(bdxcdy1, bdxcdy0, cdxbdy1, cdxbdy0);

    let bc = [bc0, bc1, bc2, bc3];

    let mut axbc = [0f64; 8];

    let axbclen = scale_expansion_zeroelim(&bc, adx, &mut axbc);

    let mut axxbc = [0f64; 16];

    let axxbclen = scale_expansion_zeroelim(&axbc[..axbclen], adx, &mut axxbc);

    let mut aybc = [0f64; 8];

    let aybclen = scale_expansion_zeroelim(&bc, ady, &mut aybc);

    let mut ayybc = [0f64; 16];

    let ayybclen = scale_expansion_zeroelim(&aybc[..aybclen], ady, &mut ayybc);

    let mut adet = [0f64; 32];

    let alen = fast_expansion_sum_zeroelim(&axxbc[0..axxbclen], &ayybc[0..ayybclen], &mut adet);

    let (cdxady1, cdxady0) = two_product(cdx, ady);

    let (adxcdy1, adxcdy0) = two_product(adx, cdy);

    let (c3, c2, c1, c0) = two_two_diff(cdxady1, cdxady0, adxcdy1, adxcdy0);

    let ca = [c0, c1, c2, c3];

    let mut bxca = [0f64; 8];

    let bxcalen = scale_expansion_zeroelim(&ca, bdx, &mut bxca);

    let mut bxxca = [0f64; 16];

    let bxxcalen = scale_expansion_zeroelim(&bxca[..bxcalen], bdx, &mut bxxca);

    let mut byca = [0f64; 8];

    let bycalen = scale_expansion_zeroelim(&ca, bdy, &mut byca);

    let mut byyca = [0f64; 16];

    let byycalen = scale_expansion_zeroelim(&byca[..bycalen], bdy, &mut byyca);

    let mut bdet = [0f64; 32];

    let blen = fast_expansion_sum_zeroelim(&bxxca[..bxxcalen], &byyca[0..byycalen], &mut bdet);

    let (adxbdy1, adxbdy0) = two_product(adx, bdy);

    let (bdxady1, bdxady0) = two_product(bdx, ady);

    let (ab3, ab2, ab1, ab0) = two_two_diff(adxbdy1, adxbdy0, bdxady1, bdxady0);

    let ab = [ab0, ab1, ab2, ab3];

    let mut cxab = [0f64; 8];

    let cxablen = scale_expansion_zeroelim(&ab, cdx, &mut cxab);

    let mut cxxab = [0f64; 16];

    let cxxablen = scale_expansion_zeroelim(&cxab[..cxablen], cdx, &mut cxxab);

    let mut cyab = [0f64; 8];

    let cyablen = scale_expansion_zeroelim(&ab, cdy, &mut cyab);

    let mut cyyab = [0f64; 16];

    let cyyablen = scale_expansion_zeroelim(&cyab[..cyablen], cdy, &mut cyyab);

    let mut cdet = [0f64; 32];

    let clen = fast_expansion_sum_zeroelim(&cxxab[..cxxablen], &cyyab[..cyyablen], &mut cdet);

    let mut abdet = [0f64; 64];

    let ablen = fast_expansion_sum_zeroelim(&adet[..alen], &bdet[..blen], &mut abdet);

    let mut fin1 = [0f64; 1152];

    let mut finlength = fast_expansion_sum_zeroelim(&abdet[..ablen], &cdet[..clen], &mut fin1);

    let mut det = estimate(&fin1[..finlength]);

    let errbound = ICCERRBOUND_B * permanent;

    if det >= errbound || -det >= errbound {

        return det;

    }

    let adxtail = two_diff_tail(pa.x, pd.x, adx);

    let adytail = two_diff_tail(pa.y, pd.y, ady);

    let bdxtail = two_diff_tail(pb.x, pd.x, bdx);

    let bdytail = two_diff_tail(pb.y, pd.y, bdy);

    let cdxtail = two_diff_tail(pc.x, pd.x, cdx);

    let cdytail = two_diff_tail(pc.y, pd.y, cdy);

    if adxtail == 0.0

        && bdxtail == 0.0

        && cdxtail == 0.0

        && adytail == 0.0

        && bdytail == 0.0

        && cdytail == 0.0

    {

        return det;

    }

    let errbound = ICCERRBOUND_C * permanent + RESULTERRBOUND * abs(det);

    det += ((adx * adx + ady * ady)

        * ((bdx * cdytail + cdy * bdxtail) - (bdy * cdxtail + cdx * bdytail))

        + 2.0 * (adx * adxtail + ady * adytail) * (bdx * cdy - bdy * cdx))

        + ((bdx * bdx + bdy * bdy)

            * ((cdx * adytail + ady * cdxtail) - (cdy * adxtail + adx * cdytail))

            + 2.0 * (bdx * bdxtail + bdy * bdytail) * (cdx * ady - cdy * adx))

        + ((cdx * cdx + cdy * cdy)

            * ((adx * bdytail + bdy * adxtail) - (ady * bdxtail + bdx * adytail))

            + 2.0 * (cdx * cdxtail + cdy * cdytail) * (adx * bdy - ady * bdx));

    if det >= errbound || -det >= errbound {

        return det;

    }

    let mut fin2 = [0f64; 1152];

    let aa = if bdxtail != 0.0 || bdytail != 0.0 || cdxtail != 0.0 || cdytail != 0.0 {

        let (adxadx1, adxadx0) = square(adx);

        let (adyady1, adyady0) = square(ady);

        let (aa3, aa2, aa1, aa0) = two_two_sum(adxadx1, adxadx0, adyady1, adyady0);

        [aa0, aa1, aa2, aa3]

    } else {

        [0f64; 4]

    };

    let bb = if cdxtail != 0.0 || cdytail != 0.0 || adxtail != 0.0 || adytail != 0.0 {

        let (bdxbdx1, bdxbdx0) = square(bdx);

        let (bdybdy1, bdybdy0) = square(bdy);

        let (bb3, bb2, bb1, bb0) = two_two_sum(bdxbdx1, bdxbdx0, bdybdy1, bdybdy0);

        [bb0, bb1, bb2, bb3]

    } else {

        [0f64; 4]

    };

    let cc = if adxtail != 0.0 || adytail != 0.0 || bdxtail != 0.0 || bdytail != 0.0 {

        let (cdxcdx1, cdxcdx0) = square(cdx);

        let (cdycdy1, cdycdy0) = square(cdy);

        let (cc3, cc2, cc1, cc0) = two_two_sum(cdxcdx1, cdxcdx0, cdycdy1, cdycdy0);

        [cc0, cc1, cc2, cc3]

    } else {

        [0f64; 4]

    };