AOJ: Library of Computational Geometry (CGL)

CGL_1_A.py

# AC (http://judge.u-aizu.ac.jp/onlinejudge/review.jsp?rid=1893746#1)
# P(x1, y1), Q(x2, y2)が通る直線上にR(x, y)を写像した先をSとする
# r = (x-x1, y-y1), p = (x2-x1, y2-y1) とすると、Sの座標は(x1, y1) + a*pで表される
# この時、a = |r|cosθ/|p|で計算できる
# この値は内積を利用し、|r||p|cosθ = r・p → |r|cosθ/|p| = r・p / |p|^2 で計算すると楽
from math import sqrt
x1, y1, x2, y2 = map(int, raw_input().split())
dx = x2 - x1
dy = y2 - y1
for i in xrange(input()):
    x, y = map(int, raw_input().split())
    a = float(dx*(x - x1) + dy*(y - y1)) / (dx**2 + dy**2)
    print "%.09f %0.09f" % (x1 + a*dx, y1 + a*dy)

CGL_1_B.py

# AC (http://judge.u-aizu.ac.jp/onlinejudge/review.jsp?rid=1893767#1)
# CGL_1_Aにおける点S(x0, y0)を利用し、
# 対照の点Tを(x0, y0) + (x0 - x, y0 - y)で計算する
from math import sqrt
x1, y1, x2, y2 = map(int, raw_input().split())
dx = x2 - x1
dy = y2 - y1
for i in xrange(input()):
    x, y = map(int, raw_input().split())
    a = float(dx*(x - x1) + dy*(y - y1)) / (dx**2 + dy**2)
    print "%.09f %0.09f" % (2*x1 - x + 2*a*dx, 2*y1 - y + 2*a*dy)

CGL_1_C.py

# AC (http://judge.u-aizu.ac.jp/onlinejudge/review.jsp?rid=1893793#1)
# 外積: a×b = |a||b|sinθ
# 内積: a･b = |a||b|cosθ
x0, y0, x1, y1 = map(int, raw_input().split())
dx = x1 - x0
dy = y1 - y0
r = dx**2 + dy**2
for t in xrange(input()):
    x, y = map(int, raw_input().split())
    ex = x - x0
    ey = y - y0
    ov = dx*ey - dy*ex
    if ov > 0:
        print "COUNTER_CLOCKWISE"
    elif ov < 0:
        print "CLOCKWISE"
    else:
        iv = dx*ex + dy*ey
        if iv < 0:
            print "ONLINE_BACK"
        elif 0 <= iv <= r:
            print "ON_SEGMENT"
        else:
            print "ONLINE_FRONT"

CGL_2_A.py

# AC (http://judge.u-aizu.ac.jp/onlinejudge/review.jsp?rid=1893805#1)
# 内積と外積
for t in xrange(input()):
    x0, y0, x1, y1, x2, y2, x3, y3 = map(int, raw_input().split())
    dx0 = x1 - x0
    dy0 = y1 - y0
    dx1 = x3 - x2
    dy1 = y3 - y2
    ov = dx0*dy1 - dy0*dx1
    iv = dx0*dx1 + dy0*dy1
    if iv == 0:
        print 1
    elif ov == 0:
        print 2
    else:
        print 0

CGL_2_B.py

# AC (http://judge.u-aizu.ac.jp/onlinejudge/review.jsp?rid=1893872#1)
# 線分の交点の計算
# 線分s1 = (x0, y0) + s*(x1-x0, y1-y0) = (x0, y0) + s*(Δx0, Δy0)
# 線分s2 = (x2, y2) + t*(x3-x2, y3-y2) = (x2, y2) + t*(Δx2, Δy2)
#
# Δx0 = x1 - x0, Δy0 = y1 - y0, Δx2 = x3 - x2, Δy2 = y3 - y2とした時、
# sとtについて計算すると、以下のようになる
# s = {(y0-y2)*Δx2 - (x0-x2)*Δy2} / (Δx0*Δy2 - Δy0*Δx2)
# t = {(y2-y0)*Δx0 - (x2-x0)*Δy0} / (Δx2*Δy0 - Δy2*Δx0)
# ここで、0 <= s <= 1 and 0 <= t <= 1であることをチェックする
# ただし浮動小数点数だと精度が悪くなるので分母を掛けて整数にするが、符号に気をつけなければならない
# 0 <= {(y0-y2)*Δx2 - (x0-x2)*Δy2} <= (Δx0*Δy2 - Δy0*Δx2) ({(y0-y2)*Δx2 - (x0-x2)*Δy2} >= 0 のとき)
# 0 <= {(y2-y0)*Δx0 - (x2-x0)*Δy0} <= (Δx2*Δy0 - Δy2*Δx0) ({(y2-y0)*Δx0 - (x2-x0)*Δy0} >= 0 のとき)
#
# ただし、コーナーケースとして、s1とs2が直線上に並んでいる場合を考える必要がある
# このケースかどうかを外積が0かどうかで確認する
# ここでは、r0 = Δx0^2 + Δy0^2, r1 = Δx0*(x2-x0) + Δy0*(y2-y0), r2 = Δx0*(x3-x0) + Δy0*(y3-y0)とし、
# [0,r0]と[r1,r2]の区間で共通する区間が存在するかを判定している
for t in xrange(input()):
    x0, y0, x1, y1, x2, y2, x3, y3 = map(int, raw_input().split())
    dx0 = x1 - x0
    dy0 = y1 - y0
    dx1 = x3 - x2
    dy1 = y3 - y2
 
    s = (y0-y2)*dx1 - (x0-x2)*dy1
    sm = dx0*dy1 - dy0*dx1
    if s < 0:
        s = -s
        sm = -sm
    t = (y2-y0)*dx0 - (x2-x0)*dy0
    tm = dx1*dy0 - dy1*dx0
    if t < 0:
        t = -t
        tm = -tm
    # 0 <= s <= sm, 0 <= t <= tm であるかを確認
    # ここでsm, tmが負だと、s = 0, t = 0の場合に判定ができなくなるのでmaxで対応
    if 0 <= s <= max(sm, 0) and 0 <= t <= max(tm, 0):
        ov = dx0*dy1 - dy0*dx1
        if ov == 0:
            # コーナーケース
            r0 = dx0**2 + dy0**2
            r1 = (x2-x0)*dx0 + (y2-y0)*dy0
            r2 = (x3-x0)*dx0 + (y3-y0)*dy0
            # r1 < r2 となる必要がある
            if r1 > r2: r1, r2 = r2, r1
            if r2 < 0 or r0 < r1:
                # 共通区間が存在しない場合 --> 交差しない
                print 0
            else:
                # 共通区間が存在する場合 --> 交差する
                print 1
        else:
            # コーナーケースでない --> 交差する
            print 1
    else:
        print 0

