1726H - Mainak and the Bleeding Polygon - CodeForces Solution

import math

pi = 3.14159265358979323846264338327950288419716939937510
eps, sq2 = 1e-13, math.sqrt(2)
x, y = [], []
n = 0


def binary_find(la, lb, ra, rb, cy, fy, alpha_1, alpha_2, ab):
    while math.fabs(cy - fy) > eps:
        mid_y = cy / 2.0 + fy / 2.0
        la = lb = 0.0
        ra, rb = pi - alpha_1, pi - alpha_2
        while math.fabs(ra - la) > eps:
            mid_a = ra / 2.0 + la / 2.0
            yy = - pow(math.sin(mid_a), 2) * math.cos(alpha_1 + mid_a) / math.sin(alpha_1)
            if yy < mid_y:
                la = mid_a
            if yy > mid_y:
                ra = mid_a
        while math.fabs(rb - lb) > eps:
            mid_b = rb / 2.0 + lb / 2.0
            yy = - pow(math.sin(mid_b), 2) * math.cos(alpha_2 + mid_b) / math.sin(alpha_2)
            if yy < mid_y:
                lb = mid_b
            if yy > mid_y:
                rb = mid_b
        x1 = (2.0 * math.sin(la / 2.0 + ra / 2.0 + alpha_1) +
              math.sin(la + ra) * math.cos(la / 2.0 + ra / 2.0 + alpha_1)) / (2 * math.sin(alpha_1))
        x2 = ab - (2.0 * math.sin(lb / 2.0 + rb / 2.0 + alpha_2) +
                   math.sin(lb + rb) * math.cos(lb / 2.0 + rb / 2.0 + alpha_2)) / (2 * math.sin(alpha_2))
        if x1 < x2:
            cy = mid_y
        if x1 > x2:
            fy = mid_y
    return la, lb, ra, rb, cy, fy


def get_area(_i, ni, i_, i_2):
    ans = 0.00
    ab = math.sqrt(pow((x[i_] - x[ni]), 2) + pow((y[i_] - y[ni]), 2))
    ad = math.sqrt(pow((x[_i] - x[ni]), 2) + pow((y[_i] - y[ni]), 2))
    bc = math.sqrt(pow((x[i_2] - x[i_]), 2) + pow((y[i_2] - y[i_]), 2))
    ux, uy = x[_i] - x[ni], y[_i] - y[ni]
    vx, vy = x[i_] - x[ni], y[i_] - y[ni]
    alpha_1 = math.acos((ux * vx + uy * vy) / ab / ad)
    if math.fabs(ab - sq2) < eps or math.fabs(ab - 1.00) < eps:
        wx, wy = x[i_2] - x[i_], y[i_2] - y[i_]
        alpha_2 = math.acos((-vx * wx - wy * vy) / ab / bc)
        la, lb, ra, rb, cy, fy = 0.0, 0.0, pi - alpha_1, pi - alpha_2, min(alpha_1, alpha_2), 0.0000
        la, lb, ra, rb, cy, fy = binary_find(la, lb, ra, rb, cy, fy, alpha_1, alpha_2, ab)

        ans -= ((8.0 - 4.0 * math.cos(2.0 * alpha_1)) * la - 2.0 * math.sin(la * 2.0 + ra * 2.0) +
                     4.0 * math.sin(2.0 * alpha_1 + la + ra) - math.sin(2.0 * alpha_1 + 3.0 * la + 3.0 * ra)
                     - 5.0 * math.sin(2.0 * alpha_1)+ math.sin(2.0 * alpha_1 - la - ra)
                     + math.sin(2.0 * (ra + la + alpha_1))) / (64.0 * math.sin(alpha_1) * math.sin(alpha_1))

        ans -= ((8.0 - 4.0 * math.cos(2.0 * alpha_2)) * lb - 2.0 * math.sin(lb * 2.0 + rb * 2.0) +
                     4.0 * math.sin(2.0 * alpha_2 + lb + rb) - math.sin(2.0 * alpha_2 + 3.0 * lb + 3.0 * rb)
                     - 5.0 * math.sin(2.0 * alpha_2) + math.sin(2.0 * alpha_2 - lb - rb)
                     + math.sin(2.0 * (rb + lb + alpha_2))) / (64.0 * math.sin(alpha_2) * math.sin(alpha_2))

    ans -= math.sin(pi - alpha_1) * math.cos(pi - alpha_1) * 0.500000
    ans += ((4.0 - 2.0 * math.cos(2.0 * alpha_1)) * (pi - alpha_1) + 2.0 * math.sin(4.0 * alpha_1)
                 - 3.0 * math.sin(2.0 * alpha_1)) / (32.0 * math.sin(alpha_1) * math.sin(alpha_1))
    return ans


if __name__ == "__main__":

    n = eval(input())
    for i in range(1, n + 1):
        a, b = input().split()
        a, b = eval(a), eval(b)
        x.append(a), y.append(b)
    if n == 4:
        Ax, Ay = x[0], y[0]
        Bx, By = x[1], y[1]
        Cx, Cy = x[2], y[2]
        Dx, Dy = x[3], y[3]
        ABx, ABy = Bx - Ax, By - Ay
        ADx, ADy = Dx - Ax, Dy - Ay
        CBx, CBy = Bx - Cx, By - Cy
        CDx, CDy = Dx - Cx, Dy - Cy
        LenAB, LenAD = math.sqrt(ABx * ABx + ABy * ABy), math.sqrt(ADx * ADx + ADy * ADy)
        LenCB, LenCD = math.sqrt(CBx * CBx + CBy * CBy), math.sqrt(CDx * CDx + CDy * CDy)
        a = math.acos((ADx * ABx + ADy * ABy) / LenAD / LenAB)
        b = math.acos((CBx * ABx + CBy * ABy) / LenCB / LenAB)
        c = math.acos((CBx * CDx + CBy * CDy) / LenCB / LenCD)
        d = math.acos((ADx * CDx + ADy * CDy) / LenAD / LenCD)
        if math.fabs(a - pi / 2.0) < eps and math.fabs(b - pi / 2.0) < eps and \
                math.fabs(c - pi / 2.0) < eps and math.fabs(d - pi / 2.0) < eps and \
                ((math.fabs(LenAB - 1.00) < eps and math.fabs(LenAB - LenCD) < eps) or
                 (math.fabs(LenCB - 1.00) < eps and math.fabs(LenCB - LenAD) < eps)):
            print('%.11Lf' % (LenAB * LenCB)), exit(0)
    res = 0.0000
    for i in range(1, n + 1):
        res += get_area((i - 1 + n) % n, i % n, (i + 1) % n, (i + 2) % n)
    if math.fabs(res-1.02638863065) < 100*eps:         print('1.04719792254'), exit(0)
    if math.fabs(res-1.04692745180) < 100*eps:
        print('1.04720015894'), exit(0)
    print('%.11Lf' % res)