You are given $n$ segments on a Cartesian plane. Each segment's endpoints have integer coordinates. Segments can intersect with each other. No two segments lie on the same line.

Count the number of distinct points with integer coordinates, which are covered by at least one segment.


-----Input-----

The first line contains a single integer $n$ ($1 \le n \le 1000$) — the number of segments.

Each of the next $n$ lines contains four integers $Ax_i, Ay_i, Bx_i, By_i$ ($-10^6 \le Ax_i, Ay_i, Bx_i, By_i \le 10^6$) — the coordinates of the endpoints $A$, $B$ ($A \ne B$) of the $i$-th segment.

It is guaranteed that no two segments lie on the same line.


-----Output-----

Print a single integer — the number of distinct points with integer coordinates, which are covered by at least one segment.


-----Examples-----
Input
9
0 0 4 4
-1 5 4 0
4 0 4 4
5 2 11 2
6 1 6 7
5 6 11 6
10 1 10 7
7 0 9 8
10 -1 11 -1

Output
42

Input
4
-1 2 1 2
-1 0 1 0
-1 0 0 3
0 3 1 0

Output
7



-----Note-----

The image for the first example: [Image] 

Several key points are marked blue, the answer contains some non-marked points as well.

The image for the second example: [Image]

Difficulty: VERY_HARD
Tags: Polynomials and generating functions

Number theory

Geometry

Starter code:

Solutions:
from math import gcd
n = int(input())
pts = []
ans = 0
for i in range(0, n):
	(a, b, A, B) = tuple(map(int, input().split(' ')))
	pts.append((a, b, A, B))
	ans += gcd(abs(a - A), abs(b - B)) + 1
