Let's call a non-empty sequence of positive integers a_1, a_2... a_{k} coprime if the greatest common divisor of all elements of this sequence is equal to 1.

Given an array a consisting of n positive integers, find the number of its coprime subsequences. Since the answer may be very large, print it modulo 10^9 + 7.

Note that two subsequences are considered different if chosen indices are different. For example, in the array [1, 1] there are 3 different subsequences: [1], [1] and [1, 1].


-----Input-----

The first line contains one integer number n (1 ≤ n ≤ 100000).

The second line contains n integer numbers a_1, a_2... a_{n} (1 ≤ a_{i} ≤ 100000).


-----Output-----

Print the number of coprime subsequences of a modulo 10^9 + 7.


-----Examples-----
Input
3
1 2 3

Output
5

Input
4
1 1 1 1

Output
15

Input
7
1 3 5 15 3 105 35

Output
100



-----Note-----

In the first example coprime subsequences are:   1  1, 2  1, 3  1, 2, 3  2, 3 

In the second example all subsequences are coprime.

Difficulty: HARD
Tags: Number theory

Bit manipulation

Combinatorics
Skill types: Bit manipulation

Starter code:

Solutions:
import sys, math, queue
MOD = 10 ** 9 + 7
sys.setrecursionlimit(1000000)

def getMul(x):
	a = 1
	for xi in x:
		a *= xi
	return a
n = int(input())
a = list(map(int, input().split()))
d = {}
for ai in a:
	if ai in d:
		d[ai] += 1
	else:
		d[ai] = 1
f = [[] for i in range(max(a) + 10)]
for i in range(1, len(f)):
	for j in range(i, len(f), i):
		f[j].append(i)
seq = [0 for i in range(max(a) + 10)]
for ai in d:
	for fi in f[ai]:
		seq[fi] += d[ai]
for i in range(len(seq)):
	seq[i] = (pow(2, seq[i], MOD) - 1 + MOD) % MOD
pf = [[] for i in range(max(a) + 10)]
pf[0] = None
pf[1].append(1)
for i in range(2, len(f)):
	if len(pf[i]) == 0:
		for j in range(i, len(pf), i):
			pf[j].append(i)
for i in range(1, len(pf)):
	mul = getMul(pf[i])
	if mul == i:
		if len(pf[i]) & 1 == 1:
			pf[i] = -1
		else:
			pf[i] = 1
	else:
		pf[i] = 0
pf[1] = 1
ans = 0
for i in range(1, len(seq)):
	ans += seq[i] * pf[i]
	ans = (ans + MOD) % MOD
print(ans)


n = int(input())
r = list(map(int, input().split()))
dp = [0] * (10 ** 5 + 1)
cnt = [0] * (10 ** 5 + 1)
tmp = [0] * (10 ** 5 + 1)
mod = 10 ** 9 + 7
for i in range(n):
	cnt[r[i]] += 1
for i in range(1, 10 ** 5 + 1):
	for j in range(2 * i, 10 ** 5 + 1, i):
		cnt[i] += cnt[j]
	tmp[i] = pow(2, cnt[i], mod) - 1
for i in range(10 ** 5, 0, -1):
	for j in range(2 * i, 10 ** 5 + 1, i):
		tmp[i] = (tmp[i] - tmp[j]) % mod
print(tmp[1] % mod)


mod = int(1000000000.0 + 7)
n = int(input())
a = [int(_) for _ in input().split()]
freq = {i: 0 for i in range(100001)}
power = {0: 1}
for i in range(1, 100001):
	power[i] = 2 * power[i - 1] % mod
for v in a:
	freq[v] += 1
dp = {i: 0 for i in range(100001)}
for gcd in range(100000, 0, -1):
	mult = 2
	total = freq[gcd]
	complement = 0
	while mult * gcd <= 100000:
		total += freq[mult * gcd]
		complement += dp[mult * gcd]
		mult += 1
	dp[gcd] = (power[total] - 1 - complement + mod) % mod
print(dp[1])


N = 10 ** 5 + 5
MOD = 10 ** 9 + 7
freq = [0 for i in range(N)]
p2 = [0 for i in range(N)]
p2[0] = 1
for i in range(1, N):
	p2[i] = p2[i - 1] * 2
	p2[i] %= MOD

