Description：

给一个函数： $$ f(x)=\prod_{i=0}^k(x-i^2) $$ $k$ 为满足 $x-k^2>0$ 的最大自然数。

给出 $L,R,P$ ，求有多少个正整数 $x\in[L,R]$ 满足 $P|f(x)$ 。

$1\leqslant L,R\leqslant 10^{18},1\leqslant P\leqslant 10^6$

Solution：

设： $$ f_k(x)\prod_{i=0}^k(x-i^2) $$ 首先由于 $f$ 是一个多项式，因此有： $f(x+P)=f(x)\bmod P$ 。那么我么就可以对于每个 $i\in[0,P)$ 统计一下 $\bmod P=i$ 的 $x$ 中合法的个数，这个可以这样做，我们可以求出一个 $maxk[i]$ 表示 $\bmod P=i$ 的 $x$ 中 $f(x)\bmod P=0$ 的最小的 $k$ ，因为由于 $f(x)$ 单调而且是一个多项式，所以如果 $f_k(x)$ 合法，那么 $k+1$ 也一定是合法的，也就是说 $x>k^2,x\bmod P=i$ 的 $x$ 都合法，因此我们如果可以计算出 $maxk[i]$ 就可以计算 $x\bmod P=i,P|f(x)$ 的 $x$ 个数了，我们可以再设一个 $dp[i][j]$ 表示模 $P$ 为 $i$ 并且 $f(x)$ 是质数 $j$ 在 $P$ 中的幂 $j^q$ 的倍数的 $k$ 最小是几，那么有 $maxk[i]=\max_j(dp[i][j])$ ，那么我们就可以不停找直到找够了 $q$ 个为止，也就是说从小往大枚举 $x-k^2=0\bmod j$ 的 $k$ ，每次消掉 $x-k^2$ 中 $j$ 的幂，会发现这个是一个二次剩余，但是我们可以把它提前预处理出来放进一个 $vector$ 里。

Code：

#include<algorithm>

#include<iostream>

#include<cstdlib>

#include<cstdio>

#include<cmath>

#include<vector> 

#include<cctype>

#include<cstring>

using namespace std;

typedef long long ll;

ll L,R,P;

#define MAXN 1000010

ll ans[MAXN];

vector<ll> g[MAXN];

#define INF 0x3f3f3f3f3f3f3f3f

ll calc(ll n)

{

	if(n == 0)return 0;

	ll res = 0;

	for(int i = 0;i < P;++i)

	{

		if(ans[i] != INF)

		{

			if(ans[i] > n)continue;

			if(n >= i)res += (n - i) / P + 1;

			if(ans[i] >= i)res -= (ans[i] - i) / P + 1;

		}

		//cout << ans[i] << " "; 

	}//cout << endl;

	return res;

}

void solve(ll p,ll k)

{

	for(ll i = 0;i < P;++i)g[i].clear();

	for(ll i = 0;i < P;++i)g[i * i % p].push_back(i);

	for(ll i = 0;i < P;++i)

	{

		vector<ll> v = g[i % p];

		int cnt = k;

		ll x = 0;

		while(cnt > 0)

		{

			for(vector<ll>::iterator it = v.begin();it != v.end();++it)

			{

				ll val = (i - *it * *it % P + P) % P;//cout << val << endl;

				while(cnt && val % p == 0)val /= p,--cnt;

				if(cnt == 0){ans[i] = max(ans[i],(*it + x) * (*it + x));break;}

			}

			if(cnt == k){ans[i] = INF;break;}

			x += P;

		}

	}

	return;

}

int main()

{

	scanf("%lld%lld%lld",&L,&R,&P);

	ll pr = P;

	for(int i = 2;i * i <= P;++i)

	{

		if(P % i == 0)

		{

			int cnt = 0;

			while(pr % i == 0)++cnt,pr /= i;

			solve(i,cnt);

		}

	}

	if(pr != 1)solve(pr,1);

	cout << (calc(R) - calc(L - 1)) << endl;

	return 0;

}