The Golomb sequence $G_1, G_2, \ldots$ is a non-decreasing integer sequence such that for each positive integer $n$, $G_n$ is the number of occurrences of $n$ in this sequence. The first few elements of $G$ are $[1, 2, 2, 3, 3, 4, 4, 4, 5, \ldots]$. Do you know the recurrence relation for the Golomb sequence? It is $G_1 = 1$ and $G_{n+1} = 1+G_{n+1-G_{G_n}}$ for each $n \ge 1$. A self-describing sequence, isn't it?
Mr. Strange wants to learn CP, so he asked Chef, who is one of the best competitive programmers in the world, to teach him. Chef decided to test his ability by giving him the following task.
Find the sum of squares of the $L$-th through $R$-th term of the Golomb sequence, i.e. $S = \sum_{i=L}^R G_i^2$. Since the sum can be quite large, compute it modulo $10^9+7$.
Can you help Mr. Strange carry out this task given to him by his teacher?

-----Input-----
- The first line of the input contains a single integer $T$ denoting the number of test cases. The description of $T$ test cases follows.
- The first and only line of each test case contains two space-separated integers $L$ and $R$.

-----Output-----
For each test case, print a single line containing one integer $S$ modulo $10^9+7$.

-----Constraints-----
- $1 \le T \le 10^5$
- $1 \le L \le R \le 10^{10}$

-----Subtasks-----
Subtask #1 (50 points):
- $T \le 10^2$
- $R \le 10^9$
Subtask #2 (50 points): original constraints

-----Example Input-----
3
1 5
2 4
100 100

-----Example Output-----
27
17
441

-----Explanation-----
Example case 1: $1^2 + 2^2 + 2^2 + 3^2 + 3^2 = 27$
Example case 2: $2^2 + 2^2 + 3^2 = 17$
Example case 3: $21^2 = 441$

Difficulty: VERY_HARD
Tags: Sorting

Range queries

Greedy algorithms

Complete search

Mathematics
Skill types: Complete search

Sorting

Range queries

Greedy algorithms

Starter code:

Solutions:
import sys
n = 2000000
M = 10 ** 9 + 7
ans = 1
dp = [0] * (n + 1)
a = [0, 1]
s = [0, 1]
dp[1] = 1
su = 1
for i in range(2, n + 1):
	dp[i] = 1 + dp[i - dp[dp[i - 1]]]
	ans += dp[i]
	a.append(ans)
	su = (su + i ** 2 * dp[i]) % M
	s.append(su)
	if ans > 10 ** 10:
		break
import bisect as bs
from sys import stdin, stdout, setrecursionlimit
T = int(stdin.readline())
for _ in range(T):
	(l, r) = list(map(int, stdin.readline().split()))
	l -= 1
	lb = bs.bisect_right(a, l)
	lb -= 1
	left = s[lb]
	left = (left + (l - a[lb]) * (lb + 1) * (lb + 1)) % M
	rb = bs.bisect_right(a, r)
	rb -= 1
	right = s[rb]
	right = (right + (r - a[rb]) * (rb + 1) * (rb + 1)) % M
	print((right - left) % M)


import bisect
MAXR = 100
MAXN = 20000000
T = int(input().strip())
queries = []
for t in range(T):
	(L, R) = map(int, input().strip().split())
	queries.append((L, R))
	MAXR = max(MAXR, R + 1)
MOD = 10 ** 9 + 7
g = [0, 1, 2]
p = [0, 1, 3]
s = [0, 1, 9]
for n in range(3, MAXN):
	gg = 1 + g[n - g[g[n - 1]]]
	pp = p[n - 1] + gg
	ss = (s[n - 1] + gg * n * n) % MOD
	g.append(gg)
	p.append(pp)
	s.append(ss)
	if pp > MAXR:
		break

def process(m):
	n = bisect.bisect_right(p, m)
	return (s[n - 1] + (m - p[n - 1]) * n * n) % MOD
for (L, R) in queries:
	print((process(R) - process(L - 1)) % MOD)


import bisect as bi
temp_g = [0, 1, 2]
temp_s = [0, 1, 9]
temp_p = [0, 1, 3]
max_mod = 10 ** 9 + 7

def find_gn(n):
	return 1 + temp_g[n - temp_g[temp_g[n - 1]]]

def find_sn(n, gn):
	return (temp_s[n - 1] + gn * n * n) % max_mod

def find_pn(n, gn):
	return temp_p[n - 1] + gn

def bisection(x):
	n = bi.bisect_right(temp_p, x)
	return (temp_s[n - 1] + (x - temp_p[n - 1]) * n * n) % max_mod
