Ma5termind is going through the Euler's book of knowledge when he finds the following function: 

def function(L,R):
    sum=0
    for i in range(L,R+1):              #[L,R]
        for j in range(1,i+1):          #[1,i]
            cnt=0
            if (j)*(j+1) == 2*i :
                for k in range(1,i+1):  #[1,i]
                    if __gcd(k,i)==1:   #__gcd Greatest Common Divisor 
                        cnt+=1
            sum=sum+cnt
    return sum

The function seems weird as well as interesting. Don't you think ?? 

Fascinated by the function, Ma5termind codes it and successfully compiles it in the first go itself. After all, he is a good coder. Ma5termind checks the output computed by the function for all the values of $\mathbf{L}$ and $\mathbf{R}$ in the range $\mathbf{[1,100]}$ and finds that the results are correct.

But as we all know, Ma5termind loves large numbers - very large numbers. When he inputs $\mathbf{L=1}$ & $\mathbf{R=10^{12}}$, the code does not produce any output #TimeLimitExceeded. He still wants to know the answer for such large ranges. Can you help him by writing a code that can also produce an answer for such large ranges?   

Input Format 

First line of input contains a single integer $\mathbf{T}$ denoting the number of test cases. First and the only line of each test case contains two space separated integers $\mathbf{L}$ and $\mathbf{R}$ denoting inputs to the above function.  

Output Format 

Output consists of $\mathbf{T}$ lines each containing result computed by the function when corresponding range is inputted to the function.

Constraints 

$\mathbf{1\leq T\leq10^5}$ 

$\mathbf{1\leq L\leq R\leq10^{12}}$

Sample Input  

2
1 10
2 5

Sample Output  

9
2

Explanation 

Pass Parameters to the function and check the result.

Source: hackerrank | Difficulty: MEDIUM
Tags: Number theory


Solutions:
import sys, random
from math import sqrt

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

def function(L, R):
	sum = 0
	for i in range(L, R + 1):
		for j in range(1, i + 1):
			cnt = 0
			if j * (j + 1) == 2 * i:
				for k in range(1, i + 1):
					if gcd(k, i) == 1:
						cnt += 1
			sum = sum + cnt
	return sum

def eulerSieve(N):
	if N == 1:
		return [0, 1]
	limit = int(N ** 0.5) + 1
	sieve = [0] * (N + 1)
	for i in range(3, limit, 2):
		if sieve[i] == 0:
			for j in range(i * i, N + 1, 2 * i):
				sieve[j] = i
	sieve[0:3] = [0, 1, 1]
	for i in range(3, N + 1):
		p = sieve[i]
		if p != 0:
			j = i // p
			sieve[i] = sieve[j] * (p if j % p == 0 else p - 1)
		elif i & 3 == 0:
			sieve[i] = sieve[i >> 1] << 1
		elif i & 1 == 0:
			sieve[i] = sieve[i >> 1]
		else:
			sieve[i] = i - 1
	return sieve

def solve():

	def f(n):
		return int((-1 + (1 + 8 * n) ** 0.5) / 2.0)
	N = 1420000
	phis = eulerSieve(N + 1)
	phi2 = [0] * (N + 1)
	for i in range(1, N + 1):
		j = i * (i + 1)
		phi2[i] = phis[i] * phis[i + 1]
		if not j & 3:
			phi2[i] //= 2
	for j in range(1, N + 1):
		phi2[j] += phi2[j - 1]
	stdin = sys.stdin
	MOD = 10 ** 9 + 7
	T = int(stdin.readline())
	for _ in range(T):
		(L, R) = map(int, stdin.readline().split())
		L = f(L - 1)
		R = f(R)
		print(phi2[R] - phi2[L])
solve()


from math import sqrt
from sys import stderr
import itertools as itool
MAX = 10 ** 12
MAXTRI = int(sqrt(8 * MAX + 1) - 1) // 2
primefactor = [0, 1] + [None] * (MAXTRI - 1)
for p in range(2, len(primefactor)):
	if primefactor[p] is None:
		primefactor[p] = p
		skip = p if p == 2 else 2 * p
		for n in range(p * p, len(primefactor), skip):
			primefactor[n] = p
