PE556

考虑推广sum(i in Z){mu^2(i)}的做法.

#include"roundCount.cpp"
#include<cstdio>
#include<vector>
#include<cstring>
#include<algorithm>
#include<primesieve.hpp>
namespace GPG{
	typedef long long ll;
	typedef bool B;
	namespace Gauss{
	#define opr(x,op) inline x operator op (
	#define ret(...) ){return ({__VA_ARGS__;});}
	#define AL2(x,y,z) x y,x z
	#define vopr(x,op,y,z,...) opr(x,op) AL2(x,y,z) ret(__VA_ARGS__)
	#define lbd(x,n) inline x n (
		struct GI{
			int a,b; // a+bi
			GI():a(0),b(0){}
			GI(int a,int b=0):a(a),b(b){}
		};
		 lbd(ll,norm)GI a ret(a.a*a.a+a.b*a.b)
		vopr(GI,+,a,b,GI(a.a+b.a,a.b+b.b))
		vopr(GI,-,a,b,GI(a.a-b.a,a.b-b.b))
		vopr(GI,*,a,b,GI(a.a*b.a-a.b*b.b,a.a*b.b+a.b*b.a))
		 opr(B ,<)AL2(GI,a,b)ret(norm(a)<norm(b))
	
	} using namespace Gauss;
	#define maxn 10000011
	using namespace std;
	char v[maxn];
	vector<int> prs;
	int ppp[10000010];
	int nextPP[10001001];
	int main(){
		primesieve::generate_primes(10000000,&prs);
		for(int i=0,_=prs.size();i<_;++i) v[prs[i]]=1;
		for(int i=1,j=3162;i*i<=10000000;++i){
			while(j && i*i+j*j>10000000)--j;
			if(v[i]&&((i&3)==3))ppp[i*i]++;
			for(int k=1;k<=j;++k) if(v[i*i+k*k]) ppp[i*i+k*k]++;
		}
		for(int i=1,last=0;i<=10000000;++i){
			if(ppp[i]) nextPP[last]=i,last=i;
		}
		return 0;
	}
} using namespace GPG;
namespace Combine{
	int C[22][22];
	inline void init_combine(){
		C[0][0]=1;
		for(int i=1;i<=20;++i){
			C[i][0]=1;
			for(int j=1;j<=i;++j) C[i][j]=C[i-1][j]+C[i-1][j-1];
		}
	}
} using namespace Combine;
#define threshold 100000000000000ll
#define thresholdsqrt 10000000ll

/*
#define threshold 10000ll
#define thresholdsqrt 100ll
*/
inline ll dfs(int now,int nowexp,ll mul,ll tot,int mu){
//	printf("%d %d %-4lld %-4lld %-3d\n",now,nowexp,mul,tot,mu);
	ll ans=nowexp?tot*mu:0;
	if(ans){
		ans*=RC(threshold/mul/mul);
//		printf("+ RC(%lld)=%lld\n",threshold/mul/mul,RC(threshold/mul/mul));
	}
	int v=nextPP[now];
	if(!v || v*mul>thresholdsqrt) return ans;
	ans+=dfs(v,0,mul,tot,mu);
	for(int i=1;i<=ppp[v];++i){
		mu=-mu;
		if(double(mul)*v>double(thresholdsqrt))break;
		mul*=v;
		if(mul>thresholdsqrt)break;
		ans+=dfs(v,i,mul,tot*C[ppp[v]][i],mu);
	}
	return ans;
}
int main(){
	RCCC::init();
	GPG::main();
	Combine::init_combine();
	printf("%lld\n",dfs(0,0,1,1,1)+RC(threshold));
	return 0;
}