How many infinite sequences a_1, a_2, ... consisting of {{1, ... ,n}} satisfy the following conditions?
 - The n-th and subsequent elements are all equal. That is, if n \leq i,j, a_i = a_j.
 - For every integer i, the a_i elements immediately following the i-th element are all equal. That is, if i < j < k\leq i+a_i, a_j = a_k.
Find the count modulo 10^9+7.

-----Constraints-----
 - 1 \leq n \leq 10^6

-----Input-----
Input is given from Standard Input in the following format:
n

-----Output-----
Print how many sequences satisfy the conditions, modulo 10^9+7.

-----Sample Input-----
2

-----Sample Output-----
4

The four sequences that satisfy the conditions are:
 - 1, 1, 1, ...
 - 1, 2, 2, ...
 - 2, 1, 1, ...
 - 2, 2, 2, ...

Difficulty: competition

Starter code:

Solutions:
# dp[i] := i 項目以降を見たときの場合の数
#       = dp[i+1] + (Σ_{3<=d<=N-1} dp[i+d]) + (n-1)**2

# i 項目と i+1 項目について、
#   (1) i 項目が 1 の場合
#   (2) i 項目が 1 で無く、i+1 項目が 1 の場合
#   (3) i 項目も i+1 項目も 1 で無い場合
# で場合分けしてそれぞれを足している
# これを累積和で高速化する

# 初期化がややこしい

N = int(input())
mod = 10**9+7
dp = [N] * (N+1)
dp[N-1] = N
dp[N-2] = N*N
c = N * (N+1) + N-1
const = (N-1)**2
for i in range(N-3, -1, -1):
    ans = c - dp[i+2] + const
    dp[i] = ans
    c = (c + ans - 1) % mod
print((dp[0]%mod))


n = int(input())
mod = 10**9+7

dp = [0] * n
dp[0] = 1
tot = 1
for i in range(1,n):
    if(i==1):
        dp[i] = tot
    else:
        dp[i] = tot - dp[i-2]
    dp[i] %= mod
    tot += dp[i]
    tot %= mod

ans = 0
for i in range(n-1):
    ans += dp[i] * ((n-1)**2 + min(n-1, i+2))
    ans %= mod

ans += dp[-1]*n
ans %= mod
print(ans)

n=int(input())
dp=[0 for i in range(n+1)]
dp[0]=1
x=(n-1)**2
dp[1]=n
pref=dp[1]
mod=10**9+7
x%=mod
for i in range(2,n+1):
    dp[i]+=x+pref-dp[i-2]+(n-i+2)
    pref+=dp[i]
    pref%=mod
    dp[i]%=mod
        
print((dp[-1]))
    
    
    


# dp[i] := i 項目以降を見たときの場合の数
#       = dp[i+1] + (Σ_{3<=d<=N-1} dp[i+d]) + (n-1)**2

# i 項目と i+1 項目について、
#   (1) i 項目が 1 の場合
#   (2) i 項目が 1 で無く、i+1 項目が 1 の場合
#   (3) i 項目も i+1 項目も 1 で無い場合
# で場合分けしてそれぞれを足している
# これを累積和で高速化する

# 初期化がややこしい

N = int(input())
mod = 10**9+7
dp = [0] * N
dp[N-1] = N
dp[N-2] = N*N
c = N * (N+1) + N-1
const = (N-1)**2
for i in range(N-3, -1, -1):
    ans = c - dp[i+2] + const
    dp[i] = ans
    c = (c + ans - 1) % mod
print((dp[0]%mod))


import sys
def I(): return int(sys.stdin.readline().rstrip())


N = I()
mod = 10**9+7
A = [0]*(N+1)
SA = [0]*(N+1)

A[1] = N
SA[1] = N
if N >= 2:
    A[2] = N**2
    SA[2] = N+N**2

for i in range(3,N+1):
    A[i] = A[i-1]+SA[i-3]+(N-i+2)+(N-1)**2
    A[i] %= mod
    SA[i] = SA[i-1]+A[i]
    SA[i] %= mod

print((A[-1]))