mutpri = [0, 1]
while len(mutpri) < len(primefactor):
	n = len(mutpri)
	p = primefactor[n]
	(rest, pcount) = (n // p, 1)
	while rest % p == 0:
		(rest, pcount) = (rest // p, pcount + 1)
	mutpri.append((p - 1) * p ** (pcount - 1) * mutpri[rest])
mutpritri = [mutpri[i // 2] * mutpri[i + 1] if i & 1 == 0 else mutpri[i] * mutpri[(i + 1) // 2] for i in range(MAXTRI)]
csmutpritri = list(itool.accumulate(itool.chain([0], mutpritri)))
ncases = int(input())
for _ in range(ncases):
	(L, R) = [int(f) for f in input().split()]
	start = int(sqrt(8 * L + 1) - 2) // 2
	stop = int(sqrt(8 * R + 1) - 1) // 2
	while L > start * (start + 1) // 2:
		start += 1
	assert R >= stop * (stop + 1) // 2
	print(csmutpritri[stop + 1] - csmutpritri[start])


from math import sqrt

def xmax(R):
	a = int(sqrt(2 * R))
	if a * (a + 1) <= 2 * R:
		return a
	else:
		return a - 1

def xmin(L):
	a = int(sqrt(2 * L))
	if a * (a + 1) >= 2 * L:
		return a
	return a + 1
(MIN_L, MAX_R) = (1, 10 ** 12)
MIN_X = xmin(MIN_L)
MAX_X = xmax(MAX_R)

def primesTo(n):
	potential = [0, 0, 1] + n // 2 * [1, 0]
	if len(potential) > n + 1:
		potential.pop()
	for p in range(3, int(n ** 0.5) + 1, 2):
		if p:
			for pq in range(3 * p, n + 1, 2 * p):
				potential[pq] = 0
	return [p for p in range(n + 1) if potential[p]]
assert primesTo(101) == [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97, 101]

def phisTo(n):
	primes = primesTo(n)
	phis = list(range(n + 1))
	for p in primes:
		for pq in range(p, n + 1, p):
			phis[pq] = phis[pq] // p * (p - 1)
	return phis
assert phisTo(10) == [0, 1, 1, 2, 2, 4, 2, 6, 4, 6, 4]

def phiProdHalvesTo(maxX):
	phis = phisTo(maxX + 1)
	phiPHT = []
	for x in range(maxX + 1):
		if x % 2:
			phiPHT.append(phis[x] * phis[(x + 1) // 2])
		else:
			phiPHT.append(phis[x // 2] * phis[x + 1])
	return phiPHT

def partialSums(things):
	if not things:
		return []
	ps = [things[0]]
	for i in range(1, len(things)):
		ps.append(ps[-1] + things[i])
	return ps
phiPHTsums = partialSums(phiProdHalvesTo(MAX_X))
assert len(phiPHTsums) == MAX_X + 1

def query(L, R):
	return phiPHTsums[xmax(R)] - phiPHTsums[xmin(L) - 1]
num_tests = int(input())
for test in range(num_tests):
	(L, R) = map(int, input().split())
	print(query(L, R))


from math import floor, ceil, sqrt
max_n = 1414219
totients = [1] * max_n
sums = [0] * max_n
for i in range(2, max_n):
	if totients[i] == 1:
		for j in range(i, max_n, i):
			totients[j] *= i - 1
			k = j // i
			while k % i == 0:
				totients[j] *= i
				k /= i
for i in range(1, max_n - 2, 2):
	sums[i] = sums[i - 1] + totients[i] * totients[(i + 1) // 2]
	sums[i + 1] = sums[i] + totients[(i + 1) // 2] * totients[i + 2]
for _ in range(int(input())):
	(l, r) = map(int, input().split())
	l = ceil((sqrt(8 * l + 1) - 1) / 2)
	r = floor((sqrt(8 * r + 1) - 1) / 2)
	print(sums[r] - sums[l - 1])


def isqrt(x):
	if x < 0:
		raise ValueError('negative')
	b = x.bit_length() - 1 & -2
	y = 0
	while b >= 0:
		y <<= 1
		z = y + 1
		if z * z <= x >> b:
			y = z
		b -= 2
	return y
lim = 1500000
phi = [1] * lim
j = 4
while j < lim:
	for k in range(j, lim, j):
		phi[k] <<= 1
	j <<= 1
for i in range(3, lim, 2):
	if phi[i] > 1:
		continue
	for k in range(i, lim, i):
		phi[k] *= i - 1
	j = i * i
	while j < lim:
		for k in range(j, lim, j):
			phi[k] *= i
		j *= i
sum = [1] * lim
for i in range(2, lim - 2, 2):
	sum[i] = sum[i - 1] + phi[i >> 1] * phi[i + 1]
	sum[i + 1] = sum[i] + phi[i + 1] * phi[i + 2 >> 1]
sum[0] = 0
T = int(input())
for t in range(T):
	(L, R) = map(int, input().split())
	i = (isqrt(1 + 8 * L) - 1) // 2
	j = (isqrt(1 + 8 * R) - 1) // 2
	while i * (i + 1) >> 1 < L:
		i += 1
	print(sum[j] - sum[i - 1])


def phi(n):
	res = 1
	i = 2
	while i * i <= n:
		if n % i == 0:
			res -= res // i
		while n % i == 0:
			n //= i
		i += 1
	if n > 1:
		res -= res // n
	return res

def norm(n):
	return n if n % 2 else n // 2
n = 2 * 10 ** 6
phis = [i for i in range(n + 1)]
isprime = [True for _ in range(n + 1)]
for i in range(2, n + 1):
	if not isprime[i]:
		continue
	for j in range(i * i, n + 1, i):
		isprime[j] = False
	for j in range(i, n + 1, i):
		phis[j] -= phis[j] // i
vals = [0]
for i in range(1, n):
	vals.append(vals[-1] + phis[norm(i)] * phis[norm(i + 1)])

def tri(n):
	return (int((8 * n + 1) ** 0.5) - 1) // 2
for _ in range(int(input())):
	(l, r) = [int(x) for x in input().split()]
	print(vals[tri(r)] - vals[tri(l - 1)])


import sys
from time import time
import bisect

def gcd(a, b):
	while b:
		(a, b) = (b, a % b)
	return a

def solve_bf(L, R):
	sm = 0
	for i in range(L, R + 1):
		for j in range(1, i + 1):
			cnt = 0
			if j * (j + 1) == 2 * i:
				for k in range(1, i + 1):
					if gcd(k, i) == 1:
						print('i=%d j=%d k=%d' % (i, j, k), file=sys.stderr)
						cnt += 1
			sm = sm + cnt
	return sm

def prime_gen(n):
	sieve = [1] * n
	for p in range(2, n):
		if sieve[p]:
			yield p
			for i in range(p + p, n, p):
				sieve[i] = 0

def pfs_gen(n, prms):
	for p in prms:
		if p > n:
			break
		cnt = 0
		while n % p == 0:
			cnt += 1
			n //= p
		if cnt:
			yield (p, cnt)
	if n > 1:
		yield (n, 1)

def totient(n, prms):
	t = 1
	for (p, cnt) in pfs_gen(n, prms):
		t *= pow(p, cnt - 1) * (p - 1)
	return t

def a_gen():
	(a, i) = (1, 1)
	while True:
		yield a
		i += 1
		a += i

def totient_cache(n):
	cache = [1] * n
	for p in prime_gen(n):
		for i in range(p, n, p):
			cache[i] *= p - 1
		q = p * p
		while q < n:
			for i in range(q, n, q):
				cache[i] *= p
			q *= p
	return cache

def create_cache(n):
	tcache = totient_cache(2 * n + 2)
	(aa, stt) = ([0], [0])
	sm = 0
	for i in range(1, n):
		(u, v, w) = (i, 2 * i - 1, 2 * i + 1)
		for (u, v) in ((i, 2 * i - 1), (i, 2 * i + 1)):
			a = u * v
			t = tcache[u] * tcache[v]
			sm += t
			aa.append(a)
			stt.append(sm)
	return (aa, stt)

def solve_slow(l, r, prms):
	it = a_gen()
	a = next(it)
	while a < l:
		a = next(it)
	sm = 0
	while a <= r:
		t = totient(a, prms)
		sm += t
		a = next(it)
	return sm

def solve(l, r, aa, stt):
	i = bisect.bisect_left(aa, l) - 1
	j = bisect.bisect_right(aa, r) - 1
	print('l=%d r=%d i=%d j=%d' % (l, r, i, j), file=sys.stderr)
	return stt[j] - stt[i]
tm = time()
n = pow(10, 12)
k = int(pow(-1 + (1 + 8 * n), 0.5) / 2)
(aa, stt) = create_cache(k)
print('cache generated: %d' % (len(aa),), file=sys.stderr)
print('aa=%s' % (aa[:20],), file=sys.stderr)
print('stt=%s' % (stt[:20],), file=sys.stderr)
t = int(input())
for _ in range(t):
	(l, r) = map(int, input().split())
	print(solve(l, r, aa, stt))
print('time elapsed: %s' % (time() - tm,), file=sys.stderr)


import sys, random
from math import sqrt

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

def function(L, R):
	sum = 0
	for i in range(L, R + 1):
		for j in range(1, i + 1):
			cnt = 0
			if j * (j + 1) == 2 * i:
				for k in range(1, i + 1):
					if gcd(k, i) == 1:
						cnt += 1
			sum = sum + cnt
	return sum

def eulerSieve(N):
	if N == 1:
		return [0, 1]
	limit = int(N ** 0.5) + 1
	sieve = [0] * (N + 1)
	for i in range(3, limit, 2):
		if sieve[i] == 0:
			for j in range(i * i, N + 1, 2 * i):
				sieve[j] = i
	sieve[0:3] = [0, 1, 1]
	for i in range(3, N + 1):
		p = sieve[i]
		if p != 0:
			j = i // p
			sieve[i] = sieve[j] * (p if j % p == 0 else p - 1)
		elif i & 3 == 0:
			sieve[i] = sieve[i >> 1] << 1
		elif i & 1 == 0:
			sieve[i] = sieve[i >> 1]
		else:
			sieve[i] = i - 1
	return sieve

def solve():

	def f(n):
		return int((-1 + (1 + 8 * n) ** 0.5) / 2.0)
	N = 1420000
	phis = eulerSieve(N + 1)
	phi2 = [0] * (N + 1)
	for i in range(1, N + 1):
		j = i * (i + 1)
		phi2[i] = phis[i] * phis[i + 1]
		if not j & 3:
			phi2[i] //= 2
	for j in range(1, N + 1):
		phi2[j] += phi2[j - 1]
	stdin = sys.stdin
	MOD = 10 ** 9 + 7
	T = int(stdin.readline())
	for _ in range(T):
		(L, R) = map(int, stdin.readline().split())
		L = f(L - 1)
		R = f(R)
		print(phi2[R] - phi2[L])
solve()


from math import *

def computeTotient(n):
	phi = [i for i in range(n + 1)]
	for p in range(2, n + 1):
		if phi[p] == p:
			phi[p] = p - 1
			for i in range(2 * p, n + 1, p):
				phi[i] = phi[i] // p * (p - 1)
	return phi
ph = computeTotient(2 * 10 ** 6)
a = [0]

def f(r):
	n = 1
	s = n * (n + 1) // 2
	res = 0
	while s <= r:
		if n & 1:
			res += ph[n] * ph[(n + 1) // 2]
		else:
			res += ph[n // 2] * ph[n + 1]
		a.append(res)
		n = n + 1
		s = s + n
	return
f(10 ** 12)
for t in range(int(input())):
	(l, r) = map(int, input().split())
	ll = floor(sqrt(2 * l))
	if ll * (ll + 1) >= 2 * l:
		ll = ll - 1
	rr = floor(sqrt(2 * r))
	if rr * (rr + 1) > 2 * r:
		rr = rr - 1
	print(a[rr] - a[ll])


import math
import bisect
M = 1500000
p = [i for i in range(M + 1)]
vis = [False] * (M + 1)
for i in range(2, M + 1):
	if not vis[i]:
		for j in range(i, M + 1, i):
			p[j] = p[j] * (i - 1) // i
			vis[j] = True
l = [0]
r = [0]
for i in range(M):
	if i % 2:
		a = p[(i + 1) // 2]
		b = p[i]
		c = math.gcd((i + 1) // 2, i)
	else:
		a = p[i // 2]
		b = p[i + 1]
		c = math.gcd(i // 2, i + 1)
	d = p[c]
	l.append(a * b * c // d)
	r.append(i * (i + 1) // 2)
for i in range(1, M):
	l[i] += l[i - 1]
t = int(input())
for _ in range(t):
	lr = input().split()
	L = int(lr[0])
	R = int(lr[1])
	i = bisect.bisect_left(r, L)
	j = bisect.bisect_left(r, R)
	if r[j] != R:
		j -= 1
	print(l[j] - l[i - 1])


import sys

def prepare(N):
	g = [1] * N
	for i in range(2, N):
		if g[i]:
			for x in range(i << 1, N, i):
				g[x] = 0
	g = [i for (i, c) in enumerate(g) if c]
	z = [i for i in range(N + 1)]
	for n in g[2:]:
		for i in range(n, N, n):
			z[i] = z[i] * (n - 1) // n
	g = []
	p1 = p = 0
	for i in range(N):
		(p0, p1) = (p1, i + 1 >> 1 if i & 1 else i + 1)
		p += z[p0] * z[p1]
		g.append((i * (i + 1) >> 1, p))
	return g

def get_less_or_equal(nbr, array):
	(lo, hi) = (0, len(array))
	while lo < hi:
		mid = lo + hi >> 1
		if array[mid][0] < nbr:
			lo = mid + 1
		else:
			hi = mid
	return lo if array[lo][0] == nbr else lo - 1

def weird_function(line):
	r = prepare(1414215)
	q = []
	for w in range(int(next(line))):
		(a, b) = map(int, next(line).split())
		q.append(r[get_less_or_equal(b, r)][1] - r[get_less_or_equal(a - 1, r)][1])
	print('\n'.join(map(str, q)))
weird_function(sys.stdin)


def isqrt(x):
	if x < 0:
		raise ValueError('negative')
	b = x.bit_length() - 1 & -2
	y = 0
	while b >= 0:
		y <<= 1
		z = y + 1
		if z * z <= x >> b:
			y = z
		b -= 2
	return y
lim = 1500000
phi = [1] * lim
j = 4
while j < lim:
	for k in range(j, lim, j):
		phi[k] <<= 1
	j <<= 1
for i in range(3, lim, 2):
	if phi[i] > 1:
		continue
	for k in range(i, lim, i):
		phi[k] *= i - 1
	j = i * i
	while j < lim:
		for k in range(j, lim, j):
			phi[k] *= i
		j *= i
sum = [1] * lim
for i in range(2, lim - 2, 2):
	sum[i] = sum[i - 1] + phi[i >> 1] * phi[i + 1]
	sum[i + 1] = sum[i] + phi[i + 1] * phi[i + 2 >> 1]
sum[0] = 0
T = int(input())
for t in range(T):
	(L, R) = map(int, input().split())
	i = (isqrt(1 + 8 * L) - 1) // 2
	j = (isqrt(1 + 8 * R) - 1) // 2
	while i * (i + 1) >> 1 < L:
		i += 1
	print(sum[j] - sum[i - 1])


import sys, random
from math import sqrt

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

def function(L, R):
	sum = 0
	for i in range(L, R + 1):
		for j in range(1, i + 1):
			cnt = 0
			if j * (j + 1) == 2 * i:
				for k in range(1, i + 1):
					if gcd(k, i) == 1:
						cnt += 1
			sum = sum + cnt
	return sum

def eulerSieve(N):
	if N == 1:
		return [0, 1]
	limit = int(N ** 0.5) + 1
	sieve = [0] * (N + 1)
	for i in range(3, limit, 2):
		if sieve[i] == 0:
			for j in range(i * i, N + 1, 2 * i):
				sieve[j] = i
	sieve[0:3] = [0, 1, 1]
	for i in range(3, N + 1):
		p = sieve[i]
		if p != 0:
			j = i // p
			sieve[i] = sieve[j] * (p if j % p == 0 else p - 1)
		elif i & 3 == 0:
			sieve[i] = sieve[i >> 1] << 1
		elif i & 1 == 0:
			sieve[i] = sieve[i >> 1]
		else:
			sieve[i] = i - 1
	return sieve

def solve():

	def f(n):
		return int((-1 + (1 + 8 * n) ** 0.5) / 2.0)
	N = 1420000
	phis = eulerSieve(N + 1)
	phi2 = [0] * (N + 1)
	for i in range(1, N + 1):
		j = i * (i + 1)
		phi2[i] = phis[i] * phis[i + 1]
		if not j & 3:
			phi2[i] //= 2
	for j in range(1, N + 1):
		phi2[j] += phi2[j - 1]
	stdin = sys.stdin
	MOD = 10 ** 9 + 7
	T = int(stdin.readline())
	for _ in range(T):
		(L, R) = map(int, stdin.readline().split())
		L = f(L - 1)
		R = f(R)
		print(phi2[R] - phi2[L])
solve()