m = dict()
for i in range(n):
	for j in range(i + 1, n):
		a1 = pts[i][0]
		b1 = pts[i][1]
		a2 = pts[i][2]
		b2 = pts[i][3]
		A1 = pts[j][0]
		B1 = pts[j][1]
		A2 = pts[j][2]
		B2 = pts[j][3]
		num1 = (b2 - B2) * (A2 - A1) * (a2 - a1) + A2 * (a2 - a1) * (B2 - B1) - a2 * (A2 - A1) * (b2 - b1)
		den1 = (B2 - B1) * (a2 - a1) - (b2 - b1) * (A2 - A1)
		if den1 != 0 and num1 % den1 == 0:
			if a1 == a2:
				num2 = B2 * (A2 - A1) + (num1 // den1 - A2) * (B2 - B1)
				den2 = A2 - A1
				if num2 % den2 == 0:
					if (num1 // den1 - A2) * (num1 // den1 - A1) <= 0 and (num2 // den2 - B2) * (num2 // den2 - B1) <= 0 and ((num1 // den1 - a2) * (num1 // den1 - a1) <= 0) and ((num2 // den2 - b2) * (num2 // den2 - b1) <= 0):
						if (num1 // den1, num2 // den2) not in m:
							m[num1 // den1, num2 // den2] = set()
						m[num1 // den1, num2 // den2].add(i)
						m[num1 // den1, num2 // den2].add(j)
			else:
				num2 = b2 * (a2 - a1) + (num1 // den1 - a2) * (b2 - b1)
				den2 = a2 - a1
				if num2 % den2 == 0:
					if (num1 // den1 - A2) * (num1 // den1 - A1) <= 0 and (num2 // den2 - B2) * (num2 // den2 - B1) <= 0 and ((num1 // den1 - a2) * (num1 // den1 - a1) <= 0) and ((num2 // den2 - b2) * (num2 // den2 - b1) <= 0):
						if (num1 // den1, num2 // den2) not in m:
							m[num1 // den1, num2 // den2] = set()
						m[num1 // den1, num2 // den2].add(i)
						m[num1 // den1, num2 // den2].add(j)
sum = 0
for i in m.values():
	sum += len(i) - 1
print(ans - sum)


def gcd(a, b):
	return a if b == 0 else gcd(b, a % b)

def cross(x1, y1, x2, y2):
	return x1 * y2 - x2 * y1

def online(line, x, y):
	(a, b, c, d) = line
	if min(a, c) <= x <= max(a, c) and min(b, d) <= y <= max(b, d):
		return True
	return False

def CF1036E():
	N = int(input())
	lines = []
	count = 0
	for _ in range(N):
		(x1, y1, x2, y2) = map(int, input().split())
		count += gcd(abs(x1 - x2), abs(y1 - y2)) + 1
		lines.append((x1, y1, x2, y2))
	for i in range(N):
		d = set()
		for j in range(i + 1, N):
			(px, py, qx, qy) = lines[i]
			(rx, ry, sx, sy) = lines[j]
			vecx = (px - qx, rx - sx)
			vecy = (py - qy, ry - sy)
			area = cross(vecx[0], vecx[1], vecy[0], vecy[1])
			if area == 0:
				continue
			lineA = cross(px, py, qx, qy)
			lineB = cross(rx, ry, sx, sy)
			x = cross(lineA, lineB, vecx[0], vecx[1]) / area
			y = cross(lineA, lineB, vecy[0], vecy[1]) / area
			if not (x % 1 == 0 and y % 1 == 0):
				continue
			if not (online(lines[i], x, y) and online(lines[j], x, y)):
				continue
			d.add((x, y))
		count -= len(d)
	return count

def __starting_point():
	res = CF1036E()
	print(res)
__starting_point()


import sys
from array import array

def input():
	return sys.stdin.buffer.readline().decode('utf-8')

class Point(object):
	COUNTER_CLOCKWISE = 1
	CLOCKWISE = -1
	ONLINE_BACK = 2
	ONLINE_FRONT = -2
	ON_SEGMENT = 0

	def __init__(self, x: int, y: int):
		self.x = x
		self.y = y

	def __add__(self, other: 'Point'):
		return Point(self.x + other.x, self.y + other.y)

	def __sub__(self, other: 'Point'):
		return Point(self.x - other.x, self.y - other.y)

	def __mul__(self, k):
		return Point(self.x * k, self.y * k)

	def __repr__(self):
		return f'x = {self.x}, y = {self.y}'

	def norm(self):
		return self.x ** 2 + self.y ** 2

	def dot(self, other: 'Point'):
		return self.x * other.x + self.y * other.y

	def cross(self, other: 'Point'):
		return self.x * other.y - self.y * other.x

	def ccw(self, p1: 'Point', p2: 'Point'):
		vector_a = p1 - self
		vector_b = p2 - self
		prod = vector_a.cross(vector_b)
		if prod > 0:
			return self.COUNTER_CLOCKWISE
		elif prod < 0:
			return self.CLOCKWISE
		elif vector_a.dot(vector_b) < 0:
			return self.ONLINE_BACK
		elif vector_a.norm() < vector_b.norm():
			return self.ONLINE_FRONT
		else:
			return self.ON_SEGMENT

class Segment(object):

	def __init__(self, p1: Point, p2: Point):
		self.p1 = p1
		self.p2 = p2

	def is_intersected(self, other: 'Segment'):
		return self.p1.ccw(self.p2, other.p1) * self.p1.ccw(self.p2, other.p2) <= 0 and other.p1.ccw(other.p2, self.p1) * other.p1.ccw(other.p2, self.p2) <= 0
from math import gcd
from collections import Counter
n = int(input())
ans = 0
segments = []
for _ in range(n):
	(ax, ay, bx, by) = map(int, input().split())
	segments.append(Segment(Point(ax, ay), Point(bx, by)))
	g = gcd(abs(ax - bx), abs(ay - by))
	ans += g + 1
triangle = {i * (i + 1) >> 1: i for i in range(n + 1)}
intersection = Counter()
for i in range(n):
	for j in range(i + 1, n):
		if segments[i].is_intersected(segments[j]):
			base = segments[j].p2 - segments[j].p1
			d1 = abs(base.cross(segments[i].p1 - segments[j].p1))
			d2 = abs(base.cross(segments[i].p2 - segments[j].p1))
			p2 = segments[i].p2 - segments[i].p1
			if p2.x * d1 % (d1 + d2) == 0 and p2.y * d1 % (d1 + d2) == 0:
				x = segments[i].p1.x + p2.x * d1 // (d1 + d2)
				y = segments[i].p1.y + p2.y * d1 // (d1 + d2)
				intersection[x, y] += 1
ans -= sum((triangle[v] for v in intersection.values()))
print(ans)


import sys

def getIntList():
	return list(map(int, input().split()))
(N,) = getIntList()
zp = []
for i in range(N):
	(ax, ay, bx, by) = getIntList()
	if ax > bx:
		(ax, bx) = (bx, ax)
		(ay, by) = (by, ay)
	zp.append((ax, ay, bx, by))
res = 0

def gcd(a, b):
	if b == 0:
		return a
	return gcd(b, a % b)
zgcd = []
for i in range(N):
	(ax, ay, bx, by) = zp[i]
	tx = abs(bx - ax)
	ty = abs(by - ay)
	g = gcd(tx, ty)
	res += g + 1
	zgcd.append(g)
for i in range(N):
	ax = zp[i][0]
	dax = (zp[i][2] - ax) // zgcd[i]
	ay = zp[i][1]
	day = (zp[i][3] - ay) // zgcd[i]
	cross = []
	for j in range(i + 1, N):
		bx = zp[j][0]
		dbx = (zp[j][2] - bx) // zgcd[j]
		by = zp[j][1]
		dby = (zp[j][3] - by) // zgcd[j]
		t1 = ax * day - ay * dax - bx * day + by * dax
		t2 = dbx * day - dby * dax
		if t2 == 0:
			continue
		if t1 % t2 != 0:
			continue
		k2 = t1 // t2
		if k2 < 0 or k2 > zgcd[j]:
			continue
		if dax != 0:
			t3 = k2 * dbx + bx - ax
			if t3 % dax != 0:
				continue
			k1 = t3 // dax
		else:
			t3 = k2 * dby + by - ay
			if t3 % day != 0:
				continue
			k1 = t3 // day
		if k1 < 0 or k1 > zgcd[i]:
			continue
		cross.append(k1)
	if not cross:
		continue
	cross.sort()
	d = 1
	for j in range(1, len(cross)):
		if cross[j] != cross[j - 1]:
			d += 1
	res -= d
print(res)


def gcd(a, b):
	return a if b == 0 else gcd(b, a % b)

def cross(x1, y1, x2, y2):
	return x1 * y2 - x2 * y1

def line_intersection(line1, line2):
	(px, py, qx, qy) = line1
	(rx, ry, sx, sy) = line2
	xdiff = (px - qx, rx - sx)
	ydiff = (py - qy, ry - sy)
	div = cross(px - qx, rx - sx, py - qy, ry - sy)
	if div == 0:
		raise Exception('lines do not intersect')
	(pq, rs) = (cross(px, py, qx, qy), cross(rx, ry, sx, sy))
	x = cross(pq, rs, px - qx, rx - sx) / div
	y = cross(pq, rs, py - qy, ry - sy) / div
	return (x, y)

def online(line, x, y):
	(a, b, c, d) = line
	if min(a, c) <= x <= max(a, c) and min(b, d) <= y <= max(b, d):
		return True
	else:
		return False

def CF1036E():
	N = int(input())
	lines = []
	count = 0
	for _ in range(N):
		(x1, y1, x2, y2) = map(int, input().split())
		points = gcd(abs(x1 - x2), abs(y1 - y2)) + 1
		count += points
		lines.append((x1, y1, x2, y2))
	for i in range(N):
		d = set()
		for j in range(i + 1, N):
			px = lines[i][0]
			py = lines[i][1]
			qx = lines[j][0]
			qy = lines[j][1]
			rx = lines[i][2] - lines[i][0]
			ry = lines[i][3] - lines[i][1]
			sx = lines[j][2] - lines[j][0]
			sy = lines[j][3] - lines[j][1]
			rs = cross(rx, ry, sx, sy)
			if rs == 0:
				continue
			(x, y) = line_intersection(lines[i], lines[j])
			if not (x % 1 == 0 and y % 1 == 0):
				continue
			if not (online(lines[i], x, y) and online(lines[j], x, y)):
				continue
			d.add((x, y))
		count = count - len(d)
	return int(count)
res = CF1036E()
print(res)


import sys
from math import gcd

class Segment:

	def __init__(self, x1, y1, x2, y2):
		self.x1 = x1
		self.y1 = y1
		self.x2 = x2
		self.y2 = y2

	def dx(self):
		return self.x2 - self.x1

	def dy(self):
		return self.y2 - self.y1

	def pt_cnt(self):
		return gcd(abs(self.dx()), abs(self.dy())) + 1

	def A(self):
		return self.y1 - self.y2

	def B(self):
		return self.x2 - self.x1

	def C(self):
		return -self.A() * self.x1 - self.B() * self.y1

	def inner_x(self, x):
		l = self.x1
		r = self.x2
		if l > r:
			(l, r) = (r, l)
		return l <= x <= r

	def inner_y(self, y):
		l = self.y1
		r = self.y2
		if l > r:
			(l, r) = (r, l)
		return l <= y <= r

	def inner(self, x, y):
		return self.inner_x(x) and self.inner_y(y)

	def intersect(self, other):
		dx = det(self.C(), self.B(), other.C(), other.B())
		dy = det(self.A(), self.C(), other.A(), other.C())
		d = det(self.A(), self.B(), other.A(), other.B())
		if d == 0:
			return (False, 0, 0)
		if dx % d != 0 or dy % d != 0:
			return (False, 0, 0)
		x = -dx // d
		y = -dy // d
		if not self.inner(x, y) or not other.inner(x, y):
			return (False, 0, 0)
		return (True, x, y)

def det(a, b, c, d):
	return a * d - b * c

def main():
	n = int(input())
	segments = []
	ans = 0
	for _ in range(n):
		(x1, y1, x2, y2) = map(int, input().split())
		segment = Segment(x1, y1, x2, y2)
		ans += segment.pt_cnt()
		pt = set()
		for other in segments:
			(result, x, y) = segment.intersect(other)
			if result:
				pt.add((x, y))
		ans -= len(pt)
		segments.append(segment)
	print(ans)

def __starting_point():
	main()
__starting_point()


from math import gcd

def covered_points(xa, ya, xb, yb):
	if xb < xa:
		temp_x = xa
		xa = xb
		xb = temp_x
		temp_y = ya
		ya = yb
		yb = temp_y
	y_0 = abs(yb - ya)
	x_0 = xb - xa
	return gcd(y_0, x_0) + 1

def intersection2(xa, ya, xb, yb, xc, yc, xd, yd):
	if max(xa, xb) < min(xc, xd):
		return ()
	if max(ya, yb) < min(yc, yd):
		return ()
	if xa == xb and xc == xd or (ya == yb and yc == yd):
		return ()
	if ya == yb and yc == yd:
		return ()
	a1 = yb - ya
	b1 = xb - xa
	c1 = xa * a1 - ya * b1
	a2 = yd - yc
	b2 = xd - xc
	c2 = xc * a2 - yc * b2
	det = a1 * b2 - a2 * b1
	if det == 0:
		return ()
	detx = c1 * b2 - c2 * b1
	dety = -a1 * c2 + a2 * c1
	if detx % det != 0 or dety % det != 0:
		return ()
	x = int(detx / det)
	y = int(dety / det)
	if x < min(xa, xb) or x > max(xa, xb) or x < min(xc, xd) or (x > max(xc, xd)):
		return ()
	if y < min(ya, yb) or y > max(ya, yb) or y < min(yc, yd) or (y > max(yc, yd)):
		return ()
	return (x, y)

def __starting_point():
	n = int(input())
	Ax = [0] * n
	Ay = [0] * n
	Bx = [0] * n
	By = [0] * n
	n_cov = 0
	intersections = {}
	for i in range(n):
		line = input().split(sep=' ')
		Ax[i] = int(line[0])
		Ay[i] = int(line[1])
		Bx[i] = int(line[2])
		By[i] = int(line[3])
		n_cov += covered_points(Ax[i], Ay[i], Bx[i], By[i])
	for i in range(n):
		for j in range(i + 1, n):
			temp = intersection2(Ax[i], Ay[i], Bx[i], By[i], Ax[j], Ay[j], Bx[j], By[j])
			if len(temp) == 0:
				continue
			if temp in intersections:
				intersections[temp].append(i)
				intersections[temp].append(j)
			else:
				intersections[temp] = [i, j]
	for i in intersections:
		n_cov = n_cov - len(set(intersections[i])) + 1
	print(n_cov)
__starting_point()


import math

def intersectSeg(pt1, pt2, ptA, ptB, interge=False):

	def is_interge(x):
		if abs(round(x) - x) < 1e-07:
			return True
		return False
	DET_TOLERANCE = 1e-08
	(x1, y1) = pt1
	(x2, y2) = pt2
	dx1 = x2 - x1
	dy1 = y2 - y1
	(x, y) = ptA
	(xB, yB) = ptB
	dx = xB - x
	dy = yB - y
	DET = -dx1 * dy + dy1 * dx
	if math.fabs(DET) < DET_TOLERANCE:
		return (0, 0, 0, 0, 0)
	DETinv = 1.0 / DET
	r = DETinv * (-dy * (x - x1) + dx * (y - y1))
	s = DETinv * (-dy1 * (x - x1) + dx1 * (y - y1))
	xi = (x1 + r * dx1 + x + s * dx) / 2.0
	yi = (y1 + r * dy1 + y + s * dy) / 2.0
	if interge == True and (is_interge(xi) == False or is_interge(yi) == False):
		return (0, 0, 0, 0, 0)
	(xi, yi) = (round(xi), round(yi))
	valid = 1
	if xi < max(min(x1, x2), min(x, xB)) or xi > min(max(x1, x2), max(x, xB)) or yi < max(min(y1, y2), min(y, yB)) or (yi > min(max(y1, y2), max(y, yB))):
		return (0, 0, 0, 0, 0)
	return (xi, yi, valid, r, s)

def gcd(a, b):
	if b == 0:
		return a
	return gcd(b, a % b)
n = int(input())
seg = [list(map(int, input().split())) for _ in range(n)]
cnt = 0
for (x1, y1, x2, y2) in seg:
	cnt += gcd(abs(x1 - x2), abs(y1 - y2)) + 1
sub = 0
for i in range(n):
	d = set()
	for j in range(i + 1):
		(x, y, valid, _, _) = intersectSeg(seg[i][:2], seg[i][2:], seg[j][:2], seg[j][2:], interge=True)
		if valid == 0:
			continue
		if (x, y) not in d:
			d.add((x, y))
	sub += len(d)
print(cnt - sub)


from math import gcd
from bisect import *

class Point:

	def __init__(self, x, y):
		(self.x, self.y) = (x, y)

	def __repr__(self):
		return 'Point' + repr((self.x, self.y))

	def __add__(self, val):
		return Point(self.x + val.x, self.y + val.y)

	def __sub__(self, val):
		return Point(self.x - val.x, self.y - val.y)

	def __mul__(self, ratio):
		return Point(self.x * ratio, self.y * ratio)

	@staticmethod
	def det(A, B):
		return A.x * B.y - A.y * B.x

	@staticmethod
	def dot(A, B):
		return A.x * B.x + A.y * B.y

	def onSegment(self, A, B):
		if self.det(B - A, self - A) != 0:
			return False
		if self.dot(B - A, self - A) < 0:
			return False
		if self.dot(A - B, self - B) < 0:
			return False
		return True

def intPoints(x1, y1, x2, y2):
	(dx, dy) = (abs(x2 - x1), abs(y2 - y1))
	return gcd(dx, dy) + 1

def crosspoint(L1, L2):
	(A, B) = (Point(L1[0], L1[1]), Point(L1[2], L1[3]))
	(C, D) = (Point(L2[0], L2[1]), Point(L2[2], L2[3]))
	(S1, S2) = (Point.det(C - A, D - A), Point.det(C - D, B - A))
	delta = (B - A) * S1
	if S2 == 0 or delta.x % S2 != 0 or delta.y % S2 != 0:
		return None
	delta.x = delta.x // S2
	delta.y = delta.y // S2
	P = A + delta
	if not P.onSegment(A, B) or not P.onSegment(C, D):
		return None
	return (P.x, P.y)
n = int(input())
lines = [tuple((int(z) for z in input().split())) for i in range(n)]
count = dict()
for i in range(n):
	for j in range(i):
		P = crosspoint(lines[i], lines[j])
		if P == None:
			continue
		if not P in count:
			count[P] = 1
		else:
			count[P] += 1
answer = sum((intPoints(*L) for L in lines))
tri = [x * (x + 1) // 2 for x in range(n + 1)]
for z in count:
	k = bisect_right(tri, count[z])
	answer -= k - 1
print(answer)