import math
#import numpy as np
import queue
from collections import deque,defaultdict
import heapq as hpq
from sys import stdin,setrecursionlimit
#from scipy.sparse.csgraph import dijkstra
#from scipy.sparse import csr_matrix
ipt = stdin.readline
setrecursionlimit(10**7)

def main():
    n = int(ipt())
    mod = 10**9+7
    a = [n,n*n%mod]
    sa = [0,0]
    for i in range(2,n):
        a.append((a[-1]+sa[-1]+(n-i+1)*n%mod)%mod)
        sa.append((sa[-1]+a[i-2]+n-1)%mod)

    print((a[-1]))



    return

def __starting_point():
    main()

__starting_point()

mod = 1000000007
eps = 10**-9


def main():
    import sys
    input = sys.stdin.readline

    N = int(input())
    ans = [0] * (N+1)
    ans[N] = N
    ans[N-1] = N*N
    S = 0
    for i in range(N-2, 0, -1):
        ans[i] = (ans[i+1] + S + N*(i+1) + (N-1) * (N-i-2))%mod
        S = (S + ans[i+2])%mod
    print((ans[1]))
    #print(ans)


def __starting_point():
    main()

__starting_point()

N=int(input())
mod=10**9+7
DP=[0]*(N*2+5)
X=[0]*(N*2+5)
Y=[0]*(N*2+5)
C=0
Y[0]=-1
for i in range(N):
  C=(C+X[i])%mod
  DP[i]=(C-Y[i])%mod
  Y[i+2]=DP[i]
  X[i+1]=(X[i+1]+DP[i])%mod
  X[i+N+2]=(X[i+N+2]-DP[i])%mod
for i in range(N,N*2+2):
  C=(C+X[i])%mod
  DP[i]=(C-Y[i])%mod
P=sum(DP[N:])
for i in range(N-1):
  P=(P+DP[i]*(N-1)*(N-1))%mod
print(P)

import sys
MOD = 10**9 + 7
N = int(input())
if N == 1:
  print((1)); return;
if N == 2:
  print((4)); return;

DP = [0 for _ in range(N)]
DP[0] = 1
DP[1] = 1
DP_sum = 2
for i in range(2,N):
  DP[i] =  DP_sum - DP[i-2] # DP[0] + DP[1] + ... + DP[i-1] - DP[i-2]
  DP[i] %= MOD
  DP_sum += DP[i]
  DP_sum %= MOD
#print(DP)

ans = 0
fac = (N-1)**2 % MOD
for i in range(N-1):
  ans += DP[i]*(fac+i+2)
  if i >= N-2:
    ans -= DP[i]
  ans %= MOD
ans += DP[N-1] * N
ans %= MOD
  
print((ans % MOD))


n=int(input())
dp=[0 for i in range(n+1)]
dp[0]=1
x=(n-1)**2
dp[1]=n
pref=dp[1]
mod=10**9+7
x%=mod
for i in range(2,n+1):
    dp[i]+=x+pref-dp[i-2]+(n-i+2)
    pref+=dp[i]
    pref%=mod
    dp[i]%=mod
        
print((dp[-1]))


mod = 10**9+7

N = int(input())

X = [1, 1, 1]
Y = [1, 2, 3]
for i in range(3, N+2):
    x = (X[-1] + Y[i-3]) % mod
    X.append(x)
    Y.append((Y[-1]+x)%mod)

ans = ((Y[N-2] * (N-2) + (Y[N-1])) * (N-1) + (N-1) * Y[N-2] + 1) % mod
print(ans)


# from itertools import product

# A = product(range(1, N+1), repeat=N)
# P = []
# for B in A:
#     ok = True
#     B = list(B) + [B[-1]]*(N+1)
#     for i in range(N):
#         b = B[i]
#         if len(set(B[i+1:i+1+b])) != 1:
#             ok = False
#             break
#     if ok:
#         P.append(B)
# print(len(P))
# for B in P:
#     print(B)


n = int(input())
if n <= 3:
  print([n,n**2,n*(2*n-1)][n-1])
  return
mod = 10**9+7
dp = [0]*(n+1)
dp[1] = n
dp[2] = n**2
dp[3] = n*(2*n-1)
acc = sum(dp)
for i in range(4,n+1):
  dp[i] = acc-dp[i-2]+(n-1)**2+n-i+2
  dp[i] %= mod
  acc = (acc+dp[i])%mod