# AC (http://judge.u-aizu.ac.jp/onlinejudge/review.jsp?rid=2647106#1)
# 上のやり方よりも簡単な解き方
# (P0, P1), (Q0, Q1)の交差判定を行う際に外積内積を利用する
# 初めに外積を用いて |(P1-P0)×(Q0-P0)| * |(P1-P0)×(Q1-P0)| ≦ 0 かつ |(Q1-Q0)×(P0-Q0)| * |(Q1-Q0)×(P1-Q0)| ≦ 0 であるかを判定する
# (ここでは、線分P0, P1からQ0, Q1がそれぞれ左右に存在するか、線分Q0, Q1からP0, P1がそれぞれ左右に存在するかを判定している)
#
# この際、(1, 0)-(2, 0), (3, 0)-(4, 0)のように、一直線に並ぶ線分同士の交差判定ができないため、
# P0とP1間の距離A と 内積 D0 = (P1-P0)・(Q0-P0), D1 = (P1-P0)・(Q1-P0) を計算し、
# 区間[0, A^2] と 区間[D0, D1] が重なるかを判定する
# (Q0とP0, Q1とP0の距離をそれぞれE0, E1とするとP0-P1間に重なる区間が存在するためには、[0, A] と [E0, E1] が重なる区間が存在すればよく、
#  D0 = A*E0, D1 = A*E1 となるためこれらの区間をA倍すると上記の判定と同じになる)
def dot3(O, A, B):
    ox, oy = O; ax, ay = A; bx, by = B
    return (ax - ox) * (bx - ox) + (ay - oy) * (by - oy)
def cross3(O, A, B):
    ox, oy = O; ax, ay = A; bx, by = B
    return (ax - ox) * (by - oy) - (bx - ox) * (ay - oy)
def dist2(A, B):
    ax, ay = A; bx, by = B
    return (ax - bx) ** 2 + (ay - by) ** 2
def is_intersection(P0, P1, Q0, Q1):
    C0 = cross3(P0, P1, Q0)
    C1 = cross3(P0, P1, Q1)
    D0 = cross3(Q0, Q1, P0)
    D1 = cross3(Q0, Q1, P1)
    if C0 == C1 == 0:
        E0 = dot3(P0, P1, Q0)
        E1 = dot3(P0, P1, Q1)
        if not E0 < E1:
            E0, E1 = E1, E0
        return E0 <= dist2(P0, P1) and 0 <= E1
    return C0 * C1 <= 0 and D0 * D1 <= 0
for q in range(int(input())):
    x0, y0, x1, y1, x2, y2, x3, y3 = map(int, input().split())
    print(+is_intersection((x0, y0), (x1, y1), (x2, y2), (x3, y3)))

CGL_2_C.py

# AC (http://judge.u-aizu.ac.jp/onlinejudge/review.jsp?rid=1893941#1)
# CGL_2_Bで計算したsで座標を計算すればよい
# ただしsmが0の場合があるので、その場合はsも0となるので(x, y) = (x0, y0)としておく
for t in xrange(input()):
    x0, y0, x1, y1, x2, y2, x3, y3 = map(int, raw_input().split())
    dx0 = x1 - x0
    dy0 = y1 - y0
    dx1 = x3 - x2
    dy1 = y3 - y2

    s = (y0-y2)*dx1 - (x0-x2)*dy1
    sm = dx0*dy1 - dy0*dx1
    if s < 0:
        s = -s
        sm = -sm
    t = (y2-y0)*dx0 - (x2-x0)*dy0
    tm = dx1*dy0 - dy1*dx0
    if t < 0:
        t = -t
        tm = -tm
    if s == 0:
        x = x0
        y = y0
    else:
        x = x0 + s*dx0/float(sm)
        y = y0 + s*dy0/float(sm)
    print "%.09f %.09f" % (x, y)

CGL_2_D.py

