Project Euler 326

题目

Modulo Summations

Let \(a_n\) 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);
}