首页 > 算法优化 > Euler Project之Retraction【算法优化向】

Euler Project之Retraction【算法优化向】

2013年12月8日 发表评论 阅读评论

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

题目来源还是是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

  • Nima

    时间复杂度还是太大,lz最后算出吗?

    • 写这篇东西就是为了敲公式,也没怎么细想,最后也忘了这道题了。。你有啥好想法么?