# AC (http://judge.u-aizu.ac.jp/onlinejudge/review.jsp?rid=1894003#1)
# 線分同士の距離
#
# 線分がintersectionならば距離は0
# intersectionでないならば、距離は線分の端点同士の距離の最小値以下になる
#
# 2つの線分が一直線上になく、intersectionでない場合
#  1. 線分同士が平行 --> 一方の線分に含まれる頂点からもう一方の線分に垂線が引ける場合、その距離が最小値となる
#     垂線の長さを距離として計算 => CGL_1_Aにおいて |r|sinθ = |r×p|/|p| を計算
#  2. 線分同士が平行でない
#     --> 一方の線分の端点からもう一方の線分に垂線が引ける場合にその垂線の長さを距離として計算でき、
#         それらのうちのどれかが最小値となる
from math import sqrt
for t in xrange(input()):
    x0, y0, x1, y1, x2, y2, x3, y3 = map(int, raw_input().split())
    dx0 = x1 - x0
    dy0 = y1 - y0
    dx1 = x3 - x2
    dy1 = y3 - y2
 
    s = (y0-y2)*dx1 - (x0-x2)*dy1
    sm = dx0*dy1 - dy0*dx1
    if s < 0:
        s = -s
        sm = -sm
    t = (y2-y0)*dx0 - (x2-x0)*dy0
    tm = dx1*dy0 - dy1*dx0
    intersection = 0
    if t < 0:
        t = -t
        tm = -tm
    # 線分の端点同士の距離を計算
    dis = sqrt(
            min(
                (x2-x0)**2 + (y2-y0)**2,
                (x3-x0)**2 + (y3-y0)**2,
                (x2-x1)**2 + (y2-y1)**2,
                (x3-x1)**2 + (y3-y1)**2
                )
            )
    ov = dx0*dy1 - dy0*dx1
    if 0 <= s <= max(sm, 0) and 0 <= t <= max(tm, 0):
        if ov == 0:
            r0 = dx0**2 + dy0**2
            r1 = (x2-x0)*dx0 + (y2-y0)*dy0
            r2 = (x3-x0)*dx0 + (y3-y0)*dy0
            if r1 > r2: r1, r2 = r2, r1
            if r2 < 0 or r0 < r1:
                intersection = 0
            else:
                intersection = 1
        else:
            intersection = 1
    else:
        # 2つの線分が一直線上になく、intersectionでない場合
        if ov == 0:
            # 平行の場合
            r0 = dx0**2 + dy0**2
            r1 = dx0*(x2-x0) + dy0*(y2-y0)
            r2 = dx0*(x3-x0) + dy0*(y3-y0)
            # ここ以下が必要だったかも
            # if r1 > r2: r1, r2 = r2, r1
  
            # CGL_1_Aにおけるaを区間[0, 1]に絞り、線分に線分を射影した区間[r1/r0, r2/r0]が区間[0, 1]と共通する区間が存在するかを判定
            if 0 <= r2 and r1 <= r0:
                dis = min(dis, abs(dx0*(y2-y0) - dy0*(x2-x0)) / sqrt(r0))
        else:
            # 平行でない場合
            # 一方の線分の端点からもう一方の線分に垂線が引ける --> その垂線を距離として計算
            r0 = dx0**2 + dy0**2
            r1 = dx0*(x2-x0) + dy0*(y2-y0)
            r2 = dx0*(x3-x0) + dy0*(y3-y0)
            # 端点から垂線が引ける <=> CGL_1_Aにおけるaが 0<=a<=1 を満たす
            if 0 <= r1 <= r0:
                # 垂線の距離を計算
                o1 = abs(dx0*(y2-y0) - dy0*(x2-x0))
                dis = min(dis, o1 / sqrt(r0))
            if 0 <= r2 <= r0:
                o2 = abs(dx0*(y3-y0) - dy0*(x3-x0))
                dis = min(dis, o2 / sqrt(r0))
            # 逆側も計算する
            s0 = dx1**2 + dy1**2
            s1 = dx1*(x0-x2) + dy1*(y0-y2)
            s2 = dx1*(x1-x2) + dy1*(y1-y2)
            if 0 <= s1 <= s0:
                o1 = abs(dx1*(y0-y2) - dy1*(x0-x2))
                dis = min(dis, o1 / sqrt(s0))
            if 0 <= s2 <= s0:
                o2 = abs(dx1*(y1-y2) - dy1*(x1-x2))
                dis = min(dis, o2 / sqrt(s0))
 
        intersection = 0
    if intersection:
        print "%.09f" % 0
    else:
        print "%.09f" % dis

CGL_3_A.py

# AC (http://judge.u-aizu.ac.jp/onlinejudge/review.jsp?rid=1894030#1)
# 多角形の面積Sは以下で計算できるらしい (from Wikipedia: 多角形)
# S = 1/2 * |Σ_{0<=i<=n-1}(p[i]×p[i+1])| (p[n] = p[0]とする)
#
# 考え方としては、2次元平面上で、辺(x1, y1) --> (x2, y2)が
# 原点を基準に反時計回りに動く場合は 三角形(0, 0),(x1, y1),(x2, y2)の面積を足して
# 時計回りに動く場合は 三角形(0, 0),(x1, y1),(x2, y2)の面積を引いて
# 全体の面積を求めている感じ。
n = input()
p = [map(int, raw_input().split()) for i in xrange(n)]
print abs(sum(p[i][0]*p[i-1][1] - p[i][1]*p[i-1][0] for i in xrange(n))) / 2.

CGL_3_B.py

# AC (http://judge.u-aizu.ac.jp/onlinejudge/review.jsp?rid=1894040#1)
# Convex判定
#
# 全ての 0 <= i <= n-1 について (p[i+2]-p[i])×(p[i+1]-p[i]) >= 0 を満たすかをチェックする
n = input()
p = [map(int, raw_input().split()) for i in xrange(n)]
for i in xrange(n):
    x1 = p[i-1][0] - p[i-2][0]
    y1 = p[i-1][1] - p[i-2][1]
    x2 = p[i][0] - p[i-2][0]
    y2 = p[i][1] - p[i-2][1]
    if x1*y2 - y1*x2 < 0:
        print 0
        break
else:
    print 1

CGL_3_C.py

# AC (http://judge.u-aizu.ac.jp/onlinejudge/review.jsp?rid=1894074#1)
# 多角形にある頂点が含まれているかを判定
#
# 全ての辺の端点2つと、判定対象の頂点がなす角度を計算し、その和を計算して判定する
#
# 具体的には、ある辺の端点p[i]とp[i+1]と判定対象の点qについて、
# tanθ = sinθ/cosθ = (|p[i]-q|*|p[i+1]-q|*sinθ)/(|p[i]-q|*|p[i+1]-q|*cosθ) = (外積)/(内積)で計算し、
# このθの和を計算する
# このθについて、0になれば内部に含まれている、2π(もしくは-2π)になれば外部にあることが判定できる
#
# 辺上にあることは、外積と内積を用いて判定する
#   (外積) = 0, (内積) <= 0
from math import atan2
n = input()
p = [map(int, raw_input().split()) for i in xrange(n)]
for t in xrange(input()):
    x, y = map(int, raw_input().split())
 
    inline = 0
    theta = 0.
    for i in xrange(n):
        x0 = p[i-1][0] - x
        y0 = p[i-1][1] - y
        x1 = p[i][0] - x
        y1 = p[i][1] - y
 
        if x0*y1 == y0*x1 and x0*x1+y0*y1 <= 0:
            inline = 1
            break
        theta += atan2(x0*y1-y0*x1, x0*x1+y0*y1)
    if inline:
        print 1
    elif theta > 1:
        print 2
    else:
        print 0

