Project Euler 451

题目

Modular inverses

Consider the number $15$.

There are eight positive numbers less than $15$ which are coprime to $15: 1,2, 4, 7, 8, 11, 13, 14$.

The modular inverses of these numbers modulo 15 are: $1, 8, 4, 13, 2, 11, 7, 14$ because

$\begin{aligned}
1 \cdot 1 &\bmod 15=1\
2 \cdot 8=16 &\bmod 15=1\
4 \cdot 4=16 &\bmod 15=1\
7 \cdot 13=91 &\bmod 15=1\
11 \cdot 11=121 &\bmod 15=1\
14 \cdot 14=196 &\bmod 15=1
\end{aligned}$

Let $I(n)$ be the largest positive number $m$ smaller than $n-1$ such that the modular inverse of $m$ modulo $n$ equals $m$ itself.

So $I(15)=11$.

Also $I(100)=51$ and $I(7)=1$.

Find $\sum I(n)$ for $3\le n\le2\cdot10^7$

解决方案

将问题进行正式描述：求$I(n)$的值，其中$I(n)$是方程$x^2\equiv 1\pmod n$中，小于$n-1$的最大正整数解。方程也可以化为：

$$(x+1)(x-1)\equiv 0\pmod n\qquad(1)$$

假设$n$的分解质因数为$n=\prod_{i=1}^{k} p_i^{e_i}$，那么不难发现：

1. 对于任意一个$p_i$为奇数的$p_i^{e_i}$，要么$p_i^{e_i}\mid x-1$，要么$p_i^{e_i}\mid (x+1)$，因此方程$(x+1)(x-1)\equiv 0 \pmod {p_i^{e_i}}$只有两个解：$1,-1$。
2. 对于$2^{e}$，在保持上面的性质的基础上，可以发现$4$并不能整除$(x+1)$和$(x-1)$，但是$2$可以，因此对于方程$(x+1)(x-1)\equiv 0\pmod {2^e}$一共有$4$个解：$1,-1,2^{e-1}-1,2^{e-1}+1$.当$e\le 2$时，需要将这些值相同的解合并。

最终，对于每个$p_i^{e_i}$分量，方程$(x+1)(x-1)\equiv 0 \pmod {p_i^{e_i}}$都有多个解。通过中国剩余定理将它们进行合并成原方程$(1)$的解再取最大值即可。

不过由于数据范围$N=2\times 10^7$，比较大，因此实际上实现的时候，使用中国剩余定理先将多个形如$x^2\equiv 1\pmod {p_i^{e_i}}$方程的解两个两个地进行组合。具体方式如下：

假设$n’=\prod{i=1}^{k’}p_i^{e_i},k’< k,m=p{k’+1}^{e_{k’+1}}$，如果我们已经解出了$x^2\equiv 1 \pmod {n’}$的所有解，那么再方程的所有解$x_b$和方程$x^2\equiv 1 \pmod {m}$的所有解$x_a$（注意$x_a$要么为$1$，要么为$-1$）一一组合，得到如下方程组：

$$\left \{\begin{aligned} & x\equiv x_a\pmod {m}\\ & x\equiv x_b\pmod {n'}\\ \end{aligned}\right.$$

令$x=t_am+x_a=t_bn’+x_b$，那么得到一个以$t_a$和$t_b$为未知数的方程$t_am-t_bn’=x_b-x_a$.

通过扩展欧几里得算法，我们可以解得方程$t_a’m-t_b’n’=1$的解$(t_a’,t_b’)$。那么得到$t_a=(x_b-x_a)t_a’$

那么得到$x=t_am+x_a=(x_b-x_a)t_a’m+x_a=x_bt_a’m+x_a(1-t’_am)$。由于$x_a$能够分别取$1$和$-1$，因此每合并一个方程，都能够产生两个解：$x_bt_a’m\pm(1-t’_am)$

代码

#include <bits/stdc++.h>
# define lb(x) ((x)&(-x))
using namespace std;
typedef long long ll;
const int N=2e7;
int qpow(ll n,ll m,ll mod){
    ll a=1;
    for(;m;m>>=1){
        if(m&1) a=a*n%mod;
        n=n*n%mod;
    }
    return a;
}
int mxp[N+4];
int a[1004],c=0;
int r[14],m=0;
int main(){
    for(int p=3;p<=N;p++)
        if(!mxp[p])
            for(int i=p;i<=N;i+=p) mxp[i]=p;
    ll ans=0;
    for(int i=3;i<=N;i++){
        m=c=0;
        int pw2=lb(i);
        if((pw2&1)==0){
            a[c++]=1;
            if((pw2&2)==0){
                a[c++]=pw2-1;
                if((pw2&4)==0){
                    a[c++]=(pw2>>1)+1;
                    a[c++]=(pw2>>1)-1;
                }
            }
        }
        else a[c++]=0;
        int mul=pw2;
        for(int v=i/mul;v!=1;){
            int p=mxp[v],pw=1;
            for(;v%p==0;v/=p,pw*=p);
            int nmul=mul*pw;
            ll x=qpow(mul,pw/p*(p-1)-1,pw)*mul%nmul;
            ll y=1ll-x;
            for(int pc=c,j=0;j<pc;j++){
                a[c++]=(a[j]*y+x)%nmul;
                a[j]=(a[j]*y-x)%nmul;
            }
            mul=nmul;
        }
        int mx=0;
        for(int j=0;j<c;j++){
            if(a[j]<0) a[j]+=i;
            if(a[j]!=i-1) mx=max(mx,a[j]);
        }
        ans+=mx;
    }
    printf("%lld\n",ans);
}