def Calculate_Mobius(N):
	arr = [1 for i in range(N + 1)]
	prime_count = [0 for i in range(N + 1)]
	mobius_value = [0 for i in range(N + 1)]
	for i in range(2, N + 1):
		if prime_count[i] == 0:
			for j in range(i, N + 1, i):
				prime_count[j] += 1
				arr[j] = arr[j] * i
	for i in range(1, N + 1):
		if arr[i] == i:
			if prime_count[i] & 1 == 0:
				mobius_value[i] = 1
			else:
				mobius_value[i] = -1
		else:
			mobius_value[i] = 0
	return mobius_value
mobius = Calculate_Mobius(N)
n = int(input())
b = [int(i) for i in input().split()]
for i in b:
	freq[i] += 1
ans = 0
for i in range(1, N):
	cnt = 0
	for j in range(i, N, i):
		cnt += freq[j]
	total_subsequences = p2[cnt] - 1
	ans = (ans + mobius[i] * total_subsequences % MOD) % MOD
ans += MOD
print(ans % MOD)


import math
import collections

def do():
	n = int(input())
	nums = map(int, input().split(' '))
	count = collections.defaultdict(int)
	for num in nums:
		for i in range(1, int(math.sqrt(num)) + 1):
			cp = num // i
			if num % i == 0:
				count[i] += 1
			if cp != i and num % cp == 0:
				count[cp] += 1
	maxk = max(count.keys())
	freq = {k: (1 << count[k]) - 1 for k in count}
	for k in sorted(count.keys(), reverse=True):
		for kk in range(k << 1, maxk + 1, k):
			freq[k] -= freq[kk] if kk in freq else 0
	return freq[1] % (10 ** 9 + 7)
print(do())


MOD = int(1000000000.0) + 7
N = int(input())
A = list(map(int, input().split()))
pow2 = [pow(2, i, MOD) for i in range(N + 1)]
maxa = max(A)
mcnt = [0 for i in range(maxa + 1)]
mans = [0 for i in range(maxa + 1)]
for i in range(N):
	mcnt[A[i]] += 1
for i in range(1, maxa + 1):
	for j in range(i + i, maxa + 1, i):
		mcnt[i] += mcnt[j]
	mans[i] = pow2[mcnt[i]] - 1
for i in range(maxa, 0, -1):
	for j in range(i + i, maxa + 1, i):
		mans[i] = (mans[i] - mans[j]) % MOD
print(mans[1] + (mans[1] < 0) * MOD)


import sys
mod = 10 ** 9 + 7

def solve():
	n = int(input())
	a = [int(i) for i in input().split()]
	cnt = [0] * (10 ** 5 + 1)
	pat = [0] * (10 ** 5 + 1)
	p2 = [1] * (n + 1)
	for i in range(1, n + 1):
		p2[i] = 2 * p2[i - 1] % mod
	for ai in a:
		cnt[ai] += 1
	for i in range(1, 10 ** 5 + 1):
		for j in range(2 * i, 10 ** 5 + 1, i):
			cnt[i] += cnt[j]
		pat[i] = p2[cnt[i]] - 1
	for i in range(10 ** 5, 0, -1):
		for j in range(2 * i, 10 ** 5 + 1, i):
			pat[i] = (pat[i] - pat[j]) % mod
	ans = pat[1] % mod
	print(ans)

def __starting_point():
	solve()
__starting_point()


import os
import sys
from io import BytesIO, IOBase
BUFSIZE = 8192

class FastIO(IOBase):
	newlines = 0

	def __init__(self, file):
		self._fd = file.fileno()
		self.buffer = BytesIO()
		self.writable = 'x' in file.mode or 'r' not in file.mode
		self.write = self.buffer.write if self.writable else None

	def read(self):
		while True:
			b = os.read(self._fd, max(os.fstat(self._fd).st_size, BUFSIZE))
			if not b:
				break
			ptr = self.buffer.tell()
			(self.buffer.seek(0, 2), self.buffer.write(b), self.buffer.seek(ptr))
		self.newlines = 0
		return self.buffer.read()

	def readline(self):
		while self.newlines == 0:
			b = os.read(self._fd, max(os.fstat(self._fd).st_size, BUFSIZE))
			self.newlines = b.count(b'\n') + (not b)
			ptr = self.buffer.tell()
			(self.buffer.seek(0, 2), self.buffer.write(b), self.buffer.seek(ptr))
		self.newlines -= 1
		return self.buffer.readline()

	def flush(self):
		if self.writable:
			os.write(self._fd, self.buffer.getvalue())
			(self.buffer.truncate(0), self.buffer.seek(0))