(lr, max_) = ([], 100)
for _ in range(int(input())):
	(l, r) = map(int, input().split())
	lr.append((l, r))
	max_ = max(max_, r + 1)
for n in range(3, 2 * 10 ** 8, 1):
	gn = find_gn(n)
	sn = find_sn(n, gn)
	pn = find_pn(n, gn)
	temp_g.append(gn)
	temp_p.append(pn)
	temp_s.append(sn)
	if pn > max_:
		break
for pair in lr:
	print((bisection(pair[1]) - bisection(pair[0] - 1)) % max_mod)


def find_left(ar, x, l, r):
	mid = (l + r) // 2
	if l > r:
		return l
	if ar[mid] == x:
		return mid
	elif ar[mid] > x:
		return find_left(ar, x, l, mid - 1)
	else:
		return find_left(ar, x, mid + 1, r)
MOD = 10 ** 9 + 7
G = [1, 2, 2]
right = [1, 3, 5]
sq_sum = [1, 9, 27]
next = 3
count = G[next - 1]
num = 4
while right[-1] < 10 ** 10:
	while count > 0:
		G.append(next)
		count -= 1
		right.append(right[-1] + G[-1])
		sq_sum.append(sq_sum[-1] + num ** 2 * G[-1])
		sq_sum[-1] %= MOD
		num += 1
	next += 1
	count = G[next - 1]
for _ in range(int(input())):
	(l, r) = map(int, input().split())
	(a, b) = (find_left(right, l, 0, len(right)), find_left(right, r, 0, len(right)))
	ans = 0
	if a == b:
		ans += (a + 1) ** 2 * (r - l + 1)
		ans = ans % MOD
	elif a + 1 == b:
		ans += (a + 1) ** 2 * (right[a] - l + 1)
		ans = ans % MOD
		ans += (b + 1) ** 2 * (r - right[a])
		ans = ans % MOD
	else:
		ans += (a + 1) ** 2 * (right[a] - l + 1)
		ans = ans % MOD
		ans += (b + 1) ** 2 * (r - right[b - 1])
		ans = ans % MOD
		ans += sq_sum[b - 1] - sq_sum[a]
		ans = ans % MOD
	print(ans)


def find_upper_bound(arr, key):
	(low, high) = (0, len(arr) - 1)
	while low <= high:
		mid = (low + high) // 2
		if arr[mid] == key:
			return mid
		elif arr[mid] > key and mid - 1 >= 0 and (arr[mid - 1] < key):
			return mid
		elif arr[mid] > key:
			high = mid - 1
		else:
			low = mid + 1
	return mid

def get_query(l):
	nonlocal prefix_storer, bin_storer
	ind = find_upper_bound(bin_storer, l)
	surplus = abs(bin_storer[ind] - l) * ind * ind % limit
	return (prefix_storer[ind] - surplus + limit) % limit

def fire_query(l, r):
	return (get_query(r) - get_query(l - 1) + limit) % limit
(golomb, dp, prefix_storer, bin_storer) = ([], [0, 1], [0, 1], [0, 1])
limit = 10 ** 9 + 7
for i in range(2, 10 ** 6 + 100):
	dp.append(1 + dp[i - dp[dp[i - 1]]])
	bin_storer.append(dp[-1] + bin_storer[-1])
	prefix_storer.append((prefix_storer[-1] + dp[-1] * i * i % limit) % limit)
for _ in range(int(input())):
	(l, r) = map(int, input().split())
	print(fire_query(l, r))


import bisect
import sys
input = sys.stdin.readline
MAXR = 100
MAXN = 20000000
t = int(input().strip())
queries = []
for t in range(t):
	(L, R) = map(int, input().strip().split())
	queries.append((L, R))
	MAXR = max(MAXR, R + 1)
MOD = 10 ** 9 + 7
g = [0, 1, 2]
p = [0, 1, 3]
s = [0, 1, 9]
for n in range(3, MAXN):
	gn = 1 + g[n - g[g[n - 1]]]
	pn = p[n - 1] + gn
	sn = (s[n - 1] + gn * n * n) % MOD
	g.append(gn)
	p.append(pn)
	s.append(sn)
	if pn > MAXR:
		break

def process(m):
	n = bisect.bisect_right(p, m)
	return (s[n - 1] + (m - p[n - 1]) * n * n) % MOD
for (L, R) in queries:
	print((process(R) - process(L - 1)) % MOD)


import sys, math
import io, os
from bisect import bisect_left as bl, bisect_right as br, insort
from heapq import heapify, heappush, heappop
from collections import defaultdict as dd, deque, Counter

