[期望 DP || 高斯消元 KMP] BZOJ 3213 [Zjoi2013]抛硬币

21 篇文章 0 订阅
7 篇文章 0 订阅

这个其实也不复杂 

先kmp 可以发现 每个点的状态会从转移到 i+1 和 next[i] 不妨设为f

然后列出方程 直接就可以上高斯消元 大概80?

这个东西其实可以DP

E(i+1)=(E(i)-1-p[t^1]*E(f))/p[t]

E(i)=k[i]*E(0)-b[i]

k[i+1]*E(i+1)-b[i+1]=(k[i]*E(0)-b[i]-1-p[t^1]*(k[f]*E(0)-b[f]))/p[t]
k[i+1]=(k[i]-p[t^1]*k[F[i][t^1]])/p[t]
b[i+1]=(b[i]+1-p[t^1]*b[F[i][t^1]])/p[t]

然后最后用E(n)==0就解出来了

只是为什么k b这么设呢 

可以发现 k b 这样的话已知非负 

甚至可以发现 k 一直为1 因为 k0=1 k1=1 ... 然后一直归纳下去就好了

这样就只要dp b数组就行了 


#include<cstdio>
#include<cstdlib>
#include<algorithm>
#include<cstring>
using namespace std;

inline char nc(){
  static char buf[100000],*p1=buf,*p2=buf;
  if (p1==p2) { p2=(p1=buf)+fread(buf,1,100000,stdin); if (p1==p2) return EOF; }
  return *p1++;
}

inline void read(int &x){
  char c=nc(),b=1;
  for (;!(c>='0' && c<='9');c=nc()) if (c=='-') b=-1;
  for (x=0;c>='0' && c<='9';x=x*10+c-'0',c=nc()); x*=b;
}

inline int read(char *s){
  char c=nc(); int len=0;
  for (;!(c=='H' || c=='T');c=nc());
  for (;c=='H' || c=='T';s[++len]=c,c=nc()); s[++len]=0; return len-1;
}

const int con=100000000;
class Int{
public:
  long long a[305];
  Int() { }
  Int(int x){memset(a,0,sizeof(a));while (x){a[++a[0]]=x%con;x=x/con;}}
  void print(){
    if (a[0]==0||(a[0]==1&&a[1]==0)){printf("0");return;}
    printf("%lld",a[a[0]]); for (int i=a[0]-1;i;i--) printf("%08lld",a[i]);
  }
  Int operator +(const Int &X){
    Int c;memset(c.a,0,sizeof(c.a));
    for (int i=1;i<=a[0]||i<=X.a[0];i++) c.a[i]=c.a[i]+a[i]+X.a[i],c.a[i+1]+=c.a[i]/con,c.a[i]%=con;
    c.a[0]=max(a[0],X.a[0]); if (c.a[c.a[0]+1])c.a[0]++; return c;
  }
  Int operator -(const Int &X){
    Int c;memcpy(c.a,a,sizeof(c.a));
    for (int i=1;i<=a[0];i++){ c.a[i]=c.a[i]-X.a[i]; if (c.a[i]<0) c.a[i+1]--,c.a[i]+=con; }
    while (c.a[0]&&!c.a[c.a[0]]) c.a[0]--; return c;
  }
  Int operator *(const Int &X){
    Int c;memset(c.a,0,sizeof(c.a));
    for (int i=1;i<=a[0];i++)
      for (int j=1;j<=X.a[0];j++)
	c.a[i+j-1]+=a[i]*X.a[j],c.a[i+j]+=c.a[i+j-1]/con,c.a[i+j-1]%=con;
    c.a[0]=max(a[0]+X.a[0]-1,0ll); if (c.a[a[0]+X.a[0]]>0)c.a[0]++; return c;
  }
  Int operator *(int num){
    Int c;memset(c.a,0,sizeof(c.a));
    for (int i=1;i<=a[0];i++){ c.a[i]+=a[i]*num; if (c.a[i]>=con) c.a[i+1]+=c.a[i]/con,c.a[i]%=con; }
    c.a[0]=a[0]; if (c.a[c.a[0]+1]>0)c.a[0]++; return c;
  }
  Int operator /(int num){
    Int c;memset(c.a,0,sizeof(c.a));
    long long x=0; for (int i=a[0];i;i--) x=x*con+a[i],c.a[i]=x/num,x=x%num;
    c.a[0]=a[0]; if (c.a[0]&&!c.a[c.a[0]])c.a[0]--; return c;
  }
  long long operator %(int num){
    Int c;memset(c.a,0,sizeof(c.a));
    long long x=0; for (int i=a[0];i;i--) x=x*con+a[i],c.a[i]=x/num,x=x%num;
    c.a[0]=a[0]; if (c.a[0]&&!c.a[c.a[0]])c.a[0]--; return x;
  }
};

const int N=1005;

struct frac{
  Int a,b; // a/b
  void doit(){ for (int i=2;i<=100;i++) while (b%i==0 && a%i==0) a=a/i,b=b/i; }
  frac() {}
  frac(int a,int b):a(a),b(b){ doit(); }
  frac(Int ia,Int ib){ a=ia; b=ib; }
  void print(){ a.print(),putchar('/'),b.print(),putchar('\n'); }
  friend frac operator + (frac A,int p){ return frac(A.a+A.b*p,A.b); }
  friend frac operator + (frac A,frac B){ return frac(A.a*B.b+A.b*B.a,A.b*B.b); }
  friend frac operator - (frac A,frac B){ return frac(A.a*B.b-A.b*B.a,A.b*B.b); }
  friend frac operator * (frac A,frac B){ return frac(A.a*B.a,A.b*B.b); }
  friend frac operator / (frac A,frac B){ return frac(A.a*B.b,A.b*B.a); }
}p[2],k[N],b[N],Ans;

int ia,ib,n;
char S[N];
int next[N],fail[N][2];

inline void KMP(){
  next[1]=0; int k=0;
  for (int i=2;i<=n;i++){
    while (k && S[i]!=S[k+1]) k=next[k];
    if (S[i]==S[k+1]) k++;
    next[i]=k;
  }
}

int main(){
  freopen("coin.in","r",stdin);
  freopen("coin.out","w",stdout);
  read(ia); read(ib);
  n=read(S); KMP();
  p[0]=frac(ia,ib); p[1]=frac(ib-ia,ib);
  int t=S[1]=='T',f; fail[1][t]=2; fail[1][t^1]=1;
  for (int i=2;i<=n;i++)
    t=S[i]=='T',fail[i][t]=i+1,fail[i][t^1]=fail[next[i-1]+1][t^1];
  //k[1]=frac(1,1);
  b[1]=frac(0,1);
  for (int i=1;i<=n;i++){
    t=S[i]=='T',f=fail[i][t^1];
    //(k[i+1]=(k[i]-p[t^1]*k[f])/p[t]).doit();
    (b[i+1]=(b[i]+frac(1,1)-p[t^1]*b[f])/p[t]).doit();
    //k[i+1].print(); b[i+1].print();
  }
  //(Ans=b[n+1]/k[n+1]).doit();
  (Ans=b[n+1]).doit();
  Ans.print();
  return 0;
}


  • 0
    点赞
  • 0
    收藏
    觉得还不错? 一键收藏
  • 0
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值