好吧，其实我是在新博客这边练习敲公式来着。。。

题目来源还是是Euler Project,之前官方连续公布了三道类似的题目，都是基于一种叫做Retraction数的东西，它的定义是：

对于所有整数\(n>1\)，定义一族函数\(f_{n,a,b}(x)=mod(ax+b,n)\)，\(x,a,b\)均为整数，且满足0<\(n\)<\(a\),\(0\leq b < n\),\(0\leq x < n\)
如果对所有\(x\)均有：

\(f_{n,a,b}(f_{n,a,b}(x))=mod(f_{n,a,b}(x),n)\)

那么，我们管\(f_{n,a,b}(x)=mod(ax+b,n)\)叫做一个Retraction；
我们需要计算的是\(R(n)\)，定义为包含\(n\)的Retraction的个数。

虽然这么说，但是计算\(R(n)\)并不是题目的最终结果，只是之前公布的三道题目全都是跟这个\(R(n)\)有关的，所以我这里就来算一下怎么算\(R(n)\)最快，好吧，不一定是最快的，但是我觉得还是挺快的。。

如果真正运用到题目的话，倒是不会用到我这个计算\(R(n)\)的东西，因为配合题目的别的条件，可以更快的算出题目真正要求的东西。不过我就是要算\(R(n)\)，谁也别拦我在这里写公式！！以前在blog.com被不能写数学公式的烦躁感压抑太久了！！【万一哪一天Latex插件失败了，我这个页面将会显示一大堆乱码。。。想想就好可怕。。】

首先，我们先用Python写一个最简单(思维量上)的版本，就是把题目翻译一下的那个版本:

def f(n,a,b,x):
    return mod(a*x+b,n)

def IsRetraction(a,b,n):
    for x in range(n):
        if f(n,a,b,f(n,a,b,x)) != mod(f(n,a,b,x),n):
            return False
    return True

def R(n):
    counter = 0;
    for a in range(1,n):
        for b in range(0,n):
            if IsRetraction(a,b,n):
                counter += 1;
    return counter;

有兴趣的同学可以自己去跑一下这段代码，慢的不可理喻！！因为随着n的增大，计算复杂度是\(o(n^3)\)

所以我们可以对题目的条件进行一下推导：

\(\because\)\(f_{n,a,b}(f_{n,a,b}(x))=mod(f_{n,a,b}(x),n)\)

\(\therefore\)\(mod(af_{n,a,b}(x)+b,n)=mod(f_{n,a,b}(x),n)\)

\(\therefore\)\(mod((a-1)f_{n,a,b}(x)+b,n)=0\)—————————–※

\(\therefore\)\(mod((a-1)mod(ax+b,n)+b,n)=0\)

代入\(x=0\)，其中引文\(b\)<\(n\)，所以有\(mod(b,n)=b\),所以：
\(mod((a-1)b+b,n)=mod(ab,n)=0\)

所以我们得到第一个结论：\(ab=kn\)

上面的结论成立必然会有在\(x=0\)的时候成立（因为我们就是把\(x=0\)代进去推出来的嘛~），所以我们可以以此为“导火索”，采用数学归纳法！！

也就是说，我们只要可以从\(x\)成立推导出\(x+1\)成立的条件，因为\(x=0\)已经成立了，所以必然有对于一切\(x\)成立！！

上面的※式中，我们有\((a-1)mod(ax+b,n)+b=kn\)，【\(k\)是整数,下面不再累述】然后我们把\(x+1\)代进去:

\((a-1)mod(ax+b+a,n)+b=k’n\)—————————–※2

同余计算里面有以下结论：
\(\because\)\(mod(a+b,n)=mod(mod(a,n)+mod(b,n),n)\)

\(\therefore\)\(mod(ax+b+a,n)=mod(mod(ax+b,n)+mod(a,n),n)=mod(mod(ax+b,n)+a,n)\)

所以之前的※2式子就变成：\((a-1)mod(mod(ax+b,n)+a,n)+b=k’n\)

因为上面红色式子对\(x\)成立的时候必然有：\(mod(ax+b,n)=\dfrac{kn-b}{a-1}\)

所以代入之后又：\((a-1)mod(\dfrac{kn-b}{a-1}+a,n)+b=k’n\)

\(\therefore\)\((a-1)mod(\dfrac{kn-b+a^2-a}{a-1},n)+b=k’n\)

