Problem:
 Score : 1000 points    Problem Statement Takahashi and Aoki are going to together construct a sequence of integers.  First, Takahashi will provide a sequence of integers a , satisfying all of the following conditions:   The length of a is N .  Each element in a is an integer between 1 and K , inclusive.  a is a palindrome , that is, reversing the order of elements in a will result in the same sequence as the original.   Then, Aoki will perform the following operation an arbitrary number of times:   Move the first element in a to the end of a .   How many sequences a can be obtained after this procedure, modulo 10^9+7 ?      Constraints  1≤N≤10^9  1≤K≤10^9         Input The input is given from Standard Input in the following format:  N  K       Output Print the number of the sequences a that can be obtained after the procedure, modulo 10^9+7 .        Sample Input 1 4 2      Sample Output 1 6  The following six sequences can be obtained:   (1, 1, 1, 1)  (1, 1, 2, 2)  (1, 2, 2, 1)  (2, 2, 1, 1)  (2, 1, 1, 2)  (2, 2, 2, 2)        Sample Input 2 1 10      Sample Output 2 10       Sample Input 3 6 3      Sample Output 3 75       Sample Input 4 1000000000 1000000000      Sample Output 4 875699961

Solution (Python):
import sys
input = lambda : sys.stdin.readline().rstrip()
sys.setrecursionlimit(max(1000, 10**9))
write = lambda x: sys.stdout.write(x+"\n")


n,k = list(map(int, input().split()))
M = 10**9+7