print(dp[n])

P = 10 ** 9 + 7
N = int(input())
k = (N-1) ** 2 % P + 2
a, b, c = 1, 0, 0
ans = 0
for i in range(N-1):
    ans += a * (k + min(i, N-3))
    ans %= P
    a, b, c = a + b, b + c, a
    if a >= P: a -= P
    if b >= P: b -= P

ans += N * a
ans %= P
print(ans)

import sys
input = sys.stdin.readline

MOD = 10**9 + 7
N = int(input())

dp = [0] * (N+10)
dp_cum = [0] * (N+10)

dp[1] = N-1; dp_cum[1] = N-1
dp[2] = N-1; dp_cum[2] = 2*(N-1)
for n in range(3,N+1):
    dp[n] = dp[n-1] + dp_cum[n-3]
    dp_cum[n] = (dp_cum[n-1] + dp[n]) % MOD

answer = sum(dp[1:N])*N + dp[N] + 1
answer %= MOD
print(answer)

def main():
    mod = 10**9+7
    n = int(input())
    a, b, c, p = 1, 1, n, n-1
    for i in range(n-1):
        p = (p+a-1) % mod
        a, b, c = b, c, ((n-1)*(n-1)+p+c) % mod
    print(c)


main()


n = int(input())
dp = [0] * (n+1)
mod = 10**9+7
dp[0] = n
dp[1] = n*n
SUM = n+n*n
for i in range(2,n):
  dp[i] = (n-1)*(n-1) + n-i+1+SUM-dp[i-2]
  dp[i] %= mod
  SUM += dp[i]
  SUM %= mod
print(dp[n-1])

n=int(input())-1;a,b,c,p=1,1,n+1,n
for _ in[0]*n:p+=a-1;a,b,c=b,c,(n*n+p+c)%(10**9+7)
print(c)

n=int(input());a,b,c,p=1,1,n,n-1
for i in range(n-1):
 p+=a-1;a,b,c=b,c,((n-1)**2+p+c)%(10**9+7)
print(c)

n=int(input())
mod=10**9+7

dp=[0 for i in range(0,10**6+1)]
S=[0 for i in range(0,10**6+1)]

dp[1]=1
dp[2]=1
S[1]=1
S[2]=2
for i in range(3,10**6+1):
    dp[i]=(S[i-1]-dp[i-2]+1)%mod
    S[i]=(S[i-1]+dp[i])%mod

if n==1:
    print(1)
elif n==2:
    print(4)
else:
    ans=(S[n-2]*(n-1)**2+(n-1)**2+(n-1)*S[n-1]+n)%mod
    print(ans)

n=int(input())
if n==1:
  print((1))
  return
mod=10**9+7
dp=[0]*(n+1)
S=[0]*(n+1)
dp[0]=1
dp[1]=1
dp[2]=1
S[0]=1
for i in range(n):
  if i+1>=3:
    dp[i+1]=2*dp[i]-dp[i-1]+dp[i-2]
  dp[i+1]%=mod
  S[i+1]=S[i]+(n-1)*dp[i]
  S[i+1]%=mod
ans=S[n]
for i in range(1,n):
  ans+=((n-1)**2)*dp[i-1]
  ans%=mod
print(ans)
#print(S)
#print(dp)


import sys

mod = 10**9 + 7

def solve():
    n = int(input())

    if n == 1:
        print(1)
        return

    dp = [0] * (n + 1)
    dp[1] = n % mod
    dp[2] = n**2 % mod

    acc = [0] * (n + 1)
    acc[1] = dp[1] % mod
    acc[2] = (dp[1] + dp[2]) % mod

    for k in range(3, n + 1):
        dp[k] = (dp[k - 1] + (n - 1)**2 + acc[k - 3] + n - k + 2) % mod
        acc[k] = (acc[k - 1] + dp[k]) % mod

    ans = dp[n]
    print(ans)

def debug(x, table):
    for name, val in table.items():
        if x is val:
            print('DEBUG:{} -> {}'.format(name, val), file=sys.stderr)
            return None