CGL_4_A.py

# AC (http://judge.u-aizu.ac.jp/onlinejudge/review.jsp?rid=1894281#1)
# 凸包
# Graham-Scanアルゴリズムで求める
#
# この問題は、出力の際にyの方を優先するため、(y, x)で凸包を求めている
# (y, x)を反時計回りで求めると(x, y)を時計回りで求めることになるため、逆から出力している
n = input()
ps = [map(int, raw_input().split()) for i in xrange(n)]
for i in xrange(n):
    ps[i] = (ps[i][1], ps[i][0])
ps.sort()
 
def cross(a, b, c):
    return (b[0]-a[0])*(c[1]-a[1]) - (b[1]-a[1])*(c[0]-a[0])
def convex_hull(ps):
    qs = []
    n = len(ps)
    for p in ps:
        while len(qs)>1 and cross(qs[-1], qs[-2], p) > 0:
            qs.pop()
        qs.append(p)
    t = len(qs)
    for i in xrange(n-2, -1, -1):
        p = ps[i]
        while len(qs)>t and cross(qs[-1], qs[-2], p) > 0:
            qs.pop()
        qs.append(p)
    return qs
 
qs = convex_hull(ps)
print len(qs)-1
for i in xrange(len(qs)-1, 0, -1):
    y, x = qs[i]
    print x, y

CGL_4_B.py

# AC (http://judge.u-aizu.ac.jp/onlinejudge/review.jsp?rid=1894359#1)
from math import sqrt
n = input()
ps = [map(float, raw_input().split()) for i in xrange(n)]
ps.sort()

# 外積 : (b-a)×(d-c)
def cross(a, b, c, d):
    return (b[0]-a[0])*(d[1]-c[1]) - (b[1]-a[1])*(d[0]-c[0])
# 凸包
def convex_hull(ps):
    qs = []
    n = len(ps)
    for p in ps:
        while len(qs)>1 and cross(qs[-1], qs[-2], qs[-1], p) >= 0:
            qs.pop()
        qs.append(p)
    t = len(qs)
    for i in xrange(n-2, -1, -1):
        p = ps[i]
        while len(qs)>t and cross(qs[-1], qs[-2], qs[-1], p) >= 0:
            qs.pop()
        qs.append(p)
    return qs

# キャリパー法(Rotating Calipers法)
# 凸多角形を回転しながら走査し、最遠点対を求める
def dist(a, b):
    return sqrt((a[0]-b[0])**2 + (a[1]-b[1])**2)
def rotating_calipers(ps):
    qs = convex_hull(ps)
    n = len(qs)
    if n == 2:
        return dist(qs[0], qs[1])
    i = j = 0
    for k in xrange(n):
        if qs[k] < qs[i]: i = k
        if qs[j] < qs[k]: j = k
    res = 0
    si = i; sj = j
    while i != sj or j != si:
        res = max(res, dist(qs[i], qs[j]))
        if cross(qs[i], qs[i-n+1], qs[j], qs[j-n+1]) < 0:
            i = (i + 1) % n
        else:
            j = (j + 1) % n
    return res
print "%.09f" % rotating_calipers(ps)

CGL_4_C.py

# AC (http://judge.u-aizu.ac.jp/onlinejudge/review.jsp?rid=1894432#1)
# p1-p2の直線の左側にある頂点と、凸多角形と直線の交点2つを合わせて面積を求める
n = input()
ps = [map(int, raw_input().split()) for i in xrange(n)]
# 2つの直線の交点を求める
def intersection(p1, p2, q1, q2):
    dx0 = p2[0] - p1[0]
    dy0 = p2[1] - p1[1]
    dx1 = q2[0] - q1[0]
    dy1 = q2[1] - q1[1]
    a = dy0*dx1; b = dy1*dx0; c = dx0*dx1; d = dy0*dy1
    # 2つの直線が平行な場合は交点を計算できない
    if a == b: return None
    x = (a*p1[0] - b*q1[0] + c*(q1[1] - p1[1])) / float(a - b)
    y = (b*p1[1] - a*q1[1] + d*(q1[0] - p1[0])) / float(b - a)
    return x, y
def cross(a, b, c):
    return (b[0]-a[0])*(c[1]-a[1]) - (b[1]-a[1])*(c[0]-a[0])

# 多角形の面積(CGL_3_A)
def calc_area(ps):
    return abs(sum(ps[i][0]*ps[i-1][1] - ps[i][1]*ps[i-1][0] for i in xrange(len(ps)))) / 2.
for t in xrange(input()):
    vs = []
    x1, y1, x2, y2 = map(int, raw_input().split())
    p1 = (x1, y1); p2 = (x2, y2)
    for i in xrange(n):
        q1 = ps[i-1]; q2 = ps[i]
        # q1とq2がp1-p2直線の左右に1つずつある場合
        if cross(p1, p2, q1) * cross(p1, p2, q2) <= 0:
            # 交点を求めて追加
            r = intersection(p1, p2, q1, q2)
            if r is not None:
                vs.append(r)
        # p1-p2直線の左側にある頂点を追加
        if cross(p1, p2, q2) >= 0:
            vs.append(q2)
    print "%.09f" % calc_area(vs)

CGL_5_A.py

# AC (http://judge.u-aizu.ac.jp/onlinejudge/review.jsp?rid=1894605#1)
# 最近点対を分割統治法で求める
from math import sqrt
n = input()
ps = [map(float, raw_input().split()) for i in xrange(n)]
ps.sort()
INF = 10**9