def data():
	return sys.stdin.readline().strip()

def mdata():
	return list(map(int, data().split()))

def outl(var):
	sys.stdout.write('\n'.join(map(str, var)) + '\n')

def out(var):
	sys.stdout.write(str(var) + '\n')
from fractions import Fraction
INF = float('inf')
mod = int(1000000000.0) + 7
a = [(0, 1), (1, 2)]
b = [0, 1]
x = 10 ** 10
ans = 1
cnt = 1
d = dd(int)
d[2] = 2
k = 2
i = 2
while cnt <= x:
	d[i] = k
	if cnt + d[i] > x:
		ans = (ans + i ** 2 * (x - cnt)) % mod
		break
	ans = (ans + i ** 2 * d[i]) % mod
	cnt += d[i]
	d[k] -= 1
	i += 1
	if d[k] == 0:
		k += 1
	a.append((ans, i))
	b.append(cnt)
ans1 = []
inp = [int(x) for x in sys.stdin.read().split()]
ii = 0
T = inp[ii]
ii += 1
for t in range(T):
	l = inp[ii]
	ii += 1
	r = inp[ii]
	ii += 1
	i = bl(b, r)
	if i == len(b) or b[i] > r:
		i -= 1
	ans = a[i][0] + a[i][1] ** 2 * (r - b[i])
	if l != 1:
		i = bl(b, l - 1)
		if i == len(b) or b[i] > l - 1:
			i -= 1
		ans -= a[i][0] + a[i][1] ** 2 * (l - 1 - b[i])
	ans1.append((ans + mod) % mod)
outl(ans1)


import bisect
M = 1000000007
G = [0, 1, 2, 2]
currPos = 4
for i in range(3, 2000000):
	if currPos > 2000000:
		break
	j = 0
	while j < G[i] and currPos < 2000000:
		G.append(i)
		currPos += 1
		j += 1
prefixSum1 = [0]
prefixSum2 = [0]
for i in range(1, 2000000):
	prefixSum1.append(prefixSum1[i - 1] + G[i])
for i in range(1, 2000000):
	prefixSum2.append((prefixSum2[i - 1] + i * i % M * G[i] % M) % M)

def solve(n):
	nthterm = bisect.bisect_left(prefixSum1, n, lo=0, hi=len(prefixSum1)) - 0
	ans = 0
	if nthterm > 0:
		ans = prefixSum2[nthterm - 1]
	ans = (ans + nthterm * nthterm % M * (n - prefixSum1[nthterm - 1]) % M) % M
	return ans
for tc in range(int(input())):
	(l, r) = map(int, input().split())
	print((solve(r) - solve(l - 1) + M) % M)


import bisect

def f(R):
	z = bisect.bisect_right(pref, R)
	rem = R - pref[z - 1]
	k = rem * (z + 1) * (z + 1)
	return (k + ans[z - 1]) % 1000000007
n = 10 ** 6
golem = [0 for i in range(3 * n)]
pref = [0 for i in range(3 * n)]
ans = [0 for i in range(3 * n)]
golem[0] = 1
hashi = dict()
cnt = 0
pref[0] = 1
ans[0] = 1
hashi[2] = 2
for i in range(1, 3 * n):
	if cnt == 0:
		golem[i] = 1 + golem[i - 1]
		cnt = hashi[golem[i]]
		cnt -= 1
	else:
		golem[i] = golem[i - 1]
		cnt -= 1
	hashi[i + 1] = golem[i]
	pref[i] = pref[i - 1] + golem[i]
	ans[i] = ans[i - 1] + golem[i] * (i + 1) * (i + 1)
	ans[i] %= 1000000007
t = int(input())
for you in range(t):
	l = input().split()
	L = int(l[0])
	R = int(l[1])
	if L == 1:
		print(f(R))
	else:
		po = f(R) - f(L - 1)
		print(po % 1000000007)


from sys import stdin, stdout
import math
from bisect import bisect_left
mod = 1000000007

def mul(a, b):
	nonlocal mod
	return a % mod * (b % mod) % mod

def add(a, b):
	nonlocal mod
	return (a % mod + b % mod) % mod
g = [0, 1]
pre = [0, 1]
ans = [0, 1]
i = 2
while True:
	g.append(1 + g[i - g[g[i - 1]]])
	pre.append(pre[i - 1] + g[i])
	ans.append(add(ans[i - 1], mul(mul(i, i), g[i])))
	if pre[i] > 10000000000:
		break
	i += 1
