This is a harder version of the problem. In this version, $n \le 50\,000$.

There are $n$ distinct points in three-dimensional space numbered from $1$ to $n$. The $i$-th point has coordinates $(x_i, y_i, z_i)$. The number of points $n$ is even.

You'd like to remove all $n$ points using a sequence of $\frac{n}{2}$ snaps. In one snap, you can remove any two points $a$ and $b$ that have not been removed yet and form a perfectly balanced pair. A pair of points $a$ and $b$ is perfectly balanced if no other point $c$ (that has not been removed yet) lies within the axis-aligned minimum bounding box of points $a$ and $b$.

Formally, point $c$ lies within the axis-aligned minimum bounding box of points $a$ and $b$ if and only if $\min(x_a, x_b) \le x_c \le \max(x_a, x_b)$, $\min(y_a, y_b) \le y_c \le \max(y_a, y_b)$, and $\min(z_a, z_b) \le z_c \le \max(z_a, z_b)$. Note that the bounding box might be degenerate. 

Find a way to remove all points in $\frac{n}{2}$ snaps.


-----Input-----

The first line contains a single integer $n$ ($2 \le n \le 50\,000$; $n$ is even), denoting the number of points.

Each of the next $n$ lines contains three integers $x_i$, $y_i$, $z_i$ ($-10^8 \le x_i, y_i, z_i \le 10^8$), denoting the coordinates of the $i$-th point.

No two points coincide.


-----Output-----

Output $\frac{n}{2}$ pairs of integers $a_i, b_i$ ($1 \le a_i, b_i \le n$), denoting the indices of points removed on snap $i$. Every integer between $1$ and $n$, inclusive, must appear in your output exactly once.

We can show that it is always possible to remove all points. If there are many solutions, output any of them.


-----Examples-----
Input
6
3 1 0
0 3 0
2 2 0
1 0 0
1 3 0
0 1 0

Output
3 6
5 1
2 4

Input
8
0 1 1
1 0 1
1 1 0
1 1 1
2 2 2
3 2 2
2 3 2
2 2 3

Output
4 5
1 6
2 7
3 8



-----Note-----

In the first example, here is what points and their corresponding bounding boxes look like (drawn in two dimensions for simplicity, as all points lie on $z = 0$ plane). Note that order of removing matters: for example, points $5$ and $1$ don't form a perfectly balanced pair initially, but they do after point $3$ is removed.  [Image]

Difficulty: interview

Starter code:

Solutions:
n = int(input())
a = [tuple(map(int, input().split())) for _ in range(n)]
s = {a[q]: q+1 for q in range(n)}
a.sort()
q, q1 = 0, 1
d, d1, d2 = [[[a[0]]]], [], []
while q1 < n:
    while q1 < n and a[q][0] == a[q1][0]:
        while q1 < n and a[q][1] == a[q1][1]:
            d[-1][-1].append(a[q1])
            q1 += 1
        if q1 < n and a[q][0] == a[q1][0]:
            d[-1].append([a[q1]])
            q = q1
            q1 += 1
    if q1 < n:
        d.append([[a[q1]]])
        q = q1
        q1 += 1
for q in range(len(d)):
    for q1 in range(len(d[q])):
        for q2 in range(1, len(d[q][q1]), 2):
            print(s[d[q][q1][q2-1]], s[d[q][q1][q2]])
        if len(d[q][q1]) % 2 == 1:
            d[q][q1] = d[q][q1][-1]
        else:
            d[q][q1] = -1
for q in range(len(d)):
    d1.append([])
    for q1 in range(len(d[q])):
        if d[q][q1] != -1:
            d1[-1].append(d[q][q1])
for q in range(len(d1)):
    for q1 in range(1, len(d1[q]), 2):
        print(s[d1[q][q1-1]], s[d1[q][q1]])
    if len(d1[q]) % 2 == 1:
        d2.append(d1[q][-1])
for q in range(1, len(d2), 2):
    print(s[d2[q-1]], s[d2[q]])


import sys

