Vivek initially has an empty array $a$ and some integer constant $m$.

He performs the following algorithm:  Select a random integer $x$ uniformly in range from $1$ to $m$ and append it to the end of $a$.  Compute the greatest common divisor of integers in $a$.  In case it equals to $1$, break  Otherwise, return to step $1$. 

Find the expected length of $a$. It can be shown that it can be represented as $\frac{P}{Q}$ where $P$ and $Q$ are coprime integers and $Q\neq 0 \pmod{10^9+7}$. Print the value of $P \cdot Q^{-1} \pmod{10^9+7}$.


-----Input-----

The first and only line contains a single integer $m$ ($1 \leq m \leq 100000$).


-----Output-----

Print a single integer — the expected length of the array $a$ written as $P \cdot Q^{-1} \pmod{10^9+7}$.


-----Examples-----
Input
1

Output
1

Input
2

Output
2

Input
4

Output
333333338



-----Note-----

In the first example, since Vivek can choose only integers from $1$ to $1$, he will have $a=[1]$ after the first append operation, and after that quit the algorithm. Hence the length of $a$ is always $1$, so its expected value is $1$ as well.

In the second example, Vivek each time will append either $1$ or $2$, so after finishing the algorithm he will end up having some number of $2$'s (possibly zero), and a single $1$ in the end. The expected length of the list is $1\cdot \frac{1}{2} + 2\cdot \frac{1}{2^2} + 3\cdot \frac{1}{2^3} + \ldots = 2$.

Difficulty: HARD
Tags: Number theory

Dynamic programming

Mathematics

Probability
Skill types: Dynamic programming

Starter code:

Solutions:
M = 1000000007
m = int(input())

def get_prm_fctrs(n):
	p = 2
	factors = []
	while p * p <= n:
		if n % p == 0:
			factors.append(p)
		while n % p == 0:
			n //= p
		p += 1
	if n > 1:
		factors.append(n)
	return factors

def get_n_coprimes(n, m):
	if n == 1:
		return m
	l = len(prm_factors[n])
	ans = 0
	for i in range(1, 2 ** l):
		prod = 1
		cnt = 0
		for j in range(l):
			if 1 << j & i > 0:
				prod *= prm_factors[n][j]
				cnt += 1
		if cnt & 1 == 1:
			ans -= m // prod
		else:
			ans += m // prod
	return m + ans
prm_factors = [[0], [0]]
for i in range(2, m + 1):
	prm_factors.append(get_prm_fctrs(i))