def __starting_point():
    solve()
__starting_point()

n = int(input())

if n == 1:
    print(1)
else:
    res_v = [0]*(n + 1)
    res_v_cumsum = [0]*(n + 1)
    res_v[0] = 0
    res_v[1] = 1
    res_v[2] = 1
    res_v_cumsum[0] = 0
    res_v_cumsum[1] = 1
    res_v_cumsum[2] = 2

    M = 1000000007

    for k in range(3, n):
        res_v[k] = (1 + res_v_cumsum[k-1] - res_v[k-2]) % M
        res_v_cumsum[k] = (res_v_cumsum[k-1] + res_v[k]) % M
    #print(res_v)
    print((((res_v_cumsum[n-2] * (((n-1) * (n-1)) % M)) % M) + ((res_v_cumsum[n-1] * (n-1)) % M) + n + (n-1)*(n-1)%M) % M)

MOD = 10 ** 9 + 7
n = int(input())

"""
dp = [-1] * (n + 1)
dp[0] = 1
def cnt(num):
    ans = 0
    if num < 0:
        return 1
    if dp[num] != -1:
        return dp[num]
    if num == 1:
        dp[num] = n
        return n
      
    for i in range(n):
        if i == 0:
            ans += cnt(num - 1)
            ans %= MOD
        else:
            if num - 2 - i < 0:
                ans += 1
            else:
                ans += cnt(num - 2 - i)
            ans += n - 1
            ans %= MOD
    dp[num] = ans
    return ans

dp[0] = 1
print(cnt(n))
"""

dq = [0] * (n + 1)
sum_ = [0] * (n + 1)
dq[0] = 1
sum_[0] = 1
for i in range(n):
    if i == 0:
        dq[i + 1] = n
        sum_[i + 1] = sum_[i] + dq[i + 1]
        continue
    """
    for j in range(n):
        if j == 0:
            dq[i + 1] += dq[i]
        else:
            if i - 1 - j < 0:
                dq[i + 1] += 1
            else:
                dq[i + 1] += dq[i - 1 - j]
            dq[i + 1] += n - 1
    """
    dq[i + 1] += dq[i]
    dq[i + 1] += (n - 1) * (n - 1)
    dq[i + 1] += sum_[i - 2]
    dq[i + 1] += abs(i - n)
    dq[i + 1] %= MOD
    sum_[i + 1] = sum_[i] + dq[i + 1]
    sum_[i + 1] %= MOD

print((dq[-1]))


import sys

n = int(input())
mod = 10 ** 9 + 7

if n == 1:
    print((1))
    return

lis = [0] * n
lis[0] = 1
bsum = 0
ans = min(n-1,2) + (n-1) ** 2



for i in range(n-1):

    i += 1

    lis[i] += lis[i-1]

    if i - 3 >= 0:
        bsum += lis[i-3]
        bsum %= mod

    lis[i] += bsum
    lis[i] %= mod

    if i != n-1:

        ans += lis[i] * ((n-1) ** 2)
        ans += lis[i] * min(n-1 , i+2)
        ans %= mod

    

ans += lis[-1] * n
ans %= mod

#print (lis)
print (ans)


import sys
sys.setrecursionlimit(10**7)
def main1(n):
  # DP メモ化再帰 O(N^2)
  mod=10**9+7
  ret=0
  memo={}
  memo[n-1]=n
  #memo[i]:idx=iに何を書き込むか検討しているとき、i以降の数列の種類数
  def func(i):
    if i in memo:return memo[i]
    if i>n-1:return 1
    ret=0
    # iに1を書き込む
    ret+=func(i+1)
    for k in range(2,n+1):
      # iにkを書き込み、その次に1を書き込む。
      ret+=func(i+k+1)
      # iにkを書き込み、その次に2以上を書き込む。
      ret+=n-1
      ret%=mod
    memo[i]=ret
    return ret
  ret=func(0)
  for k in memo:print((k,memo[k]))
  return ret