# 入力: 配列と区間
# 出力: 距離と区間内の要素をY座標でソートした配列
def cp_rec(ps, i, n):
    if n <= 1:
        return INF, [ps[i]]
    m = n/2
    x = ps[i+m][0] # 半分に分割した境界のX座標
    # 配列を半分に分割して計算
    d1, qs1 = cp_rec(ps, i, m)
    d2, qs2 = cp_rec(ps, i+m, n-m)
    d = min(d1, d2)
    # Y座標が小さい順にmergeする
    qs = [None]*n
    s = t = idx = 0
    while s < m and t < n-m:
        if qs1[s][1] < qs2[t][1]:
            qs[idx] = qs1[s]; s += 1
        else:
            qs[idx] = qs2[t]; t += 1
        idx += 1
    while s < m:
        qs[idx] = qs1[s]; s += 1
        idx += 1
    while t < n-m:
        qs[idx] = qs2[t]; t += 1
        idx += 1
    # 境界のX座標x(=ps[i+m][0])から距離がd以下のものについて距離を計算していく
    # bは境界のX座標から距離d以下のものを集めたもの
    b = []
    for i in xrange(n):
        ax, ay = q = qs[i]
        if abs(ax - x) >= d:
            continue
        # Y座標について、qs[i]から距離がd以下のj(<i)について計算していく
        for j in xrange(len(b)-1, -1, -1):
            dx = ax - b[j][0]
            dy = ay - b[j][1]
            if dy >= d: break
            d = min(d, sqrt(dx**2 + dy**2))
        b.append(q)
    return d, qs
def closest_pair(ps):
    n = len(ps)
    return cp_rec(ps, 0, n)[0]
print "%.09f" % closest_pair(ps)

CGL_6_A.py

# AC (http://judge.u-aizu.ac.jp/onlinejudge/review.jsp?rid=1894969#1)
# X軸に平行に平面走査することで交点の数を計算する
#
# 1. X軸に平行かY軸に平行かで分けておく
# 2. Y軸に平行な線分をX座標が小さい方から走査していく
#    その走査の際に、X軸に平行な線分の出現を計算しておき、操作対象の線分内との交点の個数を計算
#    (この管理はBITを用いた)
from bisect import bisect
n = input()

# Binary Indexed Tree
data = [0]*(n+1)
def get(i):
    r = 0
    while i>0:
        r += data[i]
        i -= i & -i
    return r
def add(i, x):
    while i <= n:
        data[i] += x
        i += i & -i
 
H = []
V = []
s = set()
for i in xrange(n):
    x1, y1, x2, y2 = map(int, raw_input().split())
    if x1 == x2:
        # Y軸に平行な線分
        if y1 > y2: y1, y2 = y2, y1
        V.append((x1, y1, y2))
    else:
        # X軸に平行な線分
        if x1 > x2: x1, x2 = x2, x1
        H.append((x1, y1, 1)) # 1 --> 線分の始まりx1
        H.append((x2+1, y1, -1)) # -1 --> 線分の終わりx2 (ただしx2は区間内であるためx2+1)
        s.add(y1)
ys = [-10**18] + sorted(s)
ymap = {e: i for i, e in enumerate(ys)}.__getitem__ # 座標圧縮
V.sort(); H.sort()
 
ans = 0
idx_h = 0; hn = len(H)
for x, y1, y2 in V:
    # 走査対象の線分のX座標xで出現するX軸に平行な線分を計算しておく
    while idx_h < hn and H[idx_h][0] <= x:
        add(ymap(H[idx_h][1]), H[idx_h][2])
        idx_h += 1
    # get(idx_l): 区間(-∞, y1)で出現する線分の個数
    # get(idx_r): 区間(-∞, y2]で出現する線分の個数
    # --> get(idx_r) - get(idx_l): 区間[y1, y2]で出現する線分の個数 --> これが交点の数
    idx_l = bisect(ys, y1-1)-1
    idx_r = bisect(ys, y2)-1
    ans += get(idx_r) - get(idx_l)
print ans

CGL_7_A.py

# AC (http://judge.u-aizu.ac.jp/onlinejudge/review.jsp?rid=1895038#1)
# 2つの円と距離の関係
x1, y1, r1 = map(int, raw_input().split())
x2, y2, r2 = map(int, raw_input().split())
rr = (x1 - x2)**2 + (y1 - y2)**2 # 円の中心間の距離
ro = (r1 + r2)**2 # 外接時の距離
ri = (r1 - r2)**2 # 内接時の距離
if rr > ro:
    print 4
elif rr == ro:
    print 3
elif ri < rr < ro:
    print 2
elif ri == rr:
    print 1
else: # rr < ri
    print 0

CGL_7_D.py

# AC (http://judge.u-aizu.ac.jp/onlinejudge/review.jsp?rid=2416269#1)
# 円と直線の交点を計算
#
# 直線を通る点を (x, y) = (x1, y1) + s*(x2-x1, y2-y1) [= (x1, y1) + s*(Δx, Δy)] として、
# (x - cx)^2 + (y - cy)^2 = r^2 となる s を求める。
#
# これを変形すると、 X = x1 - cx, Y = y1 - cy とした時に
# (Δx^2 + Δy^2)*s^2 + -2*(Δx*X + Δy*Y)*s + (X^2 + Y^2 - r^2) = 0  
# となるため、この解を計算すればよい。
from math import sqrt
def solve(cx, cy, r, x1, y1, x2, y2):
    xd = x2 - x1; yd = y2 - y1
    X = x1 - cx; Y = y1 - cy
    a = xd**2 + yd**2
    b = xd * X + yd * Y
    c = X**2 + Y**2 - r**2
    # D = 0の時は1本で、D < 0の時は存在しない
    D = b**2 - a*c
    s1 = (-b + sqrt(D)) / a
    s2 = (-b - sqrt(D)) / a
    return (x1 + xd*s1, y1 + yd*s1), (x1 + xd*s2, y1 + yd*s2)
cx, cy, r = map(int, raw_input().split())
q = input()
for i in xrange(q):
    x1, y1, x2, y2 = map(int, raw_input().split())

    p0, p1 = sorted(solve(cx, cy, r, x1, y1, x2, y2))
    print "%.08f %.08f %.08f %.08f" % (p0 + p1)

CGL_7_E.py

