Project Euler 402

题目

Integer-valued polynomials

It can be shown that the polynomial \(n^4 + 4n^3 + 2n^2 + 5n\) is a multiple of \(6\) for every integer \(n\). It can also be shown that \(6\) is the largest integer satisfying this property.

Define \(M(a, b, c)\) as the maximum \(m\) such that \(n^4 + an^3 + bn^2 + cn\) is a multiple of \(m\) for all integers \(n\). For example, \(M(4, 2, 5) = 6\).

Also, define \(S(N)\) as the sum of \(M(a, b, c)\) for all \(0 < a, b, c \le N\).

We can verify that \(S(10) = 1972\) and \(S(10000) = 2024258331114\).

Let \(F_k\) be the Fibonacci sequence:

\(F_0 = 0, F_1 = 1\) and

\(F_k = F_{k-1} + F_{k-2}\) for \(k \ge 2\).

Find the last \(9\) digits of \(\sum S(F_k)\) for \(2 \le k \le 1234567890123\).

解决方案

令\(p(n)=n^4 + an^3 + bn^2 + cn,m=10^9\)。那么可以得到\(\gcd(p(-1)+p(1),p(-2)+p(2))=\gcd(2b+2,8b+32)=2\gcd(b+1,12)\)。

也就是说，\(M\in\{1,2,3,4,6,8,12,24\}\)。

由于\(24=3\cdot 2^3\)，经过多次打表实验，可以得到如下结论：

\(\begin{aligned} &3\mid M(a,b,c)\Longleftrightarrow b\equiv 2\pmod 3 \land(a+c)\equiv 0\pmod 3;\\ &2\mid M(a,b,c)\Longleftrightarrow (a+b+c)\equiv 1\pmod 2;\\ &4\mid M(a,b,c)\Longleftrightarrow a\equiv c\equiv0\pmod 2\land (a+b+c)\equiv 3\pmod 4\\ &8\mid M(a,b,c)\Longleftrightarrow a\equiv c\equiv2\pmod 4\land (a+b+c)\equiv 7\pmod 8\\ \end{aligned}\)

那么，接下来的任务是以\(O(1)\)的时间复杂度计算\(S(N)\)。由于\(a,b,c\)的值都是模\(24\)相等的，因此我们可以将整个\(N\times N\times N\)的输入以\(24\times24\times24\)的大小进行划分。令\(q=\left\lfloor\dfrac{N}{24}\right\rfloor,r=N\%24\)，那么每一个\(24\times24\times24\)方块内部答案的总和都是相等的。

因此，对于模数\(m\)，序列\(S(n)\% m\)的周期为\(24m\)。可知序列\(F_n\%m\)也是周期性的，并且\(24m\)一定是其中一个周期，因此我们考虑枚举\(F_n\%24m'\)的循环节。

此外，由于\(m=2^9\cdot 5^9\)，因此我们可以分成两部分去求解这个问题：分别求解\(S(n)\% 2^9\)和\(S(n)\% 5^9\)的值后，再用中国剩余定理进行合并，最终得到答案。

代码

#include <bits/stdc++.h>
# include "tools.h"
using namespace std;
typedef long long ll;

const ll N=1234567890123;
const ll mod = 1000000000;

const int B=24;
int M(int a,int b,int c){
    int w=1;
    if(b%3==2&&(a+c)%3==0) w*=3;
    if(a%4==2&&c%4==2&&(a+b+c)%8==7) w*=8;
    else if(a%2==0&&c%2==0&&(a+b+c)%4==3) w*=4;
    else if((a+b+c)%2==1) w*=2;
    return w;
}
int s[B][B][B];
ll z[4][B]={0};
void init(){
    for(int a=1;a<=B;a++)
        for(int b=1;b<=B;b++)
            for(int c=1;c<=B;c++){
                int m=M(a,b,c);
                for(int j=0;j<B;j++){
                    z[0][j]+=m;
                    z[1][j]+=m*((j>=a)+(j>=b)+(j>=c));
                    z[2][j]+=m*((j>=max(a,b))+(j>=max(a,c))+(j>=max(b,c)));
                    if(j>=a&&j>=b&&j>=c) z[3][j]+=m;
                }
            }
}
ll S(ll n,ll mod){
    ll q=n/B,r=n%B;
    q%=mod;
    return (z[0][r]*qpow(q,3,mod)%mod+z[1][r]*qpow(q,2,mod)%mod+z[2][r]*q%mod+z[3][r])%mod;
}
ll cal(ll n,ll mod){
    ll q=n/(mod*B),r=n%(mod*B);
    ll T=min(mod*B,n);
    ll a=0,b=1;
    ll all=0,res=0;
    for(int i=0;i<T;i++){
        ll k=S(b,mod);
        all=(all+k)%mod;
        if(i<r) res=(res+k)%mod;
        ll t=(a+b)%(mod*B);
        a=b;b=t;
    }
    ll ans=(q%mod*all%mod+res-S(1,mod)+mod)%mod;
    return ans;

}
int main(){
    init();
    vector<ll>mod_list,res_list;
    for(auto &[p,e]:fact(ll_2_mpz(mod))){
        mod_list.push_back(qpow_mpz(mpz_2_ll(p),e));
    }
    for(ll p:mod_list) res_list.push_back(cal(N,p));
    ll ans=CRT(mod_list,res_list);
    printf("%lld\n",ans);
}