Project Euler 326

题目

Modulo Summations

Let $an$ be a sequence recursively defined by:$\quad a_1=1,\quad\displaystyle a_n=\biggl(\sum{k=1}^{n-1}k\cdot a_k\biggr)\bmod n$.

So the first $10$ elements of $a_n$ are: $1,1,0,3,0,3,5,4,1,9$.

Let $f(N,M)$ represent the number of pairs $(p,q)$ such that:

$$1 \le p \le q \le N \text{ and } \left(\sum_{i=p}^{q}a_i\right) \text{ mod } M=0$$

It can be seen that $f(10,10)=4$ with the pairs $(3,3), (5,5), (7,9)$ and $(9,10)$.

You are also given that $f(10^4,10^3)=97158$.

Find $f(10^{12},10^6)$.

解决方案

通过以下程序进行打表，发现序列$a_n$满足如下规律：

# include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N=200;
int s[N+4],a[N+4];
int main(){
    s[1]=a[1]=1;
    for(int n=2;n<=N;n++){
        a[n]=s[n-1]%n;
        s[n]=s[n-1]+a[n]*n;
    }
    for(int i=0;i<=N;i++)
        printf("%d%c",a[i]," \n"[i%6==5]);
}

$\begin{aligned}
&a{6i}=3i\
&a{6i+1}=4i+1\
&a{6i+2}=3i+1\
&a{6i+3}=i\
&a{6i+4}=6i+3\
&a{6i+5}=i
\end{aligned}$

那么令$\displaystyle{S(n)=\sum_{i=1}^na_i}$，那么根据这个周期性质，通过待定系数法来进行计算，不难写出：

$\begin{aligned}
&S(6i)=9i^2-i\
&S(6i+1)=9i^2+3i+1\
&S(6i+2)=9i^2+6i+2\
&S(6i+3)=9i^2+7i+2\
&S(6i+4)=9i^2+13i+5\
&S(6i+5)=9i^2+14i+5
\end{aligned}$

由于$\displaystyle{S(q)-S(p)=\sum_{i=p}^qa_i}$，因此问题可以转化成如下形式：

有多少对$(p,q)$满足如下条件：

• $0\le p< q< N$
• $S(p)-S(q)\equiv 0\pmod M$，也就是$S(p)\equiv S(q)\pmod M$

那么不难发现，我们可以统计$0\sim N$中有多少个$i$满足$S(i)\%M=j$，用数组$c$来存储个数，并统计在$c[j]$中。那么最终答案为

$$f(N,M)=\sum_{i=0}^{M-1} \dfrac{c[i]\cdot (c[i]-1)}{2}$$

对于$j\in[0,5]$，$S(6i+j)$都是二次多项式，因此不难发现$S(i+M)\equiv S(i)\pmod p$。

因此，对于$i\in[0,M),j\in[0,5]$，我们统计有多少个非负整数$k$满足$6\cdot(i+Mk)+j\le N$。并将统计的值统计到$c[S(6i+j)\%M]$即可。不难计算得到满足条件的$k$的个数为$\left\lfloor\dfrac{N-6i-j}{6M}\right\rfloor+1$。

代码

# include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const ll N=1000000000000;
const int M=1000000;
const int T=6;
ll c[M];
ll v[T];
int main(){
    for(int i=0;i<M;i++){
        v[0]=1ll*9*i*i-i;
        v[1]=v[0]+4*i+1;
        v[2]=v[1]+3*i+1;
        v[3]=v[2]+i;
        v[4]=v[3]+6*i+3;
        v[5]=v[4]+i;
        for(int j=0;j<T;j++) {
            v[j]%=M;
            c[v[j]]+=(N-j-T*i+T*M)/(T*M);
        }
    }
    ll ans=0;
    for(int i=0;i<M;i++)
        ans+=c[i]*(c[i]-1)/2;
    printf("%lld\n",ans);
}