dp = [m for i in range(m + 1)]
for i in range(2, m + 1):
	prod = 2
	while prod * prod < i:
		if i % prod == 0:
			dp[i] += dp[prod] * get_n_coprimes(i // prod, m // prod)
			dp[i] += dp[i // prod] * get_n_coprimes(prod, m // (i // prod))
		prod += 1
	if prod * prod == i:
		dp[i] += dp[prod] * get_n_coprimes(i // prod, m // prod)
	dp[i] = dp[i] * pow(m - m // i, M - 2, M) % M
ans = 0
for i in range(1, m + 1):
	ans = (ans + dp[i]) % M
print(ans * pow(m, M - 2, M) % M)


N = 100010

def gen_mobius_function():
	mu = [1] * N
	mu[0] = 0
	P = [True] * N
	P[0] = P[1] = False
	for i in range(2, N):
		if P[i]:
			j = i
			while j < N:
				P[j] = False
				mu[j] *= -1
				j += i
			j = i * i
			while j < N:
				mu[j] = 0
				j += i * i
	return mu
m = int(input())
mu = gen_mobius_function()
MOD = 10 ** 9 + 7

def mod_inv(x):
	return pow(x, MOD - 2, MOD)
E = 1
for i in range(2, N):
	p = m // i * mod_inv(m)
	E += -mu[i] * p * mod_inv(1 - p)
print(E % MOD)


from math import floor
from functools import reduce
MOD = floor(1000000000.0 + 7)
expected_len = [0, 1]
n = int(input())
factors = [[] for i in range(n + 1)]
prime_factors = [[] for i in range(n + 1)]

def ext_euclid(a, b):
	if b == 0:
		return (1, 0, a)
	(x, y, q) = ext_euclid(b, a % b)
	(x, y) = (y, x - a // b * y)
	return (x, y, q)

def inverse(num):
	return ext_euclid(MOD, num)[1] % MOD
inv = [0] * (n + 1)
for i in range(n + 1):
	inv[i] = inverse(i)
for i in range(1, n + 1):
	prime_fact = False
	if len(prime_factors[i]) < 2:
		prime_factors[i].append(i)
		prime_fact = True
	factors[i].append(i)
	for j in range(2 * i, n + 1, i):
		factors[j].append(i)
		if prime_fact:
			prime_factors[j].append(i)

def f(x, y):
	remain = 0
	new_n = n // x
	new_y = reduce(lambda x, y: x * y, prime_factors[y])
	for fact in factors[new_y]:
		if fact != 1:
			if len(prime_factors[fact]) & 1:
				remain -= new_n // fact
			else:
				remain += new_n // fact
	return new_n - remain
for i in range(2, n + 1):
	e_len = 0
	for ele in factors[i]:
		if ele != i:
			e_len += f(ele, i // ele) * expected_len[ele] * inv[n] % MOD
	e_len = (e_len + 1) * n * inv[n - f(i, 1)] % MOD
	expected_len.append(e_len)
print(sum(expected_len) * inv[n] % MOD)


from sys import *
m = int(input())
q = [0] * (m + 1)
c = 1
for i in range(m, 1, -1):
	w = m // i * pow(m, 1000000007 - 2, 1000000007)
	q[i] = w * pow(1 - w, 1000000007 - 2, 1000000007) % 1000000007
	for j in range(2 * i, m + 1, i):
		q[i] = (q[i] - q[j]) % 1000000007
	c = c + q[i]
print(c % 1000000007)


from fractions import Fraction
g = [0] * 111111
t = [0] * 111111
prime = [0] * 111111
pm = [0] * 111111
for i in range(len(prime)):
	prime[i] = {}

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

def ksm(a, b, p=int(1000000000.0 + 7)):
	a = a % p
	ans = 1
	while b > 0:
		if b & 1:
			ans = ans * a % p
		b = b >> 1
		a = a * a % p
	return ans
m = int(input())
p = int(1000000000.0 + 7)
tmp = ksm(m, p - 2)
for j in range(m, 0, -1):
	c1 = 0
	tyouce = tmp
	gyouce = 0
	if j == 1:
		c1 = 0
	else:
		c1 = m // j * tmp % p
	for k in range(2, m // j + 1):

		def get_ck(m, n):
			ans = 0

			def get_fac():
				for i in range(2, 100011):
					if pm[i] == 0:
						prime[i][i] = 1
						for j in range(2, 100001 // i + 1):
							pm[i * j] = 1
							prime[i * j][i] = 1
			if prime[m] == {}:
				get_fac()
			prime_list = list(prime[m].keys())
			for i in range(1, 1 << len(prime_list)):
				t = i
				rongchixishu = -1
				lcm = 1
				index = len(prime_list) - 1
				while t > 0:
					if t & 1:
						lcm = lcm * prime_list[index]
						rongchixishu *= -1
					t >>= 1
					index -= 1
				ans += n // lcm * rongchixishu
			return ans
		ck = (m // j - get_ck(k, m // j) + p) % p
		ck = ck * tmp % p
		tyouce = (tyouce + ck * t[k * j]) % p
		gyouce = (gyouce + ck * (t[k * j] + g[k * j])) % p
	t[j] = tyouce * ksm(1 - c1 + p, p - 2) % p
	gyouce = (gyouce + tmp + c1 * t[j]) % p
	g[j] = gyouce * ksm(1 - c1 + p, p - 2) % p
print(g[1])


max_n = 10 ** 5 + 1
spf = [i for i in range(max_n)]
prime = [True for i in range(max_n)]
mobius = [0 for i in range(max_n)]
prime[0] = prime[1] = False
mobius[1] = 1
primes = []
for i in range(2, max_n):
	if prime[i]:
		spf[i] = i
		mobius[i] = -1
		primes.append(i)
	for j in primes:
		prod = j * i
		if j > spf[i] or prod >= max_n:
			break
		spf[prod] = j
		prime[prod] = False
		if spf[i] == j:
			mobius[prod] = 0
		else:
			mobius[prod] = -mobius[i]
import sys
input = sys.stdin.readline
mod = 10 ** 9 + 7
m = int(input())
ans = 1
for i in range(2, m + 1):
	n = m // i
	ans = (ans + -mobius[i] * n * pow(m - n, mod - 2, mod)) % mod
print(ans)


big = 100010

def gen_mu():
	mu = [1] * big
	mu[0] = 0
	P = [True] * big
	P[0] = P[1] = False
	for i in range(2, big):
		if P[i]:
			j = i
			while j < big:
				P[j] = False
				mu[j] *= -1
				j += i
			j = i * i
			while j < big:
				mu[j] = 0
				j += i * i
	return mu
m = int(input())
mu = gen_mu()
MOD = 10 ** 9 + 7

def mod_inv(x):
	return pow(x, MOD - 2, MOD)
s = 1
for i in range(2, big):
	p = m // i * mod_inv(m)
	s += -mu[i] * p * mod_inv(1 - p)
print(s % MOD)


MOD = 10 ** 9 + 7
m = int(input())
f = [0] * (m + 1)
ans = 1
for i in range(m, 1, -1):
	p = m // i * pow(m, MOD - 2, MOD)
	f[i] = p * pow(1 - p, MOD - 2, MOD) % MOD
	for j in range(2 * i, m + 1, i):
		f[i] = (f[i] - f[j]) % MOD
	ans += f[i]
print(ans % MOD)


import sys, os, math, re, time, random, copy, shutil, itertools, functools, heapq
from bisect import bisect_left, bisect_right
from collections import deque, defaultdict, Counter, OrderedDict
MOD = int(1000000000.0 + 7)

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

def mpow(x, n):
	r = 1
	while n > 0:
		if n & 1 == 1:
			r = r * x % MOD
		(x, n) = (x * x % MOD, n >> 1)
	return r
m = int(input())
d = [[] for _ in range(m + 1)]
for i in range(1, m + 1):
	for j in range(i + i, m + 1, i):
		d[j].append(i)
pr = [False, False] + [True] * (m - 1)
fa = [{} for _ in range(m + 1)]
for i in range(2, m + 1):
	if not pr[i]:
		continue
	for j in range(i, m + 1, i):
		pr[j] = False
		(k, c) = (j, 0)
		while k % i == 0:
			(k, c) = (k // i, c + 1)
		fa[j][i] = c

def calc(g, x):
	p = []
	for (k, v) in fa[x].items():
		if k not in fa[g]:
			p.append(k)
		elif fa[g][k] != v:
			p.append(k)
	sz = len(p)
	n = m // g
	ret = m // g
	for i in range(1, 1 << sz):
		sign = z = 1
		for j in range(sz):
			if i & 1 << j != 0:
				z *= p[j]
				sign *= -1
		ret += sign * (n // z)
	return ret
dp = [0] * (m + 1)
mm = mpow(m, MOD - 2)
for i in range(2, m + 1):
	r = 1
	for x in d[i]:
		c = calc(x, i)
		r = (r + dp[x] * c % MOD * mm % MOD) % MOD
	c = m // i
	dp[i] = r * m % MOD * mpow(m - c, MOD - 2) % MOD
ans = 1
for i in range(1, m + 1):
	ans = (ans + dp[i] * mm % MOD) % MOD
print(ans)


def ksm(a, b, p):
	a = a % p
	b = b % (p - 1)
	res = 1
	mul = a
	while b > 0:
		if b % 2 == 1:
			res = res * mul % p
		mul = mul * mul % p
		b //= 2
	return res

def inv(a, p=int(1000000000.0 + 7)):
	a = (a % p + p) % p
	return ksm(a, p - 2, p)
m = int(input())
p = int(1000000000.0 + 7)
ans = [0] * 111111
res = 1
for i in range(m, 1, -1):
	total = m // i * inv(m, p) % p
	total = total * inv(1 - total) % p
	ans[i] = total
	for j in range(i + i, m + 1, i):
		ans[i] = ((ans[i] - ans[j]) % p + p) % p
	res = (res + ans[i]) % p
print(res)


import sys, math
input = sys.stdin.readline

def inp():
	return int(input())

def inlt():
	return list(map(int, input().split()))

def insr():
	s = input()
	return s[:len(s) - 1]

def invr():
	return map(int, input().split())
t = 1
for tt in range(t):
	m = int(input())
	q = [0] * (m + 1)
	c = 1
	for i in range(m, 1, -1):
		w = m // i * pow(m, 1000000007 - 2, 1000000007)
		q[i] = w * pow(1 - w, 1000000007 - 2, 1000000007) % 1000000007
		for j in range(2 * i, m + 1, i):
			q[i] = (q[i] - q[j]) % 1000000007
		c = c + q[i]
	print(c % 1000000007)


mod = int(1000000000.0 + 7)

def ksm(a: int, k: int) -> int:
	res = 1
	while k:
		if k & 1:
			res = res * a % mod
		a = a * a % mod
		k >>= 1
	return res
m = int(input())
dp = [0 for i in range(m + 1)]
dp[1] = 0
ans = 0
inv = ksm(m, mod - 2)
for i in range(2, m + 1):
	fac = []
	j = 1
	while j * j <= i:
		if i % j == 0:
			fac.append(j)
			if j != i // j:
				fac.append(i // j)
		j += 1
	fac.sort()
	a = [0 for j in range(len(fac))]
	a[-1] = m // fac[-1]
	for j in range(len(fac) - 2, -1, -1):
		a[j] = m // fac[j]
		for k in range(j + 1, len(fac)):
			if fac[k] % fac[j] == 0:
				a[j] -= a[k]
		dp[i] += (dp[fac[j]] + 1) * a[j] * inv % mod
	dp[i] += a[-1] * inv % mod
	dp[i] = dp[i] * m * ksm(m - a[-1], mod - 2) % mod
	ans = (ans + dp[i] * inv) % mod
print((ans + 1) % mod)


import sys
factors = []
dp = []
mod = int(1000000000.0 + 7)
inv_m = 0

def sub(a, b):
	a -= b
	a %= mod
	if a < 0:
		a += mod
	return a

def add(a, b):
	return (a + b) % mod

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

def pow(a, b):
	if b == 0:
		return 1
	if b & 1:
		return mul(a, pow(a, b - 1))
	else:
		aux = pow(a, b >> 1)
		return mul(aux, aux)

def inv(a):
	return pow(a, mod - 2)

def rec(g, m):
	if g == 1:
		return 0
	if dp[g] != -1:
		return dp[g]
	fac_len = len(factors[g])
	X = [0] * fac_len
	for i in range(fac_len - 1, -1, -1):
		X[i] = m // factors[g][i]
		for j in range(i + 1, fac_len):
			if factors[g][j] % factors[g][i] == 0:
				X[i] -= X[j]
	res = 1
	for i in range(fac_len - 1):
		prob = mul(X[i], inv_m)
		cur_res = rec(factors[g][i], m)
		new_g_res = mul(prob, cur_res)
		res = add(res, new_g_res)
	a = mul(m, inv(sub(m, X[-1])))
	res = mul(res, a)
	dp[g] = res
	return res
m = int(sys.stdin.read())
inv_m = inv(m)
factors = [[] for _ in range(m + 1)]
for i in range(1, m + 1):
	j = i
	while j <= m:
		factors[j].append(i)
		j += i
dp = [-1] * (m + 1)
res = 0
for i in range(1, m + 1):
	prob = inv_m
	cur_res = add(rec(i, m), 1)
	cur_res = mul(prob, cur_res)
	res = add(res, cur_res)
print(res)


q = 100010

def rop():
	a = [1] * q
	a[0] = 0
	s = [True] * q
	s[0] = s[1] = False
	for i in range(2, q):
		if s[i]:
			o = i
			while o < q:
				s[o] = False
				a[o] *= -1
				o += i
			o = i ** 2
			while o < q:
				a[o] = 0
				o += i ** 2
	return a
d = int(input())
a = rop()

def pro(x):
	return pow(x, 10 ** 9 + 5, 10 ** 9 + 7)
f = 1
for i in range(2, q):
	z = d // i * pro(d)
	f += -a[i] * z * pro(1 - z)
print(f % (10 ** 9 + 7))