>>103522957
Field too long, so here's my implementation. Imported the functions from sympy for conciseness, but they're pretty trivial to write yourself if you want to.
from collections import defaultdict
from functools import reduce
from itertools import combinations
from operator import mul
import re
from sympy import mod_inverse
from sympy.ntheory.modular import solve_congruence
W = 101
H = 103
N = 100
quadrants = [0, 0, 0, 0]
points = []
for l in open("input.txt").read().splitlines():
p_x, p_y, v_x, v_y = (
int(n) for n in re.search(r"(-?\d+),(-?\d+) v=(-?\d+),(-?\d+)", l).groups()
)
points.append([p_x, p_y, v_x, v_y])
p_x, p_y = (p_x + v_x * N) % W, (p_y + v_y * N) % H
quadrants[0] += p_x < W // 2 and p_y < H // 2
quadrants[1] += p_x < W // 2 and p_y > H // 2
quadrants[2] += p_x > W // 2 and p_y < H // 2
quadrants[3] += p_x > W // 2 and p_y > H // 2
print(reduce(mul, quadrants))
x_intersection_times, y_intersection_times = defaultdict(int), defaultdict(int)
for (pxa, pya, vxa, vya), (pxb, pyb, vxb, vyb) in combinations(points, 2):
if (vxa - vxb) % W != 0:
x_intersection_times[((pxb - pxa) * mod_inverse(vxa - vxb, W)) % W] += 1
if (vya - vyb) % H != 0:
y_intersection_times[((pyb - pya) * mod_inverse(vya - vyb, H)) % H] += 1
best_time_x = max(x_intersection_times.items(), key=lambda i: i[1])[0]
best_time_y = max(y_intersection_times.items(), key=lambda i: i[1])[0]
print(solve_congruence((best_time_x, W), (best_time_y, H))[0])
