Project Euler 307 题解

主要是规避误差。即求

\[\frac{k^n}{n^k} \]

微分一下得到递推式。然后根据斯特林近似（byd 这里还需要 \(1\) 后的第一项。。）

\[k!=\frac{k^k}{e^k}\sqrt{2\pi k}(1+\frac{1}{12}k^{-1}-O(k^{-2})) \]

然后在递推时记录 \(f_k=g_k\times \left(\frac{en}{k}\right)^{c_k}\) 即可转移。

注意代码里 \(n,k\) 是反的。

#include<bits/stdc++.h>
using namespace std;
#define int long long
const int maxn=1e6+5;
// #define double long double
const double pi=acos(-1.0);
double f[maxn],e=0;int cnt[maxn];
double calc(int n,int k){
	f[0]=1,f[1]=n,f[2]=(1.0*(n-1)*f[1]+n*f[0])/2.0;
	double C=e*n/k;
	while(f[2]>C)f[2]/=C,cnt[2]++;
	for(int i=2;i<k;i++){
		double T=(n-(i-1)/2.0)*f[i-1]/(i+1.0);
		T+=(n-i)/(i+1.0)*f[i]*pow(C,cnt[i]-cnt[i-1]);
		f[i+1]=T,cnt[i+1]=cnt[i-1];
		while(f[i+1]>C)f[i+1]/=C,cnt[i+1]++;
		// cout<<cnt[i+1]<<" "<<f[i+1]<<endl;
	}
	double ans=f[k]*pow(C,cnt[k]-k)*sqrt(2*pi*k)*(1+1/(12.0*k));
	return 1.0-ans;
}
signed main(){
	ios::sync_with_stdio(0);
	cin.tie(0);cout.tie(0);
	double z=1;
	for(int i=1;i<=100;i++){
		e+=z;if(i>0)z/=i;
	}
	int n,k;cin>>n>>k;
	printf("%.10lf",calc(n,k));
	return 0;
}