def fs(n):
    s = set()
    for i in range(1,int(n**0.5)+2):
        if n%i==0:
            s.add(i)
            s.add(n//i)
    l = sorted(list(s))
    return l
ds = fs(n)
m = len(ds)
nums = [None]*m
for i in range(m):
    val = pow(k, (ds[i]+1)//2, M)
    for j in range(i):
        if ds[i]%ds[j]==0:
            val -= nums[j]
            val %= M
    nums[i] = val%M
ans = 0
inv2 = pow(2, M-2, M)
for i in range(m):
    if ds[i]%2==0:
        ans += ds[i]*nums[i]*inv2
    else:
        ans += ds[i]*nums[i]
    ans %= M
print(ans%M)

Solution (Python):
import sys
def MI(): return map(int,sys.stdin.readline().rstrip().split())

N,K = MI()
mod = 10**9+7

def divisor(n):  
    res = []
    for i in range(1,int(n**.5)+1):
        if n % i == 0:
            res.append(i)
            if i != n//i:
                res.append(n//i)
    return res

A = divisor(N)

d = {}  
for a in A:
    d[a] = divisor(a)

prime = []  
for i in range(2,int(N**.5)+1):
    if N % i == 0:
        prime.append(i)
        while N % i == 0:
            N //= i
if N != 1:
    prime.append(N)

mu = {}  
for a in A:
    b = a
    r = 1
    for p in prime:
        if b % p == 0:
            r *= -1
            if b//p % p == 0:
                r = 0
                break
    mu[a] = r


ans = 0
for a in A:
    for b in d[a]:
        if a % 2 == 0:
            ans += mu[a//b]*pow(K,(b+1)//2,mod)*(a//2)
            ans %= mod
        else:
            ans += mu[a//b]*pow(K,(b+1)//2,mod)*a
            ans %= mod

print(ans)

Solution (Python):
from math import gcd

def divisors(M):
    d=[]
    i=1
    while M>=i**2:
        if M%i==0:
            d.append(i)
            if i**2!=M:
                d.append(M//i)
        i=i+1
    d.sort()
    return d

mod=10**9+7
N,K=map(int,input().split())
div=divisors(N)

res=[0 for i in range(len(div))]
ans=0
if N%2==0:
    for i in range(len(div)):
        d=div[i]
        g=gcd(2*d,N)
        res[i]+=pow(K,g//2,mod)
        res[i]%=mod
        ans+=d*res[i]
        ans%=mod
        for j in range(i+1,len(div)):
            if div[j]%d==0:
                res[j]-=res[i]
                res[j]%=mod

    print(ans%mod)
else:
    for i in range(len(div)):
        d=div[i]
        g=gcd(2*d,N)
        res[i]+=pow(K,(g+1)//2,mod)
        res[i]%=mod
        ans+=d*res[i]
        ans%=mod
        for j in range(i+1,len(div)):
            if div[j]%d==0:
                res[j]-=res[i]
                res[j]%=mod

    print(ans%mod)

Solution (Python):


import sys
read = sys.stdin.read
readline = sys.stdin.readline

n,k = map(int,readline().split())
MOD = 10**9+7

def divisor_list(N): 
    if N == 1: return [1]
    res = []
    for i in range(1,N):
        if i*i >= N: break
        if N%i == 0:
            res.append(i)
            res.append(N//i)
    if i*i == N: res.append(i)
    return sorted(res)



if n&1:
    p = [i for i in divisor_list(n)]
    r = {pi:pow(k,(pi+1)//2,MOD) for pi in p}
    
    for pi in p:
        for pj in p:
            if pj >= pi: break
            if pi%pj==0:
                r[pi] -= r[pj]
    
    ans = 0
    for pi,v in r.items():
        ans += v*(pi)
    ans %= MOD
    print(ans)
    
    

else:
    p = [1] + [i for i in divisor_list(n) if i%2==0 and i!=2]
    
    r = {pi:pow(k,(pi+1)//2,MOD) for pi in p}
    for pi in p:
        for pj in p:
            if pj >= pi: break
            if pi%pj==0:
                r[pi] -= r[pj]

    
    ans = 0
    for pi,v in r.items():
        if pi==1:
            ans += v
        else:
            ans += v*(pi//2)
    ans %= MOD
        
        
    print(ans)

Solution (Python):
def prime_factor(n):
    i = 2
    factors = []
    while i * i <= n:
        if n % i:
            i += 1
        else:
            while n % i == 0:
                n //= i
            factors.append(i)
    if n > 1:
        factors.append(n)
    return factors


def divisors(n):
    divisors = []
    for i in range(1, int(n**0.5)+1):
        if n % i == 0:
            divisors.append(i)
            if i != n // i:
                divisors.append(n//i)
    divisors.sort()
    return divisors


def mebius(n):
    res = {}
    primes = prime_factor(n)
    m = len(primes)
    for i in range(1 << m):
        mu = 1
        d = 1
        for j in range(m):
            if i >> j & 1:
                mu *= -1
                d *= primes[j]
        res[d] = mu
    return res


mod = 10**9+7
N, K = map(int, input().split())
div = divisors(N)

ans = 0

for x in div:
    div2 = divisors(x)
    mu = mebius(x)
    cnt = 0
    for y in mu.keys():
        cnt += pow(K, ((x//y)+1)//2, mod)*mu[y]
    if x % 2 == 0:
        cnt *= pow(2, mod-2, mod)
    ans += x*cnt
    ans %= mod

print(ans)

Solution (Python):
import sys
input = sys.stdin.readline

import numpy as np

MOD = 10**9 + 7

N,K = map(int,input().split())

U=int(10**4.5+10)
div = np.arange(1,U,dtype=np.int64)
div = div[N%div==0]
div = sorted(np.union1d(div,N//div))

A = {d:pow(K,int((d+1)//2),MOD) for d in div}
B = {}
for d in div:
    B[d] = A[d] - sum(B[dd] for dd in div if dd < d and d%dd == 0)

answer = sum(B[d] * (d//2 if d%2 == 0 else d) for d in div)
answer %= MOD
print(answer)

Solution (Python):
N, K = map(int, input().split())
S = []; T = []
for x in range(1, int(N**.5)+1):
    if N % x == 0:
        S.append(x)
        if x*x < N:
            T.append(N//x)
T.reverse()
S += T
M = len(S)
U = []

MOD = 10**9 + 7
ans = 0
v = 0
for i in range(M):
    x = S[i]
    v = pow(K, (x+1)//2, MOD)
    for j in range(i):
        y = S[j]
        if x % y == 0:
            v -= U[j]
    U.append(v % MOD)
    ans = (ans + (v * x if x & 1 else v * (x//2))) % MOD
print(ans)

Solution (Python):
N, K = map(int, input().split())
MOD = 10 ** 9 + 7
i = 1
M = N
div = []
while i * i <= N:
    if N % i == 0:
        div.append(i)
        if N // i != i:
            div.append(N // i)
    i += 1


div = sorted(div)
N1 = len(div)
num = [0] * (N1 + 1)
ans = 0
for i in range(N1):
    num[i] = pow(K, (div[i] + 1) // 2, MOD)
    for j in range(i):
        if div[i] % div[j] == 0:
            num[i] = num[i] - num[j] + MOD
            num[i] %= MOD

    if div[i] % 2 == 0:
        ans += div[i] * num[i] // 2
        ans %= MOD
    else:
        ans += div[i] * num[i]
        ans %= MOD

print(ans)

Solution (Python):
M=10**9+7;D=[];i=1;r=0;N,K=map(int,input().split())
while i*i<=N:
	if N%i==0:D.append(i);i*i<N and D.append(N//i)
	i+=1
D.sort();N=[]
for i in range(len(D)):
	N.append(pow(K,-~D[i]//2,M))
	for j in range(i):
		if D[i]%D[j]==0:N[i]=(N[i]-N[j])%M
	r=(r+(N[i]*D[i]*pow(2,M-2+D[i]%2,M)))%M
print(r)

Solution (Python):
M=10**9+7;D=[];i=1;r=0;N,K=map(int,raw_input().split())
while i*i<=N:
	if N%i==0:D.append(i);i*i<N and D.append(N/i)
	i+=1
D.sort();N=[]
for i in range(len(D)):
	N.append(pow(K,-~D[i]/2,M))
	for j in range(i):
		if D[i]%D[j]==0:N[i]=(N[i]-N[j])%M
	r=(r+(N[i]*D[i]*pow(2,M-2+D[i]%2,M)))%M
print r