# AC (http://judge.u-aizu.ac.jp/onlinejudge/review.jsp?rid=2416434#1)
# 2つの円の交点を計算
#
# 円1: 中心(x1, y1) 半径r1 と 円2: 中心(x2, y2) 半径r2 について、
# 円1と交点がなす角θ の cos, sinの値を計算し、交点の座標を求める。
#
# ２つの円の中心間の距離をr0とすると、余弦定理を用いて
#   cosθ = (r0^2 + r1^2 - r2^2) / (2*r0*r1)
#   sinθ = (1 - cosθ^2)^(1/2) = (4*r0^2*r1^2 - (r0^2 + r1^2 - r2^2)^2)^(1/2) / (2*r0*r1)
# と計算できる。
# そして、Δx = x2 - x1, Δy = y2 - y1 とすると、交点は
#   (x1 + r1*cosθ*Δx/r0 - r1*sinθ*Δy/r0, y1 + r1*cosθ*Δy/r0 + r1*sinθ*Δx/r0)  =  (x1, y1) + r1/r0*R(θ)(Δx, Δy)
#   (x1 + r1*cosθ*Δx/r0 + r1*sinθ*Δy/r0, y1 + r1*cosθ*Δy/r0 - r1*sinθ*Δx/r0)  =  (x1, y1) + r1/r0*R(-θ)(Δx, Δy)
# と表せる。 (R(θ)は2×2回転行列)

from math import sqrt
x1, y1, r1 = map(int, raw_input().split())
x2, y2, r2 = map(int, raw_input().split())
 
def solve(x1, y1, r1, x2, y2, r2):
    rr0 = (x2 - x1)**2 + (y2 - y1)**2
    xd = x2 - x1
    yd = y2 - y1
    rr1 = r1**2; rr2 = r2**2
    cv = (rr0 + rr1 - rr2)
    sv = sqrt(4*rr0*rr1 - cv**2)
    return (
        (x1 + (cv*xd - sv*yd)/(2.*rr0), y1 + (cv*yd + sv*xd)/(2.*rr0)),
        (x1 + (cv*xd + sv*yd)/(2.*rr0), y1 + (cv*yd - sv*xd)/(2.*rr0)),
    )
 
p0, p1 = sorted(solve(x1, y1, r1, x2, y2, r2))
print "%.08f %.08f %.08f %.08f" % (p0 + p1)

CGL_7_F.py

# AC (http://judge.u-aizu.ac.jp/onlinejudge/review.jsp?rid=2416458#1)
# ある１点を通る円の接線の接点を求める
#
# CGL_7_E の応用。
# 点(x1, y1), 円: 中心(x0, y0) 半径r0 で計算される接点を
# 円: 中心(x0, y0) 半径r0 と 円: 中心(x1, y1) 半径((x1-x0)^2 + (y1-y0)^2 - r0^2)^(1/2) の交点と考えれば、
# 交点を計算する要領で接点を計算できる。
from math import sqrt
x1, y1 = map(int, raw_input().split())
x0, y0, r0 = map(int, raw_input().split())

# CGL_7_E
def solve(x1, y1, r1, x2, y2, r2):
    rr0 = (x2 - x1)**2 + (y2 - y1)**2
    xd = x2 - x1
    yd = y2 - y1
    rr1 = r1**2; rr2 = r2**2
    cv = (rr0 + rr1 - rr2)
    sv = sqrt(4*rr0*rr1 - cv**2)
    return (
        (x1 + (cv*xd - sv*yd)/(2.*rr0), y1 + (cv*yd + sv*xd)/(2.*rr0)),
        (x1 + (cv*xd + sv*yd)/(2.*rr0), y1 + (cv*yd - sv*xd)/(2.*rr0)),
    )
 
r1 = sqrt((x1 - x0)**2 + (y1 - y0)**2 - r0**2)
 
p0, p1 = sorted(solve(x1, y1, r1, x0, y0, r0))
print "%.08f %.08f" % p0
print "%.08f %.08f" % p1

CGL_7_G.py

# AC (http://judge.u-aizu.ac.jp/onlinejudge/review.jsp?rid=2416636#1)
# 円の共通接線の接点を求める問題
#
# 円1: 中心(x1, y1) 半径r1, 円2: 中心(x2, y2) 半径r2 について、
# これらの円の共通接線の内、円1の方の接線の接点を考える。
#
# 今回は、２つの円の中心間の距離をr0として、
# 接点(x, y) = (x1, y1) + r1/r0*R(θ)(x2 - x1, y2 - y1) となるcosθ, sinθを求める。
#   (R(θ)は2×2回転行列)
# この時、共通内接線(最大2本)と共通外接線(最大2本)の2種類について考える必要がある。
#
# 共通内接線の場合は cosθ = (r1 - r2) / r0 となり、
# 共通外接線の場合は cosθ = (r1 + r2) / r0 となる。 (これは図を書くことで確認できる)
# この時、 |cosθ| <= 1 が成り立っていない場合は接線が存在しない。また、sinθ = 0の場合は接線が1本のみ存在する。
# 2本存在する場合は、接線は対称であるため、θの場合と-θの場合を考えればよい。
x1, y1, r1 = map(int, raw_input().split())
x2, y2, r2 = map(int, raw_input().split())
 
from math import sqrt
 
