jzoj 6065.【NOI2019模拟2019.3.18】One? One! fft

Description

Input

Output

Sample Input

输入1:
2 1048

输入2：
4 31415926

Sample Output

输出1：
1

输出2：
942

Data Constraint

分析：
我们设$num[i]$为$i$个$1$组成的数字，显然题目要求${\sum }_{i=2}^l\frac{n}{num[i]}$的值，其中除为整除。
对式子进行转化得到，
$$\sum _{i=2}^l\frac{n}{num[i]}=\sum _{i=2}^l\frac{9n}{10^i-1}$$

其中，我们令$N=a\ast 10^d+b$，其中$b<10^d$。
$$\frac{N}{10^d-1}=\frac{a\ast 10^d+b}{10^d-1}=\frac{a\ast (10^d-1)+a+b}{10^d-1}=a+\frac{a}{10^d-1}+\frac{b}{10^d-1}$$
其中，$\frac{a}{10^d-1}$可以递归下去。相当于每次操作使得$N$减掉最低的$d$位。
我们先不考虑求$b$的部分。对于某一个的$d$，这个数的贡献相当于去掉$d$位，去掉$2\ast d$位……后的数求和。
假设原来的数中的第$i$位，去掉若干位后是新数中的第$j$位，那么显然有$(i-j)modd=0$。贡献的次数为$i-j$的大于$1$的约数个数，设为$d[i-j]$。
设$f[i]$为原数中第$i$位的贡献和，有
$$f[i]=a[i]\ast \sum _{j=1}^{i}10^{j-1}\ast d[i-j]$$
因为$10^j$太大，考虑换一种设法，设$f[j]$为新数中第$j$位的贡献和，有
$$f[j]=10^{j-1}\ast \sum _{i=1}^{n}a[i]\ast d[i-j]$$
这是一个差卷积，把$a$翻转后FFT然后再把$f$翻转。注意到前面有$10^{j-1}$，所以求和相当于把$f$的每一位依次输出。

考虑$b$的部分，相当于每$d$位分一组的数相加。数据是随机的，因为$b<=10^d-1$，所以每一组的贡献$\in [0,1]$。而我们只需要保留整数部分，所以我们只要精确到小数点后10位，每一组误差不超过$10^{-10}$，显然不会影响到整数部分，然后直接求和即可。

因为最后要除以$10^d-1$，当$d>=10$时，我们同样可以近似处理，直接除以$10^d$即可，误差是$10^{-d}$级的。小的数可以直接暴力求。

复杂度大概是$O(10\ast l\ast logl)$。

代码：

#include <iostream>
#include <cstdio>
#include <cmath>
#include <cstring>
#include <algorithm>
#define LL long long

const int maxn=1e6+7;
const double pi=acos(-1);
const LL mod=4294967296;

using namespace std;

int n,len;
int r[maxn];
LL S,cnt;
LL a[maxn],b[maxn],c[maxn],d[maxn],f[maxn],s[maxn],bit[maxn];

struct rec{
	double x,y;
}dft[maxn],idft[maxn],w[maxn];

rec operator +(rec a,rec b)
{
	return (rec){a.x+b.x,a.y+b.y};
}

rec operator -(rec a,rec b)
{
	return (rec){a.x-b.x,a.y-b.y};
}

rec operator *(rec a,rec b)
{
	return (rec){a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x};
}

rec operator !(rec a)
{
	return (rec){a.x,-a.y};
}

void prework()
{
	for (int i=1;i<=n;i++)  //处理出读入数据，可以近似看做随机
	{
		a[i]=s[i-1]/1024%10;
		s[i]=(s[i-1]*(LL)747796405%mod+mod-1403630843)%mod;
	}
	reverse(a+1,a+n+1); //输入是从高位向低位编号，先把他反过来
	int tmp=0;
	for (int i=1;i<=n;i++) //对这个数乘9
	{
		int last=tmp;
		tmp=(a[i]*9+last)/10;
		a[i]=(a[i]*9+last)%10;
	}
	if (tmp) a[++n]=tmp;
	reverse(a+1,a+n+1); //重新反回去，做FFT
	for (int i=2;i<=n;i++) //预处理d[i]
	{
		for (int j=i;j<=n;j+=i) d[j]++; 
	}
}

void fft(rec *a,LL f)
{
    for (LL i=0;i<len;i++)
    {
        if (i<r[i]) swap(a[i],a[r[i]]);
    }
    w[0]=(rec){1,0};
    for (LL i=2;i<=len;i*=2)
    {
        rec wn=(rec){cos(2*pi/i),f*sin(2*pi/i)};
        for (LL j=i/2;j>=0;j-=2) w[j]=w[j/2];
        for (LL j=1;j<i/2;j+=2) w[j]=w[j-1]*wn;
        for (LL j=0;j<len;j+=i)
        {
            for (LL k=0;k<i/2;k++)
            {
                rec u=a[j+k],v=a[j+k+i/2]*w[k];
                a[j+k]=u+v;
                a[j+k+i/2]=u-v;
            }
        }
    }
}

void FFT(LL *a,LL *b,LL *c,int n,int m)
{
    len=1;
    int k=0;
    while (len<(n+m)) len*=2,k++;
    for (int i=0;i<len;i++)
    {
        r[i]=(r[i>>1]>>1)|((i&1)<<(k-1));
    }
    for (int i=0;i<len;i++)
    {
        if (i<n) dft[i].x=a[i]; else dft[i].x=0;
        if (i<m) dft[i].y=b[i]; else dft[i].y=0;
    }
    fft(dft,1); 
    for (int i=0;i<len;i++)
    {
        int j=(len-1)&(len-i);
        rec da,db;
        da=(dft[i]+(!dft[j]))*(rec){0.5,0};
        db=(dft[i]-(!dft[j]))*(rec){0,-0.5};
        idft[i]=da*db;
    }
    fft(idft,-1);
    for (int i=0;i<m;i++) c[i]=trunc(idft[i].x/len+0.5);
}

int main()
{
	freopen("one.in","r",stdin);
	freopen("one.out","w",stdout);
	scanf("%d%lld",&n,&s[0]);			
	prework();	
	FFT(a,d,f,n+1,n+1); //卷积得出f，然后取反
	reverse(f+1,f+n+1);
	reverse(a+1,a+n+1);	//再次把a反过来
	bit[0]=1;
	for (int i=1;i<=10;i++) bit[i]=bit[i-1]*10;		
	for (int i=2;i<=10;i++) //对于小的数，分母近似处理可能会出错，暴力处理
	{
		S=0;
		for (int j=1;j<=n;j+=i)
		{
			for (int k=1;k<=i;k++) S+=a[j+k-1]*bit[k-1];
	    }
	    f[1]+=S/(bit[i]-1);
	}
	for (int i=11;i<=n;i++)
	{
		for (int j=1;j<=12;j++) b[j]=0; //记录前面的一些位的和
		for (int j=1;j<=n;j+=i)
		{
			for (int k=1;k<=11;k++) b[k]+=a[j+i-11+k-1]; //求和记录再b中，注意高位在区间的右边
		}
		for (int j=1;j<=11;j++) //移位
		{
		    b[j+1]+=b[j]/10;
		    b[j]=b[j]%10;
		}
		f[1]+=b[12]; //整数部加入到答案中，不同的i之间是相互独立的
	}
	for (int i=1;i<=n;i++)
	{
		f[i+1]+=f[i]/10;
		f[i]=f[i]%10;
	}	
	int num=n;
	while ((!f[num]) && (num>1)) num--;
	for (int i=num;i>0;i--) printf("%d",f[i]); //移位后输出f
}