n = -1
data = []
for line in sys.stdin:
    if n == -1:
        n = int(line)
        continue
    data.append(tuple([int(i) for i in line.split()]))
    if len(data) == n:
        break

data = list(zip(data, list(range(1, n+1))))
data.sort()

data2 = []
removedLast = False
for i in range(n):
    if removedLast:
        removedLast = False
        continue
    if i+1<n and data[i][0][0] == data[i+1][0][0] and data[i][0][1] == data[i+1][0][1]:
        print("%d %d"%(data[i][1], data[i+1][1]))
        removedLast = True
    else:
        data2.append(data[i])

data3 = []
removedLast = False
for i in range(len(data2)):
    if removedLast:
        removedLast = False
        continue
    if i+1<len(data2) and data2[i][0][0] == data2[i+1][0][0]:
        print("%d %d"%(data2[i][1], data2[i+1][1]))
        removedLast = True
    else:
        data3.append(data2[i])

for i in range(0, len(data3), 2):
    print("%d %d"%(data3[i][1], data3[i+1][1]))
    


from sys import stdin,stdout
n=int(stdin.readline())
l=[]
for i in range(n):
    x=list(map(int,stdin.readline().split()))
    l.append([x,i+1])
l.sort()
#print(l)
deleted=[0]*n
for i in range(n-1):
    if deleted[i]:
        continue
    if l[i][0][0]==l[i+1][0][0] and l[i][0][1]==l[i+1][0][1]:
        print(l[i][1],l[i+1][1])
        deleted[i]=deleted[i+1]=1
remaining=[]
for i in range(n):
    if deleted[i]==0:
        remaining.append(l[i])
del deleted
n=len(remaining)
deleted=[0]*n
for i in range(n-1):
    if deleted[i]:
        continue
    if remaining[i][0][0]==remaining[i+1][0][0]:
        print(remaining[i][1],remaining[i+1][1])
        deleted[i]=deleted[i+1]=1
last=[]
for i in range(n):
    if deleted[i]==0:
        last.append(remaining[i])
for i in range(0,len(last),2):
    print(last[i][1],last[i+1][1])


import sys
from collections import defaultdict
readline = sys.stdin.readline
        
N = int(readline())

XYZ = [tuple(map(int, readline().split())) for _ in range(N) ]
XYZI = [(x, y, z, i) for i, (x, y, z) in enumerate(XYZ, 1)]


XYZI.sort()
X, Y, Z, IDX = map(list, zip(*XYZI))

Di = defaultdict(list)
for i in range(N):
    x, y, z, idx = X[i], Y[i], Z[i], IDX[i]
    Di[x].append((y, z, idx))
Ans = []
Ama = []

for L in Di.values():
    D2 = defaultdict(int)
    for y, _, _ in L:
        D2[y] += 1
    pre = None
    st = None
    for y, z, i in L:
        if st is None and D2[y] == 1:
            if pre is not None:
                Ans.append((pre, i))
                pre = None
            else:
                pre = i
        else:
            if st is not None:
                Ans.append((st, i))
                st = None
            else:
                st = i
        D2[y] -= 1
    
    if st:
        Ama.append(st)
    if pre:
        Ama.append(pre)
