Geolocation (Codeforces 1220G)
Tags: Others, geometry

You are working for the Gryzzl company, headquartered in Pawnee, Indiana.The new national park has been opened near Pawnee recently and you are to implement a geolocation system, so people won't get lost. The concept you developed is innovative and minimalistic. There will be $$$n$$$ antennas located somewhere in the park. When someone would like to know their current location, their Gryzzl hologram phone will communicate with antennas and obtain distances from a user's current location to all antennas.Knowing those distances and antennas locations it should be easy to recover a user's location... Right? Well, almost. The only issue is that there is no way to distinguish antennas, so you don't know, which distance corresponds to each antenna. Your task is to find a user's location given as little as all antennas location and an unordered multiset of distances.

Input: The first line of input contains a single integer $$$n$$$ ($$$2 \leq n \leq 10^5$$$) which is the number of antennas.The following $$$n$$$ lines contain coordinates of antennas, $$$i$$$-th line contain two integers $$$x_i$$$ and $$$y_i$$$ ($$$0 \leq x_i,y_i \leq 10^8$$$). It is guaranteed that no two antennas coincide.The next line of input contains integer $$$m$$$ ($$$1 \leq n \cdot m \leq 10^5$$$), which is the number of queries to determine the location of the user.Following $$$m$$$ lines contain $$$n$$$ integers $$$0 \leq d_1 \leq d_2 \leq \dots \leq d_n \leq 2 \cdot 10^{16}$$$ each. These integers form a multiset of squared distances from unknown user's location $$$(x;y)$$$ to antennas.For all test cases except the examples it is guaranteed that all user's locations $$$(x;y)$$$ were chosen uniformly at random, independently from each other among all possible integer locations having $$$0 \leq x, y \leq 10^8$$$.

Output: For each query output $$$k$$$, the number of possible a user's locations matching the given input and then output the list of these locations in lexicographic order.It is guaranteed that the sum of all $$$k$$$ over all points does not exceed $$$10^6$$$.

Note: As you see in the second example, although initially a user's location is picked to have non-negative coordinates, you have to output all possible integer locations.

Sample Inputs:
3
0 0
0 1
1 0
1
1 1 2

4
0 0
0 1
1 0
1 1
2
0 1 1 2
2 5 5 8

Sample Outputs:
1 1 1

4 0 0 0 1 1 0 1 1 
4 -1 -1 -1 2 2 -1 2 2

Solution (Python):
from sys import stdin, stdout
from math import isqrt
n = int(input())
(x0, y0) = (0, 0)
(sx, sy, s2) = (0, 0, 0)
qx = []
qy = []
for i in range(n):
    (px, py) = map(int, input().split())
    qx.append(px)
    qy.append(py)
    sx += px
    sy += py
    s2 += px * px + py * py
(p, q) = (None, None)
for i in range(n):
    if (n * qx[i], n * qy[i]) != (sx, sy):
        (p, q) = (qx[i], qy[i])
        break
m = int(input())
a = 2 * p * n - 2 * sx
b = 2 * q * n - 2 * sy
pc = n * (p * p + q * q) - s2
pdt = -(a * a * p * p + 2 * a * b * p * q + b * b * q * q)
sr = a * a + b * b

def gethash(x):
    return x * 13331 + 23333
