LOJ 3090 「BJOI2019」勘破神机——斯特林数+递推式求通项+扩域

题目:https://loj.ac/problem/3090

题解:https://www.luogu.org/blog/rqy/solution-p5320

1.用斯特林数把下降幂化为普通的幂次求和

2.找出通项公式,使得幂次变成二项式,进而将 [ l , r ] 的部分变成等比数列求和

3.模 998244353 下没有 \( \sqrt{5} \) ,所以“扩域”,就是把数表示成 \( a+b*\sqrt{5} \) ;\( \sqrt{3} \) 也同理

注意扩域之后,不满足费马小定理,所以快速幂的指数不能对 ( mod-1 ) 取模!!!

还是不太知道怎么求的通项。为什么是 \( f[n]=A*x_{1}^{n}+B*x_{2}^{n} \) 的形式呢?如果不是二阶怎么推?

UPD:

  设特征根是 x1,x2,...,xk,因为 x^n 是通解,又有线性性(?),所以通项可以写成 \( f(i)=A*x_1^i + B*x_2^i + ... \)

  但是有重根的话就不是这样。 k 重根的系数是次数界为 k 的多项式。这里的次数指的是 i 的几次幂。

  (k重根是针对根而言的,比如一个六次方程,x1=x2=2 , x3=5 , x4=x5=x6=1,那么 2 是2重根,5是单根,1是3重根)

  比如,\( f_i = 2*f_{i-1} - f_{i-2} \),\( f_0 = 1 , f_1 = 2 \) 

  解出特征根是 x1=x2=1 ,那么可以设通项公式为 \( f(i)=(A*i+B)x^i \) ,解得 A=B=1 。

  又如,\( f_i = 4*f_{i-1} - f_{i-2} \) , \( f_0 = 1 , f_1 = 7 \)

  解出特征根是 x1=x2=2 ,设通项是 \( f(i)=(A*i+B)x^i \) ,解得 \( A=\frac{5}{2} , B=1 \) 

  ^  ^

#include<cstdio>
#include<cstring>
#include<algorithm>
#define ll long long
using namespace std;
ll rdn()
{
  ll ret=0;bool fx=1;char ch=getchar();
  while(ch>'9'||ch<'0'){if(ch=='-')fx=0;ch=getchar();}
  while(ch>='0'&&ch<='9')ret=ret*10+ch-'0',ch=getchar();
  return fx?ret:-ret;
}
const int N=505,mod=998244353;
int upt(int x){while(x>=mod)x-=mod;while(x<0)x+=mod;return x;}
int pw(int x,ll k)
{int ret=1;while(k){if(k&1)ret=(ll)ret*x%mod;x=(ll)x*x%mod;k>>=1;}return ret;}

int k,s1[N][N],c[N][N],bs,tlen2,ans;ll tl,tlen;
struct Node{
  int x,y;
  Node(int x=0,int y=0):x(x),y(y) {}
  Node operator+ (const Node &b)const
  { return Node(upt(x+b.x),upt(y+b.y));}
  Node operator- (const Node &b)const
  { return Node(upt(x-b.x),upt(y-b.y));}
  Node operator* (const Node &b)const
  { return Node(((ll)x*b.x+(ll)bs*y%mod*b.y)%mod,((ll)y*b.x+(ll)x*b.y)%mod);}
}A[N],B[N],x1[N],x2[N],one;
Node pw(Node x,ll k)
{ Node ret=Node(1,0);
  while(k){if(k&1)ret=ret*x;x=x*x;k>>=1;}return ret;}