Le = len(Ama)
for i in range(Le//2):
    Ans.append((Ama[2*i], Ama[2*i+1]))

for a in Ans:
    sys.stdout.write('{} {}\n'.format(*a))

n = int(input())
coors =[]
for i in range(n):
    coors.append(list(map(int,input().split()))+[i+1])
# print(coors)
coors=sorted(coors,key=lambda x:(x[0],x[1],x[2]))
# for i in range(0,n,2):
#     if coors[i][0]==coors[i+1][0]:
#         if coors[i][1]==coors[i+1][1]:
#             print(coors[i][3],coors[i+1][3])
#         else:
left=[]
def oneD(co,remain):
    if len(co)%2==0:
        for i in range(0, len(co), 2):
            print(co[i][3],co[i+1][3])
    else:
        remain.append(co[-1])
        for i in range(0, len(co)-1, 2):
            print(co[i][3],co[i+1][3])
    return remain
def twoD(co):
    if len(co)==1:
        left.append(*co)
        return
    remain = []
    l = 1
    for i in range(len(co)-1):
        if co[i][1] == co[i + 1][1]:
            l += 1
        else:
            remain = oneD(co[i+1 -l:i + 1],remain)
            l=1
    remain = oneD(co[len(co)-l:],remain)
    if len(remain)%2==0:
        for i in range(0, len(remain), 2):
            print(remain[i][3],remain[i+1][3])
    else:
        left.append(remain[-1])
        for i in range(0, len(remain)-1, 2):
            print(remain[i][3],remain[i+1][3])

l=1
for i in range(n-1):
    if coors[i][0]==coors[i+1][0]:
        l+=1
    else:
        twoD(coors[i+1-l:i+1])
        l=1
twoD(coors[n-l:])
# print(left)
for i in range(0,len(left),2):
    print(left[i][3],left[i+1][3])

import sys
import math
input = sys.stdin.readline

n=int(input())

arr=[]
for i in range(n):
	arr.append(list(map(int,input().split())))
	arr[-1].append(i+1)

ans=[]
arr.sort()
#print(arr)
done=[0]*n
i=0
j=1
while i<n and j<n:
	if arr[i][0]==arr[j][0] and arr[i][1]==arr[j][1]:
		ans.append([arr[i][3], arr[j][3]])
		done[arr[i][3]-1]+=1
		done[arr[j][3]-1]+=1
		i+=2
		j+=2
	else:
		i+=1
		j+=1
#print(ans)
#print(done)
arr1=[]
for i in range(n):
	if done[arr[i][3]-1]==0:
		arr1.append(arr[i])
n=len(arr1)
#print("arr1",arr1)
i=0
j=1
while i<n and j<n:
	if arr1[i][0]==arr1[j][0]:
		ans.append([arr1[i][3], arr1[j][3]])
		done[arr1[i][3]-1]+=1
		done[arr1[j][3]-1]+=1
		i+=2
		j+=2
	else:
		i+=1
		j+=1
#print(ans)
arr2=[]
#print("arr1", arr1)
for i in range(n):
	if done[arr1[i][3]-1]==0:
		arr2.append(arr1[i])
#print(done)
n=len(arr2)
#print(arr2)
for i in range(n):
	if i%2==0:
		ans.append([arr2[i][3], arr2[i+1][3]])

for i in range(len(ans)):
	print(*ans[i])

def solve(n, p):
    p = sorted(p, key=lambda x: (x[0], x[1], x[2]))
    # print('solve', p)
    p.append([10**20, 10**20, 10**20, 10**20])
    # eliminate matching x coord ones:
    c = 1
    pn = []
    for i in range(1, len(p)):
        if p[i][0] == p[i-1][0]:
            c += 1
        else:
            if c >= 2:
                s = solve_2d(p[i-c:i])
                if s:
                    pn.append(s)
            else:
                pn.append(p[i-1])
            c = 1
    # print(pn)
    for i in range(0, len(pn)-1, 2):
        print(pn[i][3], pn[i+1][3])


def solve_2d(p1):
    # print('solve_2d', p1)
    p1.append([10**20, 10**20, 10**20, 10**20])
    c = 1
    p1n = []
    i = 1
    while i < len(p1):
        if p1[i][1] == p1[i-1][1]:
            print(p1[i][3], p1[i-1][3])
            i += 2
        else:
            p1n.append(p1[i-1])
            i += 1
    for i in range(0, len(p1n)-1, 2):
        print(p1n[i][3], p1n[i+1][3])
    if len(p1n) % 2 == 1:
        return p1n[-1]
    else:
        return 0


def main():
    n = int(input())
    p = []
    for i in range(n):
        pt = [int(i) for i in input().split()]
        p.append(pt + [i+1])
    solve(n, p)


main()


n = int(input())
coord = []
for i in range(n):
    coord.append(tuple(list(map(int,input().split()))+[i+1]))
coord.sort()
def odind(l):
    l.sort()
    for i in range(len(l)//2):
        print(l[2*i][-1],l[2*i+1][-1])
    if len(l)%2==0:
        return []
    return l[-1]
def dvad(l):
    cx = {}
    ans1 = []
    for x,y,n in l:
        if x in cx:
            cx[x].append([y,n])
        else:
            cx[x] = [[y,n]]
    for x,j in cx.items():
        v = odind(j)
        if len(v)!=0:
            ans1.append(tuple([x]+v))
    ans1.sort()
    for i in range(len(ans1)//2):
        print(ans1[2*i][-1],ans1[2*i+1][-1])
    if len(ans1)%2==0:
        return []
    return ans1[-1]
def trid(l):
    cx = {}
    ans1 = []
    for x, y,z, n in l:
        if x in cx:
            cx[x].append((y, z, n))
        else:
            cx[x] = [(y, z, n)]
    for x,j in cx.items():
        v = list(dvad(j))
        if len(v) != 0:
            ans1.append(tuple([x]+v))
    ans1.sort()
    for i in range(len(ans1) // 2):
        print(ans1[2 * i][-1], ans1[2 * i + 1][-1])
trid(coord)

from sys import stdin
input = stdin.readline
n = int(input())
l = [list(map(int,input().split())) for i in range(n)]
a = [[l[i],i+1] for i in range(n)]
used = [0] * (n+1)
a.sort()
d = {}
for i in a:
	d[i[0][0]] = []
for i in a:
	d[i[0][0]].append([i[0][1],i[0][2], i[1]])
for x in d:
	#doing for 2d
	l = d[x] # this is restricted to x-coordinate
	b = {}
	for i in l:
		b[i[0]] = []
	for i in l:
		b[i[0]].append([i[1], i[2]]) #{z, n}
	for y in b:
		pom = 0
		while pom + 1 < len(b[y]):
			print(b[y][pom][1], b[y][pom+1][1])
			used[b[y][pom][1]] = 1
			used[b[y][pom+1][1]] = 1
			pom += 2
	c = []
	for gg in l:
		if used[gg[2]] == 0:
			c.append(gg)
	cc = [[c[gg][0], c[gg][2]] for gg in range(len(c))]
	g = 0
	while g < len(cc) - 1:
		print(cc[g][-1], cc[g+1][-1])
		used[cc[g][-1]] = 1
		used[cc[g+1][-1]] = 1
		g += 2
c = []
for i in range(0, n):
	if used[a[i][1]] == 0:
		c.append([a[i][0][0],a[i][1]])
c.sort()
i = 0
while i < len(c) - 1:
	print(c[i][1], c[i+1][1])
	i += 2

n = int(input())
points = []

for i in range(n):
    points.append(list(map(int, input().split())))
    points[i].append(i+1)
points.sort()

rem = []

i = 0
while i < len(points):
    if i < len(points)-1 and points[i][:2] == points[i+1][:2]:
        print(points[i][3],points[i+1][3])
        i+=2
    else:
        rem.append(points[i])
        i+=1
points = list(rem)
rem = []
i = 0
while i < len(points):
    if i < len(points)-1 and points[i][0] == points[i+1][0]:
        print(points[i][3],points[i+1][3])
        i+=2
    else:
        rem.append(points[i])
        i+=1
#print(rem)
points = list(rem)
rem = []
i = 0
while i < len(points)-1:
    print(points[i][3], points[i+1][3])
    i+=2





"""
    > File Name: c.py
    > Author: Code_Bear
    > Mail: secret 
    > Created Time: Thu Oct 17 16:34:03 2019
"""
from collections import OrderedDict
def SortPoint(p, ids, k, D):
    if k == D:
        return ids[0]
    maps = OrderedDict()
    for i in ids:
        if p[i][k] not in maps:
            maps[p[i][k]] = []
        maps[p[i][k]].append(i)
    a = []
    for i in sorted(maps.keys()):
        cnt = SortPoint(p, maps[i], k + 1, D)
        if cnt != -1: a.append(cnt)
    for i in range(0, len(a), 2):
        if i + 1 < len(a):
            print(a[i] + 1, a[i + 1] + 1)
    return -1 if len(a) % 2 == 0 else a[-1]
            

def solver():
    n = int(input())
    p = []
    ids = [i for i in range(n)]
    for i in range(n):
        point = list(map(int, input().split()))
        p.append(point)
    SortPoint(p, ids, 0, 3)
    
def main():
    while 1:
        try:
            solver()
        except EOFError:
            break

def __starting_point():
    main()
__starting_point()

def solve(points,ans,coords,x):
    arr = []
    curr = []
    y = points[0][0]
    for p in points:
        if p[0] == y:
            curr.append((y,p[1],p[2]))
        else:
            arr.append(curr)
            y = p[0]
            curr = [(y,p[1],p[2])]

    arr.append(curr)

    arr1 = []

    for i in arr:
        while len(i) >= 2:
            ans.append((i.pop(0)[2],i.pop(0)[2]))

        if len(i) == 1:
            arr1.append((i[0][0],i[0][1],i[0][2]))

    while len(arr1) >= 2:
        ans.append((arr1.pop(0)[2],arr1.pop(0)[2]))

    coords[x] = arr1[:]

def main():
    n = int(input())
    points = []
    coords = {}
    for i in range(n):
        x,y,z = list(map(int,input().split()))
        if x not in list(coords.keys()):
            coords[x] = [(y,z,i+1)]
        else:
            coords[x].append((y,z,i+1))

    ans = []

    for i in list(coords.keys()):
        coords[i].sort()
        if len(coords[i]) > 1:
            solve(coords[i],ans,coords,i)

        if len(coords[i]) == 1:
            points.append((i,coords[i][0][0],coords[i][0][1],coords[i][0][2]))

    points.sort()

    for i in range(0,len(points),2):
        a,b = points[i][3],points[i+1][3]
        ans.append((a,b))

    for i in ans:
        print(i[0],i[1])

main()


def d(le):
    for i in range(0,len(le),2):
    #    if le[i][3]==29 and le[i+1][3]==38:
     #       print(me)
      #      return ''
        print(le[i][3],le[i+1][3])
    
def dd(it):
    nonlocal do
   # print(it)
    tt={}
    ddd=[]
    it.sort(key=lambda x:x[1])
    l=list(set([i[1] for i in it]))
    l.sort()
    for i in l:
        tt[i]=[]
    for i in it:
        tt[i[1]].append(i)
    for i in l:
        x=tt[i]
        if len(x)%2==1:
            ddd.append(x[-1][:])
            x=x[:-1]
        if len(x)>0:
            x.sort(key=lambda x:x[2])
            d(x[:])
    if len(ddd)%2==1:
        do.append(ddd[-1][:])
        ddd=ddd[:-1]
    if len(ddd)>0:
        ddd.sort(key=lambda x:x[1])
        d(ddd)
n=int(input())
it=[]
for i in range(n):
    it.append(list(map(int,input().split()))+[i+1])
ss={}
it.sort()
l=list(set([i[0] for i in it]))
l.sort()
for i in l:
    ss[i]=[]
for i in it:
    ss[i[0]].append(i)
do=[]
for i in l:
    x=ss[i]
    dd(x[:])
#do.sort(key=lambda x:x[2])
d(do)

    
    


n = int(input())
sentinel = [(10**9, 10**9, 10**9)]
points = sorted(tuple(map(int, input().split())) + (i,) for i in range(1, n+1)) + sentinel
ans = []
rem = []

i = 0
while i < n:
    if points[i][:2] == points[i+1][:2]:
        ans.append(f'{points[i][-1]} {points[i+1][-1]}')
        i += 2
    else:
        rem.append(i)
        i += 1

rem += [n]
rem2 = []

n = len(rem)
i = 0
while i < n-1:
    if points[rem[i]][0] == points[rem[i+1]][0]:
        ans.append(f'{points[rem[i]][-1]} {points[rem[i+1]][-1]}')
        i += 2
    else:
        rem2.append(points[rem[i]][-1])
        i += 1

print(*ans, sep='\n')
for i in range(0, len(rem2), 2):
    print(rem2[i], rem2[i+1])