pans = []
for _ in range(m):
    d = list(map(int, input().split()))
    sd = sum(d)
    hd = 0
    for j in range(n):
        hd ^= gethash(d[j])
    ans = []
    for j in range(n):
        if j >= 1 and d[j] == d[j - 1]:
            continue
        c = pc + sd - n * d[j]
        dt = pdt + a * a * d[j] + 2 * a * c * p + b * b * d[j] + 2 * b * c * q - c * c
        if dt < 0:
            continue
        idt = isqrt(dt)
        if idt * idt != dt:
            continue

        def check(idt):
            x = b * idt - a * b * q + a * c + b * b * p
            y = -a * idt + a * a * q + b * c - a * b * p
            if x % sr != 0:
                return
            (x, y) = (x // sr, y // sr)
            z = 0
            for k in range(n):
                z ^= gethash((x - qx[k]) ** 2 + (y - qy[k]) ** 2)
            if z == hd:
                ans.append([x, y])
        check(idt)
        if idt != 0:
            check(-idt)
    ans.sort()
    buf = [len(ans)]
    for x in ans:
        buf.extend(x)
    pans.append(' '.join(map(str, buf)))
stdout.write('\n'.join(pans))

Solution (Python):
import sys
import math
n = int(input())
x = [0] * n
y = [0] * n
for i in range(n):
    (x[i], y[i]) = map(int, input().split())
sx = sum(x)
sy = sum(y)
for i in range(n):
    x[i] = n * x[i] - sx
    y[i] = n * y[i] - sy
m = int(input())
d = [0] * n
e = [0] * n
HD = 0

def check(a, b):
    global HD
    HE = 0
    for i in range(n):
        HE ^= hash((a - x[i]) * (a - x[i]) + (b - y[i]) * (b - y[i]))
    return HD == HE

def sqrt(x):
    nn = int(x)
    if nn == 0:
        return 0
    (fa, fb) = divmod(nn.bit_length(), 2)
    x = 2 ** (fa + fb)
    while True:
        y = (x + nn // x) // 2
        if y >= x:
            return x
        x = y

def hash(x):
    return x * 9991 + 43
pans = []

def solve():
    global d
    d = list(map(int, input().split()))
    c = 0
    d = [p * n * n for p in d]
    for i in range(n):
        c += d[i] - x[i] * x[i] - y[i] * y[i]
    assert c % n == 0
    c //= n
    ans = []
    ax = x[0]
    ay = y[0]
    if ax is 0 and ay is 0:
        ax = x[1]
        ay = y[1]
    rev = 0
    if ay == 0:
        ay = ax
        ax = 0
        rev = 1
    d.sort()
    global HD
    HD = 0
    for p in d:
        HD ^= hash(p)
    old = -1
    for p in d:
        if p == old:
            continue
        old = p
        a = c + ax * ax + ay * ay - p
        if a % 2 != 0:
            continue
        a //= 2
        A = ax * ax + ay * ay
        B = a * ax
        C = a * a - ay * ay * c
        D = B * B - A * C
        if D < 0:
            continue
        sD = sqrt(D)
        if D != sD * sD:
            continue
        if (B + sD) % A == 0:
            qx = (B + sD) // A
            qy = (a - ax * qx) // ay
            if rev:
                t = qx
                qx = qy
                qy = t
            if (qx + sx) % n == 0 and (qy + sy) % n == 0 and check(qx, qy):
                qx = (qx + sx) // n
                qy = (qy + sy) // n
                ans.append([qx, qy])
        if sD == 0:
            continue
        if (B - sD) % A == 0:
            qx = (B - sD) // A
            qy = (a - ax * qx) // ay
            if rev:
                t = qx
                qx = qy
                qy = t
            if (qx + sx) % n == 0 and (qy + sy) % n == 0 and check(qx, qy):
                qx = (qx + sx) // n
                qy = (qy + sy) // n
                ans.append([qx, qy])
    ans.sort()
    buf = []
    buf.append(len(ans))
    for p in ans:
        buf.append(p[0])
        buf.append(p[1])
    global pans
    pans.append(' '.join(map(str, buf)))
while m > 0:
    m -= 1
    solve()
sys.stdout.write('\n'.join(pans))

Solution (Python):
from sys import stdin, stdout
from math import isqrt
n = int(input())
(x0, y0) = (0, 0)
(sx, sy, s2) = (0, 0, 0)
qx = []
qy = []
for i in range(n):
    (px, py) = map(int, input().split())
    qx.append(px)
    qy.append(py)
    sx += px
    sy += py
    s2 += px * px + py * py
(p, q) = (None, None)
for i in range(n):
    if (n * qx[i], n * qy[i]) != (sx, sy):
        (p, q) = (qx[i], qy[i])
        break
m = int(input())
a = 2 * p * n - 2 * sx
b = 2 * q * n - 2 * sy
pc = n * (p * p + q * q) - s2
pdt = -(a * a * p * p + 2 * a * b * p * q + b * b * q * q)
sr = a * a + b * b

def gethash(x):
    return x * 13331 + 23333
pans = []
for _ in range(m):
    d = list(map(int, input().split()))
    sd = sum(d)
    hd = 0
    for j in range(n):
        hd ^= gethash(d[j])
    ans = []
    for j in range(n):
        if j >= 1 and d[j] == d[j - 1]:
            continue
        c = pc + sd - n * d[j]
        dt = pdt + a * a * d[j] + 2 * a * c * p + b * b * d[j] + 2 * b * c * q - c * c
        if dt < 0:
            continue
        idt = isqrt(dt)
        if idt * idt != dt:
            continue

        def check(idt):
            x = b * idt - a * b * q + a * c + b * b * p
            y = -a * idt + a * a * q + b * c - a * b * p
            if x % sr != 0:
                return
            (x, y) = (x // sr, y // sr)
            z = 0
            for k in range(n):
                z ^= gethash((x - qx[k]) ** 2 + (y - qy[k]) ** 2)
            if z == hd:
                ans.append([x, y])
        check(idt)
        if idt != 0:
            check(-idt)
    ans.sort()
    buf = [len(ans)]
    for x in ans:
        buf.extend(x)
    pans.append(' '.join(map(str, buf)))
stdout.write('\n'.join(pans))