这个时候\(a-1\)不能带进去，所以两边先对\(n\)取模：

\(mod((a-1)mod(\dfrac{kn-b+a^2-a}{a-1},n)+b,n)=0\)

这个时候\(a-1\)是可以带进里面去的，所以就变成:

\(mod(mod((kn-b+a^2-a),n)+b,n)=0\)

\(mod(mod((a^2-a-b),n)+b,n)=0\)

\(mod((a^2-a-b+b),n)=mod(a^2-a)=0\)

之前我推到这里简直震惊了！！因为我们要判断\(f_{a,b,n}(x)\)是不是Retraction，只要判断下面这两个式子是否成立就可以了！！

\(\left\{\begin{array}{ll}mod(ab,n)=0\\mod(a^2-a,n)=0 \end{array}\right.\)

于是乎，代码就变成这个样子了：

def R_version2(n):
    counter = 0;
    for a in range(1,n):
        if mod(a*a-a,n) != 0:
            continue
        for b in range(0,n):
            if mod(a*b,n) == 0:
                counter += 1;
    return counter;

代码瞬间就变成只有一个函数了~
我在我的破机子上跑，计算R(1000)，最开始那个代码要20.1578052518秒，现在只要0.0142358875219，这个已经算时间计算误差了，基本就是秒出了。
但是，其实还可以优化的！！

之前的思路是对\(a\)从1扫到n-1，然后判断是否满足\(mod(a^2-a,n)=0\)，满足，再把\(b\)从0扫到\(n-1\)，判断是否满足\(mod(ab,n)=0\)，如果满足，那么这一组\(a,b,n\)就是一个Retraction。

但是我们再看看第二个条件：假设我们已知了\(a\)满足\(mod(a^2-a,n)=0\)，应该可以直接算出有多少个\(b\)满足\(mod(ab,n)=0\)的对吧！！没错，就是a和n的最大公约数！！！

可能有人想不通，下面我详细解释一下：

原因是这样子的，我们要计算满足\(mod(ab,n)=0\)的\(b\)，也就是要寻找满足\(ab=kn\)的\(b\)，也就是说，\(b\)最小为0，再其次的话，第二小的\(b\)里面一定要包含\(n\)的质约数与\(a\)的质约数的差集，比如\(n=40\)，\(a=16\)，那么最小的\(b\)为0，\(n=2^3\times 5\),\(a=2^4\)，所以第二小的\(b\)的约数是集合{2,2,2,5}减去集合{2,2,2,2}【这里说集合好像不是很对，因为集合不能有重复的元素的~嘛~别在意啦。。】，也就是5啦！而且，\(b\)取5，\(k\)就是去2，再然后每个\(b\)就不断地增加就是了，所以得到的\(b\)就是0,5,10,15,20,25,30,35，八个。

那么最小的非零\(b\)怎么算，就是用\(n,a\)的最小公约数除以\(a\)，也就是\(\dfrac{LCM(a,n)}{a}\)，最大的\(b\)就是\(n-\dfrac{LCM(a,n)}{a}\)，说白了就是\(b\)的取值是:

\(0,\dfrac{LCM(a,n)}{a},\dfrac{LCM(a,n)}{a}\times 2…n-\dfrac{LCM(a,n)}{a}\)

显然，因为公差是\dfrac{LCM(a,n)}{a}，所以最后得到的\(b\)的个数就是\(\dfrac{n}{\dfrac{LCM(a,n)}{a}}=\dfrac{an}{LCM(a,n)}\)，小学奥数告诉我们，两个数最小公倍数乘最大公约数就是两个数的乘积，所以\(\dfrac{an}{LCM(a,n)}=GCD(a,n)\)!!

证明完毕！！Q.E.D!!

最后代码就会变成：

def R_version2(n):
    counter = 0;
    for a in range(1,n):
        if mod(a*a-a,n) != 0:
            continue
        counter += gcd(a,n);
    return counter;

简洁吧，这都是数学的功劳！！现在的时间复杂度就是\(o(n)\)了，如果想要再加速，那么就要想怎么快速的计算出满足\(mod(a^2-a,n)=0\)的a来了。。。有人有好的思路么？
嗯，好吧，先写到这里吧！敲了半天公式，一本满足~

【完】

本文内容遵从CC版权协议，转载请注明出自http://www.kylen314.com