Read problems statements in Mandarin Chinese  and Russian. 
Mike likes to invent new functions. The latest one he has invented is called SuperFunction. Let's consider how it can be calculated:

You are given two integers N and K. SuperFunction of N and K equals to the sum of K'th powers of the positive numbers, which are coprime with N and also not greater than N.

I.e. SuperFuntion of 6 and 3 equals to 1^{3} + 5^{3} = 126.

Mike knows how to calculate SuperFunction in case N and K are rather small integers. Now he wants you to write him a solution, which can calculate SuperFuntion for bigger constraints. As far as SuperFunction could be extremely huge, you are asked to do all the calculation under modulo M.

Help Mike!

------ Input ------ 

The first line contains three integers N, K and M.

------ Output ------ 

The first line of contains the only integer, denoting SuperFuntion of N and K under modulo M.

------ Example ------ 

Input:
5 2 100
Output:
30

------ Explanation ------ 

In the test case N equals to 5, K equals to 2, M equals to 100. The answer is (1^{2} + 2^{2} + 3^{2} + 4^{2}) mod 100 = 30.

------ Scoring ------ 

Subtask 1 (12 points): 1 ≤ N ≤ 100, 1 ≤ K ≤ 2, 1 ≤ M ≤ 40000;

Subtask 2 (28 points): 1 ≤ N ≤ 1 000 000, 1 ≤ K ≤ 2, 1 ≤ M ≤ 10^{9};

Subtask 3 (30 points): 1 ≤ N ≤ 1 000 000 000 000, 1 ≤ K ≤ 3, 1 ≤ M ≤ 10^{9};

Subtask 4 (30 points): 1 ≤ N ≤ 1 000 000 000 000, 1 ≤ K ≤ 10, 1 ≤ M ≤ 10^{18}.

Difficulty: VERY_HARD
Tags: Number theory

Mathematics

Starter code:

Solutions:
import bisect
import copy
import decimal
import fractions
import functools
import heapq
import itertools
import math
import random
import sys
from collections import Counter, deque, defaultdict
from functools import lru_cache, reduce
from heapq import heappush, heappop, heapify, heappushpop, _heappop_max, _heapify_max

def _heappush_max(heap, item):
	heap.append(item)
	heapq._siftdown_max(heap, 0, len(heap) - 1)

def _heappushpop_max(heap, item):
	if heap and item < heap[0]:
		(item, heap[0]) = (heap[0], item)
		heapq._siftup_max(heap, 0)
	return item
from math import gcd as GCD
read = sys.stdin.read
readline = sys.stdin.readline
readlines = sys.stdin.readlines

def Comb(N, K, mod=0):
	if K < 0 or K > N:
		return 0
	K = min(K, N - K)
	s = 1
	if mod:
		for i in range(N, N - K, -1):
			s *= i
			s %= mod
		ss = 1
		for i in range(1, K + 1):
			ss *= i
			ss %= mod
		s *= MOD(mod).Pow(ss, -1)
		s %= mod
	else:
		for i in range(N - K + 1, N + 1):
			s *= i
		for i in range(1, K + 1):
			s //= i
	return s

def Bernoulli_Numbers(N, mod=0):
	bernoulli_numbers = [0] * (N + 1) if mod else [fractions.Fraction(0)] * (N + 1)
	bernoulli_numbers[0] += 1
	if mod:
		MD = MOD(mod)
		MD.Build_Fact(N + 1)
		for i in range(1, N + 1):
			bernoulli_numbers[i] = -MD.Pow(i + 1, -1) * sum((MD.Comb(i + 1, j) * bernoulli_numbers[j] for j in range(i))) % mod
	else:
		for i in range(1, N + 1):
			bernoulli_numbers[i] = -sum((Comb(i + 1, j) * bernoulli_numbers[j] for j in range(i))) / (i + 1)
	return bernoulli_numbers

