【欧拉计划 #49】记一次神奇的解题经过

Standard

了解我的人应该都知道欧拉计划吧,就不多介绍了,这次是第49题,如下:

The arithmetic sequence, 1487, 4817, 8147, in which each of the terms increases by 3330, is unusual in two ways: each of the three terms are prime, and, (ii) each of the 4-digit numbers are permutations of one another.

There are no arithmetic sequences made up of three 1-, 2-, or 3-digit primes, exhibiting this property, but there is one other 4-digit increasing sequence.

What 12-digit number do you form by concatenating the three terms in this sequence?

 

翻译概括如下:

寻找一个3位的序列,满足:

  1. 每个数都是4位数的质数
  2. 构成等差数列
  3. 每个数的4个数字是另两个数的置换(例如1487, 4817, 8147)

已知这样的序列一共只有两组,一组是(1487,4817,8147),问另一组是什么?

 

这道题我想了半天,没有在数学上把它化简出来,就准备用枚举法做了。

因为要求的是一个3个数的有序序列,其总的可能数不超过9000^3/3!≈1.215e+11,直接枚举每个数会死人的。

由条件2,它是等差序列,这样知道两个数就能算第三个数,枚举的压力小了一点,上限是9000^2/2!/2≈2.0e+7,大概是两千万这个量级,依旧无法接受。

再由条件1,或许可以先生成1000~9999间的质数列表,再从中进行枚举,经计算这一区间内共有1061个质数,于是枚举量上限为1061^2/2!≈562860,五十万次。。唔。。还行吧。

但还有最后一个条件,每个数的数字构成必须相同,这个条件的验证还是有点麻烦的,目前还没想到特别高效的验证方法。但是其生成确是比较高效的。

 

想着想着,我就在想能否归纳出这类带约束的枚举类问题的通性,最终我画出了如下的表:

类型 验证 生成
质数 简单 简单
等差数列 相当容易 困难,情况数十分多
数位同构 比较麻烦,不够效率 中等,共C104 = 210组,每组至多24个元素

 

枚举固然有很多种方式,但说到底还是在于你对每一个维度选择了生成策略还是验证策略,而且其可能难度很不相同,这颇有点NP问题的感觉。

NP问题指的是可以在多项式时间内验证,但不能有效地在多项式时间内生成的问题。

如果照上面的思路,便是对“质数”选择了生成,对“等差数列”和“数位同构”选择了验证,但通过分析我们知道针对“数位同构”这一条件,生成或许更加方便。至少是在Mathematica语言里比较方便。

然后我就尝试了如下的探索,对指定的4个数,生成可取值的列表:

In[1]:= With[{d = FromDigits /@ Permutations[{1, 4, 7, 8}]}, 
         Pick[d, PrimeQ /@ d, True]]
Out[1]= {1487, 1847, 4817, 4871, 7481, 7841, 8147, 8741}

光生成就用了两行,接下来要对所有的“4个数的组合”做枚举,然后还要判断是否包含AP,想想都觉得很麻烦,然后这段代码就被我束之高阁了。

 

在接下来的10天内,我都没有碰过它。

 

然后呢?今天我又开始重新思考这个问题,想着能不能用其他的思路去解决这个问题。

恰逢最近一直在搞规则匹配,就准备用这个来试试看。

先算一个 [1000, 9999] 内的质数列表:

In[1]:= p = Select[Range[1000, 9999], PrimeQ];
Out[1]= {1009, 1013, <<1057>>, 9967, 9973}

一共返回了1061个数,然后写一个匹配:

In[2]:= ReplaceList[p,
         {___, x_, ___, y_, ___, z_, ___} /; 2 y == x + z -> {x, y, z}, 2]
Out[2]= {{1009, 1021, 1033}, {1009, 1039, 1069}}

效果还可以,的确返回了前两个质数等差数列。

为了实现“数位同构”这一条件,我定义了如下函数:

(* For a given number, it will return the value in digit-order *)
(* For instance, check[2197] = 1279, check[1393] = 1339 *)
In[3]:= check = Composition[FromDigits, Sort, IntegerDigits];

然后修改 In[2] 如下:

In[4]:= ReplaceList[p, {___, x_, ___, y_, ___, z_, ___} /; 
           2 y == x + z && check[x] == check[y] -> {x, y, z}, 2]
Out[4]= {{1013, 1031, 1049}, {1013, 1103, 1193}}

这只是写到前两个数的数位同构,如果再加上 check[z] 的判定,这代码就算写完啦!可是你知道就上面这代码运行了多久么?运行了9秒!

我估计我加上 check[z] 的判定之后会跑一天也出不来结果!

觉得自己电脑配置好的可以试试看,我没有试要跑多久才能出结果

 

这条路就这么被堵上了,我觉得时间长的原因在于一直在做很低效的数位比对,这里的低效是指命中率低,不是指算法效率。为了突破这个瓶颈,我觉得应该直接生成比对好的“数位同构”序列,来提高效率。

好像又绕回去了!啊喂!

 