Node Inv(Node u)
{
  int tp=upt(((ll)u.x*u.x-(ll)bs*u.y%mod*u.y)%mod);
  tp=pw(tp,mod-2);
  return Node((ll)u.x*tp%mod,upt(-(ll)u.y*tp%mod));
}
void init(int lx)
{
  s1[0][0]=1;
  for(int i=1;i<=k;i++)
    for(int j=1;j<=i;j++)
      s1[i][j]=(s1[i-1][j-1]+(ll)s1[i-1][j]*(i-1))%mod;
  for(int i=0;i<=k;i++)c[i][0]=1;
  for(int i=1;i<=k;i++)
    for(int j=1;j<=i;j++)
      c[i][j]=upt(c[i-1][j-1]+c[i-1][j]);

  one=Node(1,0);
  if(lx==2)
    {
      int tp=pw(5,mod-2); A[1]=Node(0,tp); B[1]=Node(0,upt(-tp));
      tp=pw(2,mod-2); x1[1]=Node(tp,tp); x2[1]=Node(tp,upt(-tp));
    }
  else
    {
      int tp=pw(6,mod-2); A[1]=Node((ll)3*tp%mod,tp);
      B[1]=Node((ll)3*tp%mod,upt(-tp));
      x1[1]=Node(2,1); x2[1]=Node(2,upt(-1));
    }
  A[0]=one; for(int i=2;i<=k;i++)A[i]=A[i-1]*A[1];
  B[0]=one; for(int i=2;i<=k;i++)B[i]=B[i-1]*B[1];
  x1[0]=one; for(int i=2;i<=k;i++)x1[i]=x1[i-1]*x1[1];
  x2[0]=one; for(int i=2;i<=k;i++)x2[i]=x2[i-1]*x2[1];
}
Node cal(Node x)
{
  if(x.x==1&&x.y==0)return Node(tlen2,0);//tlen2 not tlen!!!
  Node d=Inv(one-x);
  d=d*pw(x,tl)*(one-pw(x,tlen));
  return d;
}
void solve2()
{
  ll l=rdn(),r=rdn(); k=rdn(); bs=5; init(2);
  int iv=pw((r-l+1)%mod,mod-2); l++; r++;
  tl=l; tlen=(r-l+1); tlen2=(r-l+1)%mod;
  for(int i=0,fx=((k&1)?upt(-1):1);i<=k;i++,fx=upt(-fx))
    {
      int tp=0;
      for(int j=0;j<=i;j++)
    {
      Node tmp=x1[j]*x2[i-j];
      Node d=cal(tmp)*A[j]*B[i-j];
      tp=(tp+(ll)c[i][j]*d.x)%mod;
    }
      ans=(ans+(ll)s1[k][i]*fx%mod*tp)%mod;
    }
  ans=(ll)ans*iv%mod;
  int ml=1; for(int i=2;i<=k;i++)ml=(ll)ml*i%mod;
  ans=(ll)ans*pw(ml,mod-2)%mod;
}
void solve3()
{
  ll l=rdn(),r=rdn(); k=rdn(); bs=3; init(3);
  int iv=pw((r-l+1)%mod,mod-2); l=(l+1)>>1; r=r>>1;
  tl=l; tlen=(r-l+1); tlen2=(r-l+1)%mod;
  for(int i=0,fx=((k&1)?upt(-1):1);i<=k;i++,fx=upt(-fx))
    {
      int tp=0;
      for(int j=0;j<=i;j++)
    {
      Node tmp=x1[j]*x2[i-j];
      Node d=cal(tmp)*A[j]*B[i-j];
      tp=(tp+(ll)c[i][j]*d.x)%mod;
    }
      ans=(ans+(ll)s1[k][i]*fx%mod*tp)%mod;
    }
  ans=(ll)ans*iv%mod;
  int ml=1; for(int i=2;i<=k;i++)ml=(ll)ml*i%mod;
  ans=(ll)ans*pw(ml,mod-2)%mod;
}
int main()
{
  int op=rdn(); op=rdn();
  if(op==2)solve2(); else solve3();
  printf("%d\n",ans);
  return 0;
}

 

转载于:https://www.cnblogs.com/Narh/p/10946142.html

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

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

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

抵扣说明:

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

余额充值