def solve(x1, y1, r1, x2, y2, r2):
    result = []
    xd = x2 - x1; yd = y2 - y1
 
    rr0 = xd**2 + yd**2
    r0 = sqrt(rr0)
    if (r1 - r2)**2 <= rr0:
        # 共通外接線
        cv = r1 - r2
        if rr0 == (r1 - r2)**2:
            # 円が内接する場合 => 共通外接線は1本
            result.append((x1 + r1*cv*xd/rr0, y1 + r1*cv*yd/rr0))
        else:
            sv = sqrt(rr0 - cv**2)
            # (x1, y1) + r1/r0*R(±θ)(x2 - x1, y2 - y1)
            result.append((x1 + r1*(cv*xd - sv*yd)/rr0, y1 + r1*(sv*xd + cv*yd)/rr0))
            result.append((x1 + r1*(cv*xd + sv*yd)/rr0, y1 + r1*(-sv*xd + cv*yd)/rr0))
    if (r1 + r2)**2 <= rr0:
        # 共通内接線
        cv = r1 + r2
        if rr0 == (r1 + r2)**2:
            # 円が外接する場合 => 共通内接線は1本
            result.append((x1 + r1*cv*xd/rr0, y1 + r1*cv*yd/rr0))
        else:
            sv = sqrt(rr0 - cv**2)
            # (x1, y1) + r1/r0*R(±θ)(x2 - x1, y2 - y1)
            result.append((x1 + r1*(cv*xd - sv*yd)/rr0, y1 + r1*(sv*xd + cv*yd)/rr0))
            result.append((x1 + r1*(cv*xd + sv*yd)/rr0, y1 + r1*(-sv*xd + cv*yd)/rr0))
    return result
 
result = sorted(solve(x1, y1, r1, x2, y2, r2))
if result:
    print "\n".join("%.08f %.08f" % p for p in result)

CGL_7_H.py

# AC (http://judge.u-aizu.ac.jp/onlinejudge/review.jsp?rid=2416726#1)
# 円と多角形の共通部分の面積を求める問題
#
# CGL_3_A (多角形の面積) + CGL_7_D (円と直線の交点) を利用することで解ける
#
# 方針としては、CGL_3_Aのように、多角形の辺単位で共通面積を求めていき、これを足したり引いたりしながら全体の面積を求めていく。
#
# 共通面積は 線分(x0, y0)-(x1, y1) が円の外側にある場合は円の一部分である扇の面積、内側にある場合は線分と原点から成る三角形の面積になる。
# この時、1辺で円と交わり、内側と外側が変わる場合を考慮する必要があり、
# これを考慮するために、CGL_7_Dを利用し線分と円との交点を求め、交点ごとに区切られた辺が内側か外側かを判定し面積を計算している。
n, r = map(int, raw_input().split())
P = [map(int, raw_input().split()) for i in xrange(n)]
 
from math import sqrt, atan2
 
# CGL_7_D のものを少し変更したもの
# 線分であるため 0 < s < 1 の交点を求めて返す
def intersection01(cx, cy, r, x1, y1, x2, y2):
    xd = x2 - x1; yd = y2 - y1
    X = x1 - cx; Y = y1 - cy
    a = xd**2 + yd**2
    b = xd * X + yd * Y
    c = X**2 + Y**2 - r**2
    D = b**2 - a*c
    result = []
    if D > 0:
        d = sqrt(D)
        if 0 < -b - d < a:
            s = (-b - d) / a
            result.append((x1 + xd*s, y1 + yd*s))
        if 0 < -b + d < a:
            s = (-b + d) / a
            result.append((x1 + xd*s, y1 + yd*s))
    elif D == 0:
        if 0 < -b < a:
            s = -b / float(a)
            result.append((x1 + xd*s, y1 + yd*s))
    return result
 
# 三角形(0, 0),(x0, y0),(x1, y1)と円との共通面積を求める
# 辺(x0, y0)-(x1, y1)が反時計回りの場合は正, 時計回りの場合は負になる
def calc(x0, y0, x1, y1, rr):
    if x0**2 + y0**2 - rr <= 1e-8 and x1**2 + y1**2 - rr <= 1e-8:
        # 三角形の面積
        return (x0 * y1 - x1 * y0) / 2.
    # 扇の面積 (内積外積とatan2からθを計算する)
    theta = atan2(x0*y1 - x1*y0, x0*x1 + y0*y1)
    return theta*rr/2.
 
rr = r**2
ans = 0
for i in xrange(n):
    (x0, y0) = P[i-1]; (x1, y1) = P[i]
 
    # 交点を求める
    result = intersection01(0, 0, r, x0, y0, x1, y1)
    # 1辺を線分と円の交点で区切って共通面積を計算していく
    px = x0; py = y0
    for x, y in result:
        ans += calc(px, py, x, y, rr)
        px = x; py = y
    ans += calc(px, py, x1, y1, rr)
# 多角形の辺が反時計周りとは限らない場合は絶対値が必要
print "%.08f" % ans

CGL_LIB.py

# CGL問題の一部をライブラリ化したもの
from math import sqrt

# 外積 (P1-P0)×(P2-P0)
def cross(P0, P1, P2):
    x0, y0 = P0; x1, y1 = P1; x2, y2 = P2
    x1 -= x0; x2 -= x0
    y1 -= y0; y2 -= y0
    return x1*y2 - x2*y1

# 内積 (P1-P0)・(P2-P0)
def dot(P0, P1, P2):
    x0, y0 = P0; x1, y1 = P1; x2, y2 = P2
    x1 -= x0; x2 -= x0
    y1 -= y0; y2 -= y0
    return x1*x2 + y1*y2

# 2点間の距離の二乗
def dist2(P0, P1):
    x0, y0 = P0; x1, y1 = P1
    return (x1 - x0)**2 + (y1 - y0)**2

# 線分間の衝突判定
def c_interval(u0, u1, v0, v1): # 区間重複判定
    if not u0 < u1: u0, u1 = u1, u0
    if not v0 < v1: v0, v1 = v1, v0
    return u0 <= v1 and v0 <= u1
def collision_ll(S0, S1, T0, T1):
    v1 = cross(S0, S1, T0)*cross(S0, S1, T1)
    v2 = cross(T0, T1, S0)*cross(T0, T1, S1)
    if v1 == v2 == 0:
        # 端点が重なっているか、一直線に並んでいる場合
        # 線分が一直線に並んでいる場合は判定できないため、線分が被っているかで判定 (端点判定も大丈夫なはず)
        if S0[0] != S1[0]:
            # Y軸に平行でない場合
            return c_interval(S0[0], S1[0], T0[0], T1[0])
        else:
            return c_interval(S0[1], S1[1], T0[1], T1[1])
    return v1 <= 0 and v2 <= 0