test = int(stdin.readline())
for _ in range(test):
	(l, r) = list(map(int, stdin.readline().split()))
	(sm1, sm2) = (0, 0)
	if l == 1:
		sm1 = 0
	else:
		l -= 1
		tl = bisect_left(pre, l)
		sm1 = add(ans[tl - 1], mul(l - pre[tl - 1], mul(tl, tl)))
	tr = bisect_left(pre, r)
	sm2 = add(ans[tr - 1], mul(r - pre[tr - 1], mul(tr, tr)))
	print((sm2 - sm1) % mod)


import bisect
m = 10 ** 9 + 7
arr = [0, 1, 2, 2]
c = 5
i = 3
while c <= 10 ** 10:
	for j in range(arr[i]):
		arr.append(i)
		c += i
		if c >= 10 ** 10:
			break
	i += 1
n = len(arr)
pref1 = [0] * n
pref2 = [0] * n
for i in range(1, n):
	pref1[i] = pref1[i - 1] + arr[i]
for i in range(1, n):
	pref2[i] = pref2[i - 1] + i * i * arr[i] % m
	pref2[i] %= m

def solve(x):
	bi = bisect.bisect_left(pref1, x)
	if pref1[bi] == x:
		return pref2[bi]
	ans = x - pref1[bi - 1]
	return (pref2[bi - 1] + bi * bi * ans) % m
for _ in range(int(input())):
	(l, r) = map(int, input().split())
	print((solve(r) - solve(l - 1)) % m)


from bisect import bisect_right
mod = 10 ** 9 + 7
(mr, mn, q) = (100, 2 * 10 ** 7, [])
(g, p, s) = ([0, 1, 2], [0, 1, 3], [0, 1, 9])
for _ in range(int(input())):
	(l, r) = map(int, input().split())
	q.append((l, r))
	mr = max(mr, r + 1)
for n in range(3, mn):
	gg = 1 + g[n - g[g[n - 1]]]
	pp = p[n - 1] + gg
	ss = (s[n - 1] + gg * n * n) % mod
	g.append(gg)
	p.append(pp)
	s.append(ss)
	if pp > mr:
		break

def pro(m):
	n = bisect_right(p, m)
	return (s[n - 1] + (m - p[n - 1]) * n * n) % mod
for (l, r) in q:
	print((pro(r) - pro(l - 1)) % mod)


from bisect import bisect_left

def solve(g, prefixGolombArr, prefixFreq, size, mod):
	for i in range(1, size):
		g.append((i + 1, 1 + g[i + 1 - g[g[i][1]][1]][1]))
		prefixGolombArr.append((prefixGolombArr[i - 1] + g[i + 1][0] * g[i + 1][0] * g[i + 1][1] % mod) % mod)
		prefixFreq.append(g[i + 1][1] + prefixFreq[i - 1])

def findAnsLimit(X, g, prefixFreq, prefixGolombArr, mod):
	K = bisect_left(prefixFreq, X)
	if prefixFreq[K] == X:
		return prefixGolombArr[K]
	else:
		return (prefixGolombArr[K - 1] + g[K + 1][0] * g[K + 1][0] * (X - prefixFreq[K - 1]) % mod) % mod

def main():
	g = [(0, 1), (1, 1)]
	size = 2000009
	mod = int(1000000000.0 + 7)
	prefixGolombArr = [g[1][0] * g[1][0] % mod]
	prefixFreq = [1]
	solve(g, prefixGolombArr, prefixFreq, size, mod)
	T = int(input())
	for _ in range(T):
		(L, R) = map(int, input().split())
		if L == 1:
			print(findAnsLimit(R, g, prefixFreq, prefixGolombArr, mod))
		else:
			print((findAnsLimit(R, g, prefixFreq, prefixGolombArr, mod) - findAnsLimit(L - 1, g, prefixFreq, prefixGolombArr, mod)) % mod)
main()


import sys
import bisect as b
input = sys.stdin.readline
li = []
mod = 10 ** 9 + 7
mxr = 100
mxn = 20000000
for _ in range(int(input())):
	(l, r) = map(int, input().split())
	li.append([l, r])
	mxr = max(mxr, r + 1)
g = [0, 1, 2]
p = [0, 1, 3]
s = [0, 1, 9]
for i in range(3, mxn):
	gg = 1 + g[i - g[g[i - 1]]]
	pp = (p[i - 1] + gg) % mod
	ss = (s[i - 1] + gg * i * i) % mod
	g.append(gg)
	p.append(pp)
	s.append(ss)
	if pp > mxr:
		break

def solve(m):
	n = b.bisect_right(p, m)
	return (s[n - 1] + (m - p[n - 1]) * n * n) % mod
for (l, r) in li:
	ans = (solve(r) - solve(l - 1)) % mod
	print(ans)