BZOJ1894 : Srm444 avoidfour

首先只有质数个$4$且个数不超过$10$的限制条件才有用,

也就是长度不能为$44,444,44444,4444444$的倍数。

考虑容斥,计算长度必须是它们$lcm$的倍数,且没有连续$4$个$4$的方案数。

将DP转移方程用矩阵表示,则长度为$L$的方案数为$G^{L-1}\times V$。

所以对于一个$lcm$,答案为$(\sum_{k=0}^{\lfloor\frac{n}{lcm}\rfloor-1} (G^{lcm})^k)\times G^{lcm-1}\times V$,分治求解即可。

 

#include<cstdio>
#define rep(i) for(int i=0;i<4;i++)
typedef long long ll;
const int P=1000000007;
struct mat{
  int v[4][4];
  mat(){rep(i)rep(j)v[i][j]=0;}
  mat operator+(const mat&b){
    mat c;
    rep(i)rep(j)c.v[i][j]=(v[i][j]+b.v[i][j])%P;
    return c;
  }
  mat operator*(const mat&b){
    mat c;
    rep(i)rep(j)rep(k)c.v[i][j]=(1LL*v[i][k]*b.v[k][j]+c.v[i][j])%P;
    return c;
  }
}G,A,B,I;
ll a[4]={44,444,44444,4444444},n,t;int S,ans,ret;
ll gcd(ll a,ll b){return b?gcd(b,a%b):a;}
ll lcm(ll a,ll b){return a/gcd(a,b)*b;}
mat pow(mat a,ll b){mat t=I;for(;b;b>>=1,a=a*a)if(b&1)t=t*a;return t;}
struct PI{mat x,y;PI(){}PI(mat _x,mat _y){x=_x,y=_y;}};
PI cal(const mat&a,ll b){
  if(b==1)return PI(I,a);
  if(b&1){
    PI c=cal(a,b-1);
    return PI(c.x*a+I,c.y*a);
  }
  PI c=cal(a,b>>1);
  return PI(c.x*(c.y+I),c.y*c.y);
}
int main(){
  rep(i)G.v[0][i]=9,I.v[i][i]=1;
  G.v[1][0]=G.v[2][1]=G.v[3][2]=1;
  A.v[0][0]=8,A.v[1][0]=1;
  scanf("%lld",&n);
  for(S=0;S<16;S++){
    t=1;
    rep(i)if(S>>i&1)t=lcm(t,a[i]);
    if(t>n)continue;
    B=cal(pow(G,t),n/t).x*pow(G,t-1)*A;
    ret=0;
    rep(i)ret=(ret+B.v[i][0])%P;
    if(__builtin_popcount(S)&1)ans=(ans+P-ret)%P;else ans=(ans+ret)%P;
  }
  return printf("%d",ans),0;
}

  

转载于:https://www.cnblogs.com/clrs97/p/5648071.html

评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

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

抵扣说明:

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

余额充值