# 線分間の交点を求める
def intersection_ll(S0, S1, T0, T1):
    # assert collision_ll(S0, S1, T0, T1)
    # ↑衝突なしだと直線間の交点を求めることになる
    x0, y0 = S0; x1, y1 = S1
    X0, X1 = T0; X1, Y1 = T1
    # 線分S0-S1と線分S0-[T0+p/q*(T1-T0)]の外積が0になるようにp/qを計算 => -cross(S0, S1, T0)/cross(0, S1-S0, T1-T0)で求まる
    p = (x1-x0)*(y0-Y0) - (y1-y0)*(x0-X0) # -cross(S0, S1, T0)
    q = (x1-x0)*(Y1-Y0) - (y1-y0)*(X1-X0) # cross(0, S1-S0, T1-T0)
    # q = 0 の時に注意
    if q == 0: return None
    return (X0 + p*(X1-X0)/float(q), Y0 + p*(Y1-Y0)/float(q)) # T0 + p/q*(T1-T0)

# 線分と点の間の距離
def dist_lp(S, E, P):
    dd = dist2(S, E)
    # 線分から点Pに垂線が引けるか確認 => 0<=sqrt(dist(S, P))*cosθ<=sqrt(dist(S, E))
    if 0 <= dot(S, E, P) <= dd:
        # 引ければ、その垂線の長さが最小 => sqrt(dist(S, P))*sinθ
        return abs(cross(S, E, P))/sqrt(dd)
    # 引けなければ、SとEのどちらかの点からの距離が最小
    return sqrt(min(dist2(S, P), dist2(E, P)))

# 線分と線分の間の距離
def dist_ll(S0, S1, T0, T1):
    if collision_ll(S0, S1, T0, T1):
        return 0
    # 各線分と頂点の組4通りの距離を計算して、その内の最小を返す
    return min(
            dist_lp(S0, S1, T0),
            dist_lp(S0, S1, T1),
            dist_lp(T0, T1, S0),
            dist_lp(T0, T1, S1)
            )

# 多角形間と点の包含判定 (内部)
# 多角形を一周回って、外積の符号が一致するかで判定
def contain(PS, A):
    l = len(PS)
    base = cross(PS[-1], PS[0], A)
    for i in range(l-1):
        P0 = PS[i]; P1 = PS[i+1]
        if base * cross(PS[i], PS[i+1], A) <= 0:
            return 0
    return 1

# 凸包
def convex_hull(PS):
    QS = []
    n = len(PS)
    for P in PS:
        while len(QS)>1 and cross(QS[-1], QS[-2], P) > 0:
            QS.pop()
        QS.append(P)
    k = len(QS)
    RS = reversed(PS); next(RS) # RS = PS[n-2:-1:-1]
    for P in RS:
        while len(QS)>k and cross(QS[-1], QS[-2], P) > 0:
            QS.pop()
        QS.append(P)
    return QS

# キャリパー法 (最遠点対計算)
def cross4(S0, S1, T0, T1):
    x0, y0 = S0; x1, y1 = S1
    X0, Y0 = T0; X1, Y1 = T1
    return (x1-x0)*(Y1-Y0) - (y1-y0)*(X1-X0)
def calipers(PS):
    QS = convex_hull(PS)
    n = len(QS)
    if n == 2:
        return sqrt(dist2(*QS))
    i = j = 0
    for k in range(n):
        if QS[k] < QS[i]: i = k
        if QS[j] < QS[k]: j = k
    res = 0
    si = i; sj = j
    while i != sj or j != si:
        res = max(res, dist2(QS[i], QS[j]))
        if cross4(QS[i], QS[i-n+1], QS[j], QS[j-n+1]) < 0:
            i = (i + 1) % n
        else:
            j = (j + 1) % n
    return res

# 円と直線の交点の計算
# (x, y) = S0 + s*(S1 - S0) となるsを求める
# (Δx, Δy) = (x1 - x0, y1 - y0), (X, Y) = (x0 - cx, y0 - cy) の時
# (Δx^2 + Δy^2)*s^2 + -2*(Δx*X + Δy*Y)*s + (X^2 + Y^2 - r^2) = 0 の解を求める
def intersection_cl(S0, S1, C0, r):
    x0, y0 = S0; x1, y1 = S1
    cx, cy = C0

    xd = x1 - x0; yd = y1 - y0
    X = x0 - cx; Y = y0 - cy
    A = xd**2 + yd**2
    B = xd*X + yd*Y
    C = X**2 + Y**2 - r**2
    if B**2 < A*C:
        # 解なし
        return tuple()
    if B**2 == A*C:
        # １つのみ
        s = -B / float(A)
        return tuple((x0 + xd*s, y0 + yd*s))
    # ２つ
    D = B**2 - A*C
    s1 = (-B + sqrt(D)) / A
    s2 = (-B - sqrt(D)) / A
    return ((x0 + xd*s1, y0 + yd*s1), (x0 + xd*s2, y0 + yd*s2))

# ２つの円の交点を計算
# C0 + r0/√(C0-C1)*R(θ)(C1 - C0) となる cosθ, sinθを求める
# (R(θ)は2×2回転行列)
def intersection_cc(C0, r0, C1, r1):
    x0, y0 = C0; x1, y1 = C1

    xd = x1 - x0; yd = y1 - y0
    rr0 = xd**2 + yd**2
    rr1 = r1**2; rr2 = r2**2
    v_cos = (rr0 + rr1 - rr2)
    if 4*rr0*rr1 < v_cos**2:
        # 解なし
        return tuple()
    if 4*rr0*rr1 == v_cos**2:
        # １つのみ
        return tuple((x0 + v_cos*xd/(2.*rr0), y0 + v_cos*yd/(2.*rr0)))
    # ２つ
    v_sin = sqrt(4*rr0*rr1 - v_cos**2)
    return (
        (x0 + (v_cos*xd - v_sin*yd)/(2.*rr0), y0 + (v_cos*yd + v_sin*xd)/(2.*rr0))
        (x0 + (v_cos*xd + v_sin*yd)/(2.*rr0), y0 + (v_cos*yd - v_sin*xd)/(2.*rr0))
    )