[概率论 DP] BZOJ 3652 大新闻

边看题解边凑公式

#include<cstdio>
#include<cstdlib>
#include<algorithm>
using namespace std;
typedef long long ll;

ll n,len;
double p;
double a[105];

inline double part2()
{
	n++;
	double ans=0;
	for (int i=len-1;i>=0;i--)
	{
		double pi=( (double)(n/(1LL<<(i+1)))*(1LL<<i) + (double)max((n%(1LL<<(i+1)))-(1LL<<i),0LL) )/n;
		ans+=pi*(1-pi)*2*(1LL<<i);
	}
	return ans;
}

inline double part1()
{
	n--;
	double ans=0; ll m=n;
	a[0]=n&1?2.0/(n+1):1.0/(n+1);
	for (int i=1;i<len-1;i++)
		if (n&(1LL<<i))
			a[i]=a[i-1]+(double)(1LL<<i)/(n+1)*(1LL<<i)*2+(double)((1LL<<i)-1)/(n+1)*(1LL<<i);
		else
			a[i]=a[i-1]*2+(double)(1LL<<i)/(n+1)*(1LL<<i);
	for (int i=len-1;i>=0;i--)
	{
		if (n&(1LL<<i)) 
		{	
			if (m&(1LL<<i))
			{
				ans+=(double)((1LL<<(i+1))-1)*(m+1-(1LL<<i))/(n+1);
				m=(1LL<<i)-1;
			}
			ans+=(double)(1LL<<i)*(m+1)/(n+1);
		}
		else if (m&(1LL<<i))
		{
			ans+=(double)(1LL<<i)*(m+1-(1LL<<i))/(n+1);
			m^=(1LL<<i);
			ans+=i-1>=0?a[i-1]:0;
		}
	}
	return ans;
}

int main()
{
	freopen("t.in","r",stdin);
	freopen("t.out","w",stdout);
	scanf("%lld%lf",&n,&p);
	while ((1LL<<len)<=n) 
		len++;
	double ans1=part1();
	double ans2=part2();
	printf("%.6lf\n",ans1*p+ans2*(1-p));
	return 0;
}

UPD:正解数位DP 2017.07.05

#include<cstdio>
#include<cstdlib>
#include<algorithm>
#include<iostream>
using namespace std;
typedef double ld;
typedef long long ll;

const int N=66;

ll n;
int a[N],m;

namespace Work1{
  ld f[N][2][2];
  ld ans=0;
  inline ld Solve(){
    f[m+1][0][0]=1;
    for (int i=m+1;i>1;i--)
      for (int t1=0;t1<2;t1++)
	for (int t2=0;t2<2;t2++)
	  if (f[i][t1][t2]){
	    int xlim=t1?2:a[i-1]+1,ylim=t2?2:a[i-1]+1;
	    for (int j=0;j<xlim;j++){
	      int k=(j^1)<ylim?(j^1):j;
	      ld p;
	      if (t1)
		p=1.0/xlim;
	      else
		if (xlim==2)
		  p=(ld)(1LL<<(i-2))/((n&((1LL<<(i-1))-1))+1),p=j==0?p:1-p;
		else
		  p=1;
	      ans+=((ll)(j^k)<<(i-2))*(f[i][t1][t2]*p);
	      f[i-1][t1|(j<xlim-1)][t2|(k<ylim-1)]+=f[i][t1][t2]*p;
	    }
	  }
    return ans;
  }
}

namespace Work2{
  ld f[N][2][2];
  ld ans=0;
  inline ld Solve(){
    f[m+1][0][0]=1;
    for (int i=m+1;i>1;i--)
      for (int t1=0;t1<2;t1++)
	for (int t2=0;t2<2;t2++)
	  if (f[i][t1][t2]){
	    int xlim=t1?2:a[i-1]+1,ylim=t2?2:a[i-1]+1;
	    for (int j=0;j<xlim;j++)
	      for (int k=0;k<ylim;k++){
		ld p1,p2;
		if (t1)
		  p1=1.0/xlim;
		else
		  if (xlim==2)
		    p1=(ld)(1LL<<(i-2))/((n&((1LL<<(i-1))-1))+1),p1=j==0?p1:1-p1;
		  else
		    p1=1;
		if (t2)
		  p2=1.0/ylim;
		else
		  if (ylim==2)
		    p2=(ld)(1LL<<(i-2))/((n&((1LL<<(i-1))-1))+1),p2=k==0?p2:1-p2;
		  else
		    p2=1;
		ans+=((ll)(j^k)<<(i-2))*(f[i][t1][t2]*p1*p2);
		f[i-1][t1|(j<xlim-1)][t2|(k<ylim-1)]+=f[i][t1][t2]*p1*p2;
	      }
	  }
    return ans;
  }
}

int main(){
  ld p,ans1,ans2;
  freopen("t.in","r",stdin);
  freopen("t.out","w",stdout);
  cin>>n>>p; n--;  
  ll t=n;
  while (t)
    a[++m]=t&1,t>>=1;
  ans1=Work1::Solve();
  ans2=Work2::Solve();
  printf("%.6lf\n",(double)(ans1*p+ans2*(1-p)));
  return 0;
}