class Faulhaber:

	def __init__(self, K, mod=0):
		self.K = K
		self.mod = mod
		if self.mod:
			bernoulli_numbers = Bernoulli_Numbers(self.K, self.mod)
			MD = MOD(self.mod)
			MD.Build_Fact(self.K + 1)
			inve = MD.Pow(self.K + 1, -1)
			self.coefficient = [bernoulli_numbers[i] * MD.Comb(self.K + 1, i) * inve % mod for i in range(self.K + 1)]
			for i in range(1, self.K + 1, 2):
				self.coefficient[i] *= -1
				self.coefficient[i] %= mod
		else:
			bernoulli_numbers = Bernoulli_Numbers(self.K)
			self.coefficient = [bernoulli_numbers[i] * Comb(self.K + 1, i) / (K + 1) for i in range(self.K + 1)]
			for i in range(1, self.K + 1, 2):
				self.coefficient[i] *= -1

	def __call__(self, N):
		retu = 0
		N_pow = N
		for i in range(self.K + 1):
			retu += N_pow * self.coefficient[self.K - i]
			N_pow *= N
			if self.mod:
				retu %= self.mod
				N_pow %= self.mod
		return retu

def Extended_Euclid(n, m):
	stack = []
	while m:
		stack.append((n, m))
		(n, m) = (m, n % m)
	if n >= 0:
		(x, y) = (1, 0)
	else:
		(x, y) = (-1, 0)
	for i in range(len(stack) - 1, -1, -1):
		(n, m) = stack[i]
		(x, y) = (y, x - n // m * y)
	return (x, y)

class MOD:

	def __init__(self, p, e=1):
		self.p = p
		self.e = e
		self.mod = self.p ** self.e

	def Pow(self, a, n):
		a %= self.mod
		if n >= 0:
			return pow(a, n, self.mod)
		else:
			assert math.gcd(a, self.mod) == 1
			x = Extended_Euclid(a, self.mod)[0]
			return pow(x, -n, self.mod)

	def Build_Fact(self, N):
		assert N >= 0
		self.factorial = [1]
		self.cnt = [0] * (N + 1)
		for i in range(1, N + 1):
			ii = i
			self.cnt[i] = self.cnt[i - 1]
			while ii % self.p == 0:
				ii //= self.p
				self.cnt[i] += 1
			self.factorial.append(self.factorial[-1] * ii % self.mod)
		self.factorial_inve = [None] * (N + 1)
		self.factorial_inve[-1] = self.Pow(self.factorial[-1], -1)
		for i in range(N - 1, -1, -1):
			ii = i + 1
			while ii % self.p == 0:
				ii //= self.p
			self.factorial_inve[i] = self.factorial_inve[i + 1] * ii % self.mod

	def Fact(self, N):
		return self.factorial[N] * pow(self.p, self.cnt[N], self.mod) % self.mod

	def Fact_Inve(self, N):
		if self.cnt[N]:
			return None
		return self.factorial_inve[N]

	def Comb(self, N, K, divisible_count=False):
		if K < 0 or K > N:
			return 0
		retu = self.factorial[N] * self.factorial_inve[K] * self.factorial_inve[N - K] % self.mod
		cnt = self.cnt[N] - self.cnt[N - K] - self.cnt[K]
		if divisible_count:
			return (retu, cnt)
		else:
			retu *= pow(self.p, cnt, self.mod)
			retu %= self.mod
			return retu
(N, K, M) = map(int, readline().split())
F = Faulhaber(K)
lst = []
NN = N
for p in range(2, 10 ** 6):
	if NN % p == 0:
		lst.append(p)
		while NN % p == 0:
			NN //= p
if NN != 1:
	lst.append(NN)
l = len(lst)
ans = 0
for bit in range(1 << l):
	NN = N
	cnt = 0
	for i in range(l):
		if bit >> i & 1:
			NN //= lst[i]
			cnt += 1
	if cnt % 2:
		ans -= F(NN) * pow(N // NN, K, M)
	else:
		ans += F(NN) * pow(N // NN, K, M)
	ans %= M
print(ans)