///////// Adaptive precision floating point arithmetic //////////

double sum(double a, double b, double & roundoff)
{
    double res = a + b;
    double bv = res - a;
    double av = res - bv;
    double br = b - bv;
    double ar = a - av;
    roundoff = ar + br;
    return res;
}

template <size_t s>
void split(double a, double & hi, double & lo)
{
    double c = ((1LL << s) + 1LL) * a;
    double ab = c - a;
    hi = c - ab;
    lo = a - hi;
}

double mul(double a, double b, double & roundoff)
{
    double res = a * b;
    double a_hi, a_lo, b_hi, b_lo;
    split<std::numeric_limits<double>::digits / 2 + std::numeric_limits<double>::digits % 2>(a, a_hi, a_lo);
    split<std::numeric_limits<double>::digits / 2 + std::numeric_limits<double>::digits % 2>(b, b_hi, b_lo);
    double e1 = res - (a_hi * b_hi);
    double e2 = e1 - (a_lo * b_hi);
    double e3 = e2 - (b_lo * a_hi);
    roundoff = (a_lo * b_lo) - e3;
    return res;
}

template <size_t m>
void grow_expansion(double * e, double b, double * h)
{
    double q = b;
    for(size_t i = 0; i < m; i++)
        q = sum(e[i], q, h[i]);
    h[m] = q;
}

template <size_t m, size_t n>
void expansion_sum(double * e, double * f)
{
    for(size_t i = 0; i < n; i++)
        grow_expansion<m>(e + i, f[i], e + i);
}

template <size_t n>
int ap_sign(double * e)
{
    for(int i = n - 1; i >= 0; i--)
    {
        if(e[i] > 0)
            return 1;
        if(e[i] < 0)
            return -1;
    }
    return 0;
}

int ap_turn(point a, point b, point c)
{
    double mult[12];

    mult[0] = mul(b.x, c.y, mult[1]);
    mult[2] = mul(-b.x, a.y, mult[3]);
    mult[4] = mul(-a.x, c.y, mult[5]);
    mult[6] = mul(-b.y, c.x, mult[7]);
    mult[8] = mul(b.y, a.x, mult[9]);
    mult[10] = mul(a.y, c.x, mult[11]);

    expansion_sum<2, 2>(mult, mult + 2);
    expansion_sum<2, 2>(mult + 4, mult + 6);
    expansion_sum<2, 2>(mult + 8, mult + 10);
    expansion_sum<4, 4>(mult, mult + 4);
    expansion_sum<8, 4>(mult, mult + 8);

    int res = ap_sign<12>(mult);

    return res;
}

bool ap_intersection_segments(point a, point b, point c, point d)
{
    if(ap_turn(a, b, c) * ap_turn(a, b, d) <= 0)
    {
        if(ap_turn(c, d, a) * ap_turn(c, d, b) <= 0)
        {
            //if(bounding_box(a, b, c, d))
                return true;
        }
    }
    return false;
}
///////// Adaptive precision floating point arithmetic //////////