正题

题目链接:https://gmoj.net/senior/#main/show/6065

---

题目大意

$oneness(x)$表示$x$的约数中全是$1$的数的个数，给出一个长度为$l$的随机生成的数$n$，求$$\sum _{i=1}^{n}oneness(i)$$

---

解题思路

转换一下题意就是求$$\sum _{i=2}^l\lfloor \frac{n}{\frac{10^i-1}{9}}\rfloor =\sum _{i=2}^l\lfloor \frac{9n}{10^i-1}\rfloor$$
先让$n$乘上$9$，然后把每一位的贡献拆成小数和整数部分，我们先来看整数部分

首先对于一个$\frac{k\ast 10^i}{10^j-1}$那么它的值应该是在第$i-xj+1(x\in N)$位会有一个$k$。（如$\frac{10^9}{999}=1001001.001001001001001001001001$）

然后我们看这题，$j$是$1\sim n$的每一个值，那么对于每第$i$位我们拆成$k\ast 10^i$会对第$i-x$位产生${\sigma }_0(x)$次贡献（${\sigma }_0(x)$表示$x$的约数个数）。设$a_i$表示第$i$位的数，$f_i$表示第$i$位的答案，那么有式子$$f_i=\sum _{j=0}^{\infty }{a}_{i+j}{\sigma }_0(j)$$
这个式子可以先预处理出${\sigma }_0(j)$然后把$a$反过来就是一个卷积的式子，$FFT$优化即可。

对于小数部分，因为在比较后的位数的产生贡献概率极小，又因为数字随机生成，所以我们直接将这些部分省略。我们只处理到小数的$12$位，和上面同理算出每一个小数位的贡献即可。

时间复杂度$O(n\log n)$

---

$code$

#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>
#include<cstdlib>
#define ll long long
using namespace std;
const ll N=1e6+10;
const double Pi=acos(-1);
struct complex{double x,y;complex(double xx=0,double yy=0){x=xx;y=yy;return;}
}x[N],y[N];
ll n,m,r[N],a[N],d[N],b[N],f[N];
unsigned int sed;
complex operator+(complex x,complex y)
{return complex(x.x+y.x,x.y+y.y);}
complex operator-(complex x,complex y)
{return complex(x.x-y.x,x.y-y.y);}
complex operator*(complex x,complex y)
{return complex(x.x*y.x-x.y*y.y,x.x*y.y+x.y*y.x);}
void fft(complex *f,ll op){for(ll i=0;i<m;i++)if(i<r[i])swap(f[i],f[r[i]]);for(ll p=2;p<=m;p<<=1){ll len=p>>1;complex tmp(cos(Pi/len),op*sin(Pi/len));for(ll k=0;k<m;k+=p){complex buf(1,0);for(ll i=k;i<k+len;i++){complex tt=f[i+len]*buf;f[i+len]=f[i]-tt;f[i]=f[i]+tt;buf=buf*tmp;}}}return;
}
void solve(){m=1;ll lg=0;for(;m<=2*n;m<<=1);for(ll i=0;i<=m;i++)r[i]=(r[i>>1]>>1)|((i&1)?(m>>1):0);for(ll i=1;i<=n;i++)x[i]=complex(a[n-i+1],0),y[i]=complex(d[i],0);fft(x,1);fft(y,1);for(ll i=0;i<m;i++)x[i]=x[i]*y[i];fft(x,-1);for(ll i=0;i<m;i++)f[i]=(int)(x[i].x/(double)m+0.5);return;
}
int main()
{scanf("%lld%u",&n,&sed);for(ll i=1;i<=n;i++){a[i]=sed/1024%10;sed=sed*747796405u-1403630843u;}reverse(a+1,a+n+1);ll g=0;for(ll i=1;i<=n;i++){a[i]=a[i]*9+g;g=a[i]/10;a[i]%=10;}if(g)a[++n]=g;for(ll i=2;i<=n;i++)for(ll j=i;j<=n;j+=i)d[j]++;solve();reverse(f+1,f+n+1);ll pw[11];pw[0]=1;for(ll i=1;i<=10;i++)pw[i]=pw[i-1]*10;for(ll i=2;i<=10;i++){ll S=0;for(ll j=1;j<=n;j+=i)for(ll k=1;k<=i;k++)S+=a[j+k-1]*pw[k-1];f[1]+=S/(pw[i]-1); }for(ll i=11;i<=n;i++){for(ll j=1;j<=12;j++)b[j]=0;for(ll j=1;j<=n;j+=i)for(ll k=1;k<=11;k++)b[k]+=a[j+i+k-12];for(ll j=1;j<=11;j++){b[j+1]+=b[j]/10;b[j]%=10;}f[1]+=b[12];}for(ll i=1;i<=n;i++){f[i+1]+=f[i]/10;f[i]%=10;}while(!f[n])n--;while(n)printf("%lld",f[n--]);return 0;
}