Project Euler 515

题目

Dissonant Numbers

Let \(d(p,n,0)\) be the multiplicative inverse of \(n\) modulo prime \(p\), defined as \(n\times d(p,n,0)=1\bmod p\).

Let \(d(p,n,k) = \sum_{i=1}^n d(p,i,k−1)\) for \(k\ge1\).

Let \(D(a,b,k) = \sum(d(p,p-1,k) \bmod p)\) for all primes \(a\le p<a+b\).

You are given:

• \(D(101,1,10) = 45\)

• \(D(10^3,10^2,10^2) = 8334\)

• \(D(10^6,10^3,10^3) = 38162302\)

Find \(D(10^9,10^5,10^5)\).

解决方案

根据上面的定义，可以写出如下式子：

\(\begin{aligned} &d(p,n,0)=n^{-1}\\ &d(p,n,1)=\sum_{i_1=1}^n i_1^{-1}\\ &d(p,n,2)=\sum_{i_2=1}^n\sum_{i_1=1}^{i_2} i_1^{-1}\\ &d(p,n,3)=\sum_{i_3=1}^n\sum_{i_2=1}^{i_3}\sum_{i_1=1}^{i_2} i_1^{-1}\\ &\dots \end{aligned}\)

最终有

\(\begin{aligned} d(p,n,k)&=\sum_{i_k=1}^n\sum_{i_{k-1}=1}^{i_k}\dots\sum_{i_2=1}^{i_3}\sum_{i_1=1}^{i_2} i_1^{-1}\\ &=\sum_{i_1=1}^n i_1^{-1}\cdot\left(\sum_{i_2=i_1}^n\sum_{i_3=i_2}^n\dots\sum_{i_{k-1}=i_k}^n\sum_{i_k=i_{k-1}}^n1\right) \end{aligned}\)

我们分析求和项\(\displaystyle{\sum_{i_2=i_1}^n\sum_{i_3=i_2}^n\dots\sum_{i_{k-1}=i_k}^n\sum_{i_k=i_{k-1}}^n1}\)的值。可以发现，这个求和式的结果等价于如下组合问题的结果：

求有多少个长度为\(k-1\)的序列\((i_2,i_3,\dots,i_{k-1},i_k)\)，满足如下条件：

• \(i_1\le i_2,i_k\le n\)
• \(\forall j\in[3,k],i_{j-1}\le i_j\)均成立。

将它们的条件转换成如下条件，得到的结果是等价的（想象一下，令\(i'_j=i_j+j-1\)）。

• \(i_1\le i_2,i_k\le n+k-2\)
• \(\forall j\in[3,k],i_{j-1}< i_j\)均成立。

这个问题是一个普通的组合问题，它的答案为\(\dbinom{n+k-i_1-1}{k-1}\).

因此可以得到：

\[d(p,n,k)=\sum_{i=1}^n i^{-1}\cdot\dbinom{n+k-i-1}{k-1}\]

代入\(n=p-1\)，得到

\(\begin{aligned} d(p,p-1,k)&\equiv\sum_{i=1}^{p-1} i^{-1}\cdot\dbinom{p+k-i-2}{k-1}\\ &\equiv\sum_{i=1}^{p-1} i^{-1}\cdot\dfrac{(p-i)(p-i+1)\dots(p+k-i-3)(p+k-i-2)}{(k-1)!}\\ &\equiv-\dfrac{1}{k-1}\sum_{i=1}^{p-1}\dfrac{(p-i+1)\dots(p+k-i-3)(p+k-i-2)}{(k-2)!}\\ &\equiv-\dfrac{1}{k-1}\sum_{i=1}^{p-1}\dbinom{p+k-i-2}{k-2}\\ &\equiv-\dfrac{1}{k-1}\sum_{i=1}^{p-1}\dbinom{i+k-2}{k-2}\\ &\equiv-\dfrac{1}{k-1}\left(p\cdot \dfrac{\binom{p+k-2}{k-2}}{k-1}-1\right)\\ &\equiv\dfrac{1}{k-1} &\pmod p \end{aligned}\)

最终，直接暴力计算\(D(a,b,p)\)的值即可。

代码

#include <bits/stdc++.h>
# include "tools.h"
using namespace std;
typedef long long ll;
const int A=1000000000,B=100000,K=100000;
int main(){
    ll ans=0;
    for(int p=next_prime(A-1);p<A+B;p=next_prime(p)){
        ans+=mod_inverse(K-1,p);
    }
    printf("%lld\n",ans);
}