class IOWrapper(IOBase):

	def __init__(self, file):
		self.buffer = FastIO(file)
		self.flush = self.buffer.flush
		self.writable = self.buffer.writable
		self.write = lambda s: self.buffer.write(s.encode('ascii'))
		self.read = lambda : self.buffer.read().decode('ascii')
		self.readline = lambda : self.buffer.readline().decode('ascii')
(sys.stdin, sys.stdout) = (IOWrapper(sys.stdin), IOWrapper(sys.stdout))
input = lambda : sys.stdin.readline().rstrip('\r\n')
from math import gcd, ceil

def prod(a, mod=10 ** 9 + 7):
	ans = 1
	for each in a:
		ans = ans * each % mod
	return ans

def lcm(a, b):
	return a * b // gcd(a, b)

def binary(x, length=16):
	y = bin(x)[2:]
	return y if len(y) >= length else '0' * (length - len(y)) + y
for _ in range(int(input()) if not True else 1):
	n = int(input())
	a = list(map(int, input().split()))
	mod = 10 ** 9 + 7
	twopow = [1] * (10 ** 5 + 69)
	for i in range(1, 10 ** 5 + 69):
		twopow[i] = twopow[i - 1] * 2 % mod
	count = [0] * 100069
	for i in a:
		count[i] += 1
	multiples = [0] * 100069
	for i in range(1, 10 ** 5 + 1):
		for j in range(i, 10 ** 5 + 1, i):
			multiples[i] += count[j]
	gcd_of = [0] * 100069
	for i in range(10 ** 5, 0, -1):
		gcd_of[i] = (twopow[multiples[i]] - 1) % mod
		for j in range(2 * i, 10 ** 5 + 1, i):
			gcd_of[i] -= gcd_of[j]
	print(gcd_of[1] % mod)


import sys
mod = 10 ** 9 + 7

def solve():
	n = int(input())
	a = [int(i) for i in input().split()]
	cnt = [0] * (10 ** 5 + 1)
	for ai in a:
		for d in range(1, ai + 1):
			if d * d > ai:
				break
			if ai % d == 0:
				if d != ai // d:
					cnt[d] += 1
					cnt[ai // d] += 1
				else:
					cnt[d] += 1
	ans = 0
	for i in range(1, 10 ** 5 + 1):
		ans += mobius(i) * (pow(2, cnt[i], mod) - 1)
		ans %= mod
	print(ans)

def mobius(x):
	assert x >= 1
	divcnt = 0
	for p in range(2, x + 1):
		if p * p > x:
			break
		if x % p != 0:
			continue
		x //= p
		if x % p == 0:
			return 0
		else:
			divcnt ^= 1
	if x > 1:
		divcnt ^= 1
	return (-1) ** divcnt

def __starting_point():
	solve()
__starting_point()


fact = [1]
temp = 1
MOD = 10 ** 9 + 7
for i in range(1, 10 ** 5 + 5):
	temp *= i
	temp %= MOD
	fact += [temp]

def bino(a, b):
	up = fact[a]
	down = pow(fact[b] * fact[a - b], MOD - 2, MOD)
	return up * down % MOD

def find(A):
	MOD = 10 ** 9 + 7
	dp = [0] * (10 ** 5 + 2)
	for x in A:
		dp[x] += 1
	for i in range(2, len(dp)):
		for j in range(2, len(dp)):
			if i * j > len(dp) - 1:
				break
			dp[i] += dp[i * j]
	for i in range(2, len(dp)):
		dp[i] = (pow(2, dp[i], MOD) - 1) % MOD
	for i in range(len(dp) - 1, 1, -1):
		for j in range(2, len(dp)):
			if i * j >= len(dp):
				break
			dp[i] -= dp[i * j]
			dp[i] %= MOD
	ans = 0
	for i in range(2, len(dp)):
		ans += dp[i]
		ans %= MOD
	return (pow(2, len(A), MOD) - ans - 1) % MOD
input()
print(find(list(map(int, input().strip().split(' ')))))


n = int(input())
L = list(map(int, input().split(' ')))
D = [0] * 101000
mod = 10 ** 9 + 7
itt = [0] * 101000
p = [0] * 100010
D[0] = 1
for i in range(100010):
	D[i + 1] = D[i] * 2 % mod
for i in range(n):
	itt[L[i]] += 1
for i in range(1, 100001):
	for j in range(i * 2, 100001, i):
		itt[i] += itt[j]
	p[i] = (D[itt[i]] + mod - 1) % mod
i = 100000
while i >= 1:
	for j in range(i * 2, 100001, i):
		p[i] -= p[j]
	p[i] = (p[i] % mod + mod) % mod
	i -= 1
print(p[1])