# 累積和で高速化できる。memo[i]をiの大きい順に数える
def main2(n):
  # DP O(N)
  mod=10**9+7
  # dp[i]:=idx=iに何を書き込むか検討しているとき、i以降の数列の種類数
  dp=[1]*(2*n)
  dp[n-1]=n
  dp[n-2]=n**2
  now=n-1
  for i in range(n-3,-1,-1):
    # 1を書き込む
    dp[i]=dp[i+1]
    now-=dp[i+n+1]
    now+=dp[i+3]
    dp[i]+=now
    #for k in range(2,n+1):# i+3~i+n+1の和
    #  dp[i]+=dp[i+k+1]
    #  dp[i]%=mod
    dp[i]+=(n-1)**2
    dp[i]%=mod
  return dp[0]
def __starting_point():
  n=int(input())
  ret2=main2(n)
  print(ret2)


__starting_point()

# -*- coding: utf-8 -*-
mod = 10**9 + 7

n = int(input())

a = [0]*n
x = n - 1
for i in range(n - 1, - 1, -1):
    if i == n - 1:
        a[i] = n
    elif i == n - 2:
        a[i] = n*n
    else:
        x -= 1
        x += a[i + 3] if i + 3 < n else 1
        x %= mod
        a[i] += a[i + 1] + (n - 1)*(n - 1) + x
    a[i] %= mod
print(a[0])


import sys
readline = sys.stdin.readline
MOD = 10**9+7
N = int(readline())

dp1 = [0]*(N+3)
Dp1 = [0]*(N+3)
dp2 = [0]*(N+3)
Dp2 = [0]*(N+3)
dp1[0] = 1
Dp1[0] = 1
dp2[0] = N-1
Dp2[0] = N-1
dp1[1] = N-1
dp2[1] = (N-1)**2
dp1[2] = N-1
dp2[2] = (N-1)**2
Dp1[1] = Dp1[0] + dp1[1]
Dp2[1] = Dp2[0] + dp2[1]
Dp1[2] = Dp1[1] + dp1[2]
Dp2[2] = Dp2[1] + dp2[2]
for i in range(3, N+1):
    dp1[i] = (dp1[i-1] + Dp1[i-3] - 1) % MOD
    dp2[i] = (dp2[i-1] + Dp2[i-3]) % MOD
    Dp1[i] = (Dp1[i-1] + dp1[i]) % MOD
    Dp2[i] = (Dp2[i-1] + dp2[i]) % MOD
print(((Dp1[N-1] + Dp2[N-1]) % MOD))                 


MOD = 10 ** 9 + 7

n = int(input())

if n == 1:
    print(1)
    return

dp = [0] * (n)
dp[0] = 1
dp[1] = 1
dpsum = 2
for l in range(2, n):
    dp[l] = dpsum
    dp[l] -= dp[l-2]
    dp[l] %= MOD
    dpsum += dp[l]
    dpsum %= MOD

ans = 1 + (dpsum - dp[-1]) * (n-1) # end with "1..."
ans %= MOD
for l in range(n-1):
    ans += dp[l] * (n - 1) ** 2
    ans %= MOD

ans += dp[n-1] * (n - 1)
ans %= MOD

print(ans)

# coding: utf-8
# Your code here!
import sys
read = sys.stdin.read
readline = sys.stdin.readline

def get(i):
    #if i==-n: return 0
    if i <= 0: return 1
    else: return a[i]


n, = list(map(int,readline().split()))

MOD = 10**9+7

a = [1,n,n*n]
for i in range(3,n+1):
    res = -get(i-n-2) + get(i-3) - get(i-2) + 2*get(i-1)
    a.append(res%MOD)

#print(a)
print((a[-1]))













MOD = 10**9 + 7
n = int(input())
if n == 1:
	print(1)
else:
	dp = [0 for _ in range(n+1)]
	dp[0], dp[1], dp[2] = 1, 1, 1
	cum = [0 for _ in range(n+1)]
	cum[0], cum[1], cum[2] = 1, 2, 3
	for i in range(3, n+1):
		dp[i] = (dp[i-1] + cum[i-3]) % MOD
		cum[i] = (cum[i-1] + dp[i]) % MOD
	ans = cum[n-2] * (n-1) * (n-1) + dp[n-1] * (n-1) + cum[n-2] * (n-1) + 1
	print(ans%MOD)