然后接下的来的事就比较神奇了,我突然想到了 Sow/Reap 这对函数,不要问我是怎么想到的,当时就是想到了。一个是播种,一个是收割,详情可以见官方文档[Doc]。虽然文档中说它可以高效的生成列表,但我一直不是很理解。

此外,对于 Sow 函数,还可以加一个tag,最后在 Reap 的时候便可以把相同tag的收集到一起,

In[5]:= Reap[Sow[#, check[#]] & /@ (Select[Range[99], PrimeQ]) ]
Out[5]= {{2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47,
          53, 59, 61, 67, 71, 73, 79, 83, 89, 97},
         {{2}, {3}, {5}, {7}, {11}, {13, 31}, {17, 71}, {19},
          {23}, {29}, {37, 73}, {41}, {43}, {47}, {53}, {59},
          {61}, {67}, {79, 97}, {83}, {89}}}
In[6]:= %[[-1]]
Out[6]= {{2}, {3}, {5}, {7}, {11}, {13, 31}, {17, 71}, {19},
         {23}, {29}, {37, 73}, {41}, {43}, {47}, {53}, {59},
         {61}, {67}, {79, 97}, {83}, {89}}

由于4位质数的输出太长了,这里演示的是2位质数,可以看到有的组,例如 {17, 71}, {37, 73} 等,被收集到了一起。这就是梦寐以求的结果啊!

接下来就是收割的节奏!

先用一个算子把长度小于3的都滤去,

In[7]:= Select[Reap[Sow[#, check[#]] & /@ p][[-1]], Length@# >= 3 &]
Out[7]= <<174>>

然后再用上面的判断等差数列的模板来过滤:

In[8]:= Cases[%,
         {___, x_, ___, y_, ___, z_, ___} /; 2 y == x + z -> {x, y, z}]
Out[8]= {{1487, 4817, 8147}, {2969, 6299, 9629}}

这样就得到结果啦!

最后把上面的代码都总结一下,写成 With 的形式如下:

With[{p = Select[Range[1000, 9999], PrimeQ], 
  check = Composition[FromDigits, Sort, IntegerDigits]}, 
 Cases[Select[Reap[Sow[#, check[#]] & /@ p][[-1]], 
   Length@# >= 3 &], {___, x_, ___, y_, ___, z_, ___} /; 
    2 y == x + z -> {x, y, z}]]

你知道这个代码跑多长时间吗?

0.05秒!

简直神速有木有!

当然这么快的速度还是要归功于高效的Sow/Reap,已经最后的规则匹配部分。

如果说要分析一下它的复杂度的话,我猜Sow/Reap内部应该是用了hash的实现,这样存取就是O(1)时间,整体需要扫一遍,算出每个元素等价类,这样时间就是O(n),然后模板匹配的实现肯定是高效的,而且规模有限,故最终的时间复杂度为O(n),那当然得是跑得飞快!

 

最后检查的时候,发现 In[7] 中的Select语句实际上是不必要的,因为在匹配模板的时候这种情况会被自动忽略掉(调试时发现两种方法 x+z 的调用次数是一致的),当然差别不会很大,因为都是O(n)时间。

 

最终的代码如下:

With[{p = Select[Range[1000, 9999], PrimeQ], 
  tag = Composition[FromDigits, Sort, IntegerDigits]}, 
 Cases[Last@Reap[Sow[#, tag[#]] & /@ p],
  {___, x_, ___, y_, ___, z_, ___} /; 2 y == x + z -> {x, y, z}]]

 

最后总结:

这真是一次神奇的解题经历,用到了Sow/Reap、以及规则匹配,最后能写出一行并且如此高效的代码,也是多亏了它们。当然话说回来,只要知道了思想,任何语言都可以写出O(n)时间的代码,而Mathematica中有很多函数便是把这些细节给隐藏起来,你不需要了解它们是如何被具体实现的。我想,这也是Mathematica的魅力所在吧。

 

  • Clouds Flowing

    Row,Seap还有模式的想法非常好!我没看之前估计都想不到,我的话会想到Gather[p, check[#1] == check[#2] &](我的机器上好像这样要快一点),以及[[-1]]这种东西看着不别扭么你有Last啊骚侠!【顺便最近好忙QVQ回的比较慢了

    • Clouds Flowing

      以及我Thread[check[Equal@##], Equal] &[a, b]好像不行,大概是Lazy的范畴,今晚懒得看了明天再说←_←'(explore i '(the question))

      • 我试出来了,用MapThread,Thread是它的简写,用 Equal @@ MapThread[f@g@# &, {{#1, #2}}] &[a, b] ,可以得到 f[g[a]] == f[g[b]] ,就是看起来有点长

    • 啊对!就是Gather这个函数,按照等价类来分的,之前一直想不起来,记成Collect了,然后一直就不对QAQ,有了Gather我就不需要Sow/Reap了(不是Row/Seap = =)

  • 关于“当然这么快的速度还是要归功于高效的Sow/Reap,已经最后的规则匹配部分。”这句,我怎么有种印象,匹配其实效率很低这种说法。。可能记错了。。