欧拉计划（Project Euler）解法【1~25】

我做这个练习的主要目的是为了熟悉 Mathematica 的编程和函数，而不是练习算法，因此会充分利用 Mathematica 的内置函数和特性，尽量少造轮子．

P1: Multiples of 3 and 5

判断一个数是否是 3 的倍数或者是 5 的倍数．

Select[Range[999], Mod[#, 3] == 0 || Mod[#, 5] == 0 &] // Total

后来发现，Mathematica 有内置的 Divisible 函数，比使用 Mod 函数要快 20%~30%．

Select[Range[999], Divisible[#, 3] || Divisible[#, 5] &] // Total//Timing

P2: Even Fibonacci numbers

四百万以下的偶 Fibonacci 数相加．

第一种方法是选一个适当的上界：

Select[Fibonacci /@ Range[3, 50, 3], # <= 4000000 &] // Total

第二种是常规的累加：

Module[{n = 3, res = 0, t = Fibonacci[3]},
  While[t <= 4*^6, res += t; n += 3; t = Fibonacci[n]];
  res]

第三种是先找到满足条件的最大的 Fibonacci 数的指针，再相加：

Sum[Fibonacci[i], {i, 3, NestWhile[(# + 1) &, 1, Fibonacci[#] <= 4*^6 &], 3}]

时间上，第一种 < 第二种 < 第三种．

但是，如果第一种方法的上界取 100，则变为 第二种 < 第一种 < 第三种．

在网上见到的一种做法，时间介于第二种和第三种之间：

Last[{1, 2, 3, 0} //.
  {{i_, j_, k_, n_} /; j<4000000 -> {j+k, 2k+j, 3k+2j, n+j}}]

这种方法利用了 //.（重复替换）来进行迭代，/;（条件）表示终止条件．

受此启发，我重写了一个更快的算法．

实际上，我们可以稍微对 Fibonacci 的递推公式做一下推导：

$$\begin{aligned} F_n &= F_{n-1} + F_{n-2} \\ &= 2 F_{n-2} + F_{n-3} \\ &= 3 F_{n-3} + 2 F_{n-4} \\ &= 3 F_{n-3} + F_{n-4} + F_{n-5} + F_{n-6} \\ &= 4 F_{n-3} + F_{n-6} \end{aligned}$$

于是上面的做法可以改写成：

Last[{0, 2, 0} //.
  {{i_, j_, n_} /; j<4000000 -> {j, 4j+i, n+j}}]

这个算法的速度比第二种稍快一些．

类似的想法，也可以用 NestWhile 函数实现：

NestWhile[{#[[1]] + 3, Fibonacci[#[[1]]], #[[3]] + #[[2]]} &,
  {0, 0, 0}, #[[2]] < 4000000 &] // Last

但是这种方法比上面一种方面慢了一倍，甚至比第三种还要慢，感觉是被 Part 函数拖慢了速度．

P3: Largest prime factor

找最大的素因子．

First /@ FactorInteger[600851475143] // Max

P4: Largest palindrome product

两个三位数相乘得到的最大回文数．

Table[i*j, {i, 100, 999}, {j, 100, i}] // Flatten // ReverseSort //
  SelectFirst[PalindromeQ]

前面的列表也可以用 Tuples 生成，不过速度慢了近一倍：

Times @@@ Tuples[Range[100, 999], 2] // ReverseSort //
  SelectFirst[PalindromeQ]

P5: Smallest multiple

最小公倍数．

LCM @@ Range[20]

P6: Sum square difference

平方和与和平方的差．

Sum[i, {i, 100}]^2 - Sum[i^2, {i, 100}]

如果先把求和公式展开，再带入数值，则会快一些：

Expand[Sum[n, {n, N}]^2 - Sum[n^2, {n, N}]] /. N -> 100

不用 Expand 函数会更快一些，时间能缩减到第一种算法的一半：

Sum[n, {n, N}]^2 - Sum[n^2, {n, N}] /. N -> 100

P7: 10001st prime

第 10001 个素数．

Prime[10001]

P8: Largest product in a series

连续 13 位数字乘积的最大值．

t = 7316717653133062491922511967442657474235534919493496983520312774506326239578318016984801869478851843858615607891129494954595017379583319528532088055111254069874715852386305071569329096329522744304355766896648950445244523161731856403098711121722383113622298934233803081353362766142828064444866452387493035890729629049156044077239071381051585930796086670172427121883998797908792274921901699720888093776657273330010533678812202354218097512545405947522435258490771167055601360483958644670632441572215539753697817977846174064955149290862569321978468622482839722413756570560574902614079729686524145351004748216637048440319989000889524345065854122758866688116427171479924442928230863465674813919123162824586178664583591245665294765456828489128831426076900422421902267105562632111110937054421750694165896040807198403850962455444362981230987879927244284909188845801561660979191338754992005240636899125607176060588611646710940507754100225698315520005593572972571636269561882670428252483600823257530420752963450;

Times @@@ ToExpression /@
   Partition[StringSplit[IntegerString@t, ""], 13, 1] // Max

如果先把含有 0 的部分删除，速度会快一倍以上：

Times @@@ ToExpression /@ Select[! MemberQ[#, "0", 2] &] @
  Partition[StringSplit[IntegerString@t, ""], 13, 1] // Max

P9: Special Pythagorean triplet

勾股方程．

sol = Solve[{a^2 + b^2 == c^2, a + b + c == 1000, 0 < a < b < c}, {a, b, c}, Integers];
a b c /. sol

P10: Summation of primes

素数求和．

Sum[Prime[i],{i,PrimePi[2*^6]}]

直接用 Prime 函数的话，耗时 14 秒．

改用 PrimeQ 函数进行筛选的话，则只需要 0.8 秒：

Range[2*^6] // Select[PrimeQ] // Total

P11: Largest product in a grid

方阵中横、竖、斜连续四个数字乘积的最大值．

t = "08 02 22 97 38 15 00 40 00 75 04 05 07 78 52 12 50 77 91 08 \n
  49 49 99 40 17 81 18 57 60 87 17 40 98 43 69 48 04 56 62 00 \n
  81 49 31 73 55 79 14 29 93 71 40 67 53 88 30 03 49 13 36 65 \n
  52 70 95 23 04 60 11 42 69 24 68 56 01 32 56 71 37 02 36 91 \n
  22 31 16 71 51 67 63 89 41 92 36 54 22 40 40 28 66 33 13 80 \n
  24 47 32 60 99 03 45 02 44 75 33 53 78 36 84 20 35 17 12 50 \n
  32 98 81 28 64 23 67 10 26 38 40 67 59 54 70 66 18 38 64 70 \n
  67 26 20 68 02 62 12 20 95 63 94 39 63 08 40 91 66 49 94 21 \n
  24 55 58 05 66 73 99 26 97 17 78 78 96 83 14 88 34 89 63 72 \n
  21 36 23 09 75 00 76 44 20 45 35 14 00 61 33 97 34 31 33 95 \n
  78 17 53 28 22 75 31 67 15 94 03 80 04 62 16 14 09 53 56 92 \n
  16 39 05 42 96 35 31 47 55 58 88 24 00 17 54 24 36 29 85 57 \n
  86 56 00 48 35 71 89 07 05 44 44 37 44 60 21 58 51 54 17 58 \n
  19 80 81 68 05 94 47 69 28 73 92 13 86 52 17 77 04 89 55 40 \n
  04 52 08 83 97 35 99 16 07 97 57 32 16 26 26 79 33 27 98 66 \n
  88 36 68 87 57 62 20 72 03 46 33 67 46 55 12 32 63 93 53 69 \n
  04 42 16 73 38 25 39 11 24 94 72 18 08 46 29 32 40 62 76 36 \n
  20 69 36 41 72 30 23 88 34 62 99 69 82 67 59 85 74 04 36 16 \n
  20 73 35 29 78 31 90 01 74 31 49 71 48 86 81 16 23 57 05 54 \n
  01 70 54 71 83 51 54 69 16 92 33 48 61 43 52 01 89 19 67 48";

te = StringSplit[t, "\n"] // StringSplit // Map[ToExpression, #, 2] &;
Table[te[[i,j]]*te[[i,j+1]]*te[[i,j+2]]*te[[i,j+3]], {i,20}, {j,17}]~Join~
  Table[te[[i,j]]*te[[i+1,j]]*te[[i+2,j]]*te[[i+3,j]], {i,17}, {j,20}] ~Join~
  Table[te[[i,j]]*te[[i+1,j+1]]*te[[i+2,j+2]]*te[[i+3,j+3]], {i,17}, {j,17}]~Join~
  Table[te[[i,j]]*te[[i+1,j-1]]*te[[i+2,j-2]]*te[[i+3,j-3]], {i,17}, {j,4,20}] // Max

P12: Highly divisible triangular number

三角数的因数个数．

Module[{n = 1, t},
 While[DivisorSigma[0, t = n*(n + 1)/2] <= 500, n++];
 t]

受第2题启发，也可以用 NestWhile 函数来写，不过速度稍微慢一点。

NestWhile[(# + 1) &, 1, DivisorSigma[0, # (# + 1)/2] <= 500 &] //
  # (# + 1)/2 &

P13: Large sum

一百个整数相乘．

t = "37107287533902102798797998220837590246510135740250 \n  46376937677490009712648124896970078050417018260538 \n  74324986199524741059474233309513058123726617309629 \n  91942213363574161572522430563301811072406154908250 \n  23067588207539346171171980310421047513778063246676 \n  89261670696623633820136378418383684178734361726757 \n  28112879812849979408065481931592621691275889832738 \n  44274228917432520321923589422876796487670272189318 \n  47451445736001306439091167216856844588711603153276 \n  70386486105843025439939619828917593665686757934951 \n  62176457141856560629502157223196586755079324193331 \n  64906352462741904929101432445813822663347944758178 \n  92575867718337217661963751590579239728245598838407 \n  58203565325359399008402633568948830189458628227828 \n  80181199384826282014278194139940567587151170094390 \n  35398664372827112653829987240784473053190104293586 \n 86515506006295864861532075273371959191420517255829 \n  71693888707715466499115593487603532921714970056938 \n  54370070576826684624621495650076471787294438377604 \n  53282654108756828443191190634694037855217779295145 \n  36123272525000296071075082563815656710885258350721 \n  45876576172410976447339110607218265236877223636045 \n 17423706905851860660448207621209813287860733969412 \n  81142660418086830619328460811191061556940512689692 \n  51934325451728388641918047049293215058642563049483 \n  62467221648435076201727918039944693004732956340691 \n  15732444386908125794514089057706229429197107928209 \n  55037687525678773091862540744969844508330393682126 \n  18336384825330154686196124348767681297534375946515 \n  80386287592878490201521685554828717201219257766954 \n  78182833757993103614740356856449095527097864797581 \n  16726320100436897842553539920931837441497806860984 \n  48403098129077791799088218795327364475675590848030 \n  87086987551392711854517078544161852424320693150332 \n  59959406895756536782107074926966537676326235447210 \n  69793950679652694742597709739166693763042633987085 \n  41052684708299085211399427365734116182760315001271 \n  65378607361501080857009149939512557028198746004375 \n  35829035317434717326932123578154982629742552737307 \n  94953759765105305946966067683156574377167401875275 \n  88902802571733229619176668713819931811048770190271 \n  25267680276078003013678680992525463401061632866526 \n  36270218540497705585629946580636237993140746255962 \n  24074486908231174977792365466257246923322810917141 \n  91430288197103288597806669760892938638285025333403 \n  34413065578016127815921815005561868836468420090470 \n  23053081172816430487623791969842487255036638784583 \n  11487696932154902810424020138335124462181441773470 \n  63783299490636259666498587618221225225512486764533 \n  67720186971698544312419572409913959008952310058822 \n  95548255300263520781532296796249481641953868218774 \n  76085327132285723110424803456124867697064507995236 \n  37774242535411291684276865538926205024910326572967 \n  23701913275725675285653248258265463092207058596522 \n  29798860272258331913126375147341994889534765745501 \n  18495701454879288984856827726077713721403798879715 \n  38298203783031473527721580348144513491373226651381 \n  34829543829199918180278916522431027392251122869539 \n  40957953066405232632538044100059654939159879593635 \n  29746152185502371307642255121183693803580388584903 \n  41698116222072977186158236678424689157993532961922 \n  62467957194401269043877107275048102390895523597457 \n  23189706772547915061505504953922979530901129967519 \n  86188088225875314529584099251203829009407770775672 \n  11306739708304724483816533873502340845647058077308 \n  82959174767140363198008187129011875491310547126581 \n  97623331044818386269515456334926366572897563400500 \n  42846280183517070527831839425882145521227251250327 \n  55121603546981200581762165212827652751691296897789 \n  32238195734329339946437501907836945765883352399886 \n  75506164965184775180738168837861091527357929701337 \n  62177842752192623401942399639168044983993173312731 \n  32924185707147349566916674687634660915035914677504 \n  99518671430235219628894890102423325116913619626622 \n  73267460800591547471830798392868535206946944540724 \n  76841822524674417161514036427982273348055556214818 \n  97142617910342598647204516893989422179826088076852 \n  87783646182799346313767754307809363333018982642090 \n  10848802521674670883215120185883543223812876952786 \n  71329612474782464538636993009049310363619763878039 \n  62184073572399794223406235393808339651327408011116 \n  66627891981488087797941876876144230030984490851411 \n  60661826293682836764744779239180335110989069790714 \n  85786944089552990653640447425576083659976645795096 \n  66024396409905389607120198219976047599490197230297 \n  64913982680032973156037120041377903785566085089252 \n  16730939319872750275468906903707539413042652315011 \n  94809377245048795150954100921645863754710598436791 \n  78639167021187492431995700641917969777599028300699 \n  15368713711936614952811305876380278410754449733078 \n  40789923115535562561142322423255033685442488917353 \n  44889911501440648020369068063960672322193204149535 \n  41503128880339536053299340368006977710650566631954 \n  81234880673210146739058568557934581403627822703280 \n  82616570773948327592232845941706525094512325230608 \n  22918802058777319719839450180888072429661980811197 \n  77158542502016545090413245809786882778948721859617 \n  72107838435069186155435662884062257473692284509516 \n  20849603980134001723930671666823555245252804609722 \n  53503534226472524250874054075591789781264330331690";

ToExpression /@ StringSplit[t, "\n"] // Total // IntegerString //
  StringTake[#, 10] &

P14: Longest Collatz sequence

Collatz 猜想（3n+1猜想），找到一百万以下步骤最多的数字．

Table[NestWhileList[If[EvenQ[#], #/2, 3 # + 1] &, n, # != 1 &] //
    Length, {n, 1000000 - 1}] // Ordering[#, -1] &

未经优化，计算耗时3分半．

P15: Lattice paths

格点路径．

Binomial[20 + 20, 20]

P16: Power digit sum

计算各位数字和．

DigitCount[2^1000]*{1, 2, 3, 4, 5, 6, 7, 8, 9, 0} // Total

直接使用 IntegerDigits 函数会更快一些．

Total@IntegerDigits[2^1000]

P17: Number letter counts

数字转换为英文，计算单词总长度．

(** 获取对应数字的英文 **)
NumberToEnglish[n_ /; 0 < n <= 1000] :=
 Module[{hn, tn, table1, table2, a, b, c, d, res},
  table1 = {"one", "two", "three", "four", "five",
    "six", "seven", "eight", "nine", "ten",
    "eleven", "twelve", "thirteen", "fourteen", "fifteen",
    "sixteen", "eventeen", "eighteen", "nignteen"};
  table2 = {"ten", "twenty", "thirty", "forty", "fifty",
    "sixty", "seventy", "eighty", "nighty"};
  (** 百位 **)
  {a, b} = QuotientRemainder[n, 100];
  hn = If[a != 0, table1~Part~a <> "hundred", ""] <>
    If[a != 0 && b != 0, "and", ""];
  (** 十位、个位 **)
  {c, d} = QuotientRemainder[b, 10];
  tn = If[c < 2, If[b == 0, "", table1~Part~b],
    table2~Part~c <> If[d == 0, "", table1~Part~d]];
  (** 千位 **)
  If[n == 1000, "onethousand", hn <> tn]]；

Sum[StringLength@NumberToEnglish[i], {i, 1000}]

Mathematica 有 IntegerName 函数可以用，但是没有 and 连接词．

P18: Maximum path sum I

最大路径和．

t = "75\n
95 64\n
17 47 82\n
18 35 87 10\n
20 04 82 47 65\n
19 01 23 75 03 34\n
88 02 77 73 07 63 67\n
99 65 04 28 06 16 70 92\n
41 41 26 56 83 40 80 70 33\n
41 48 72 33 47 32 37 16 94 29\n
53 71 44 65 25 43 91 52 97 51 14\n
70 11 33 28 77 73 17 78 39 68 17 57\n
91 71 52 38 17 14 91 43 58 50 27 29 48\n
63 66 04 68 89 53 67 30 73 16 69 87 40 31\n
04 62 98 27 23 09 70 98 73 93 38 53 60 04 23";

(** 数据导入数组 **)
te = ImportString[t, "Table"];

(** 每次更新1行，获得截至该行各元素的最长路径 **)
For[i = 2, i <= Length@te, i++,
 te[[i]] += MapThread[Max]@{Prepend[te[[i - 1]], 0], Append[te[[i - 1]], 0]}]
te[[-1]] // Max

P19: Counting Sundays

计算每月1日星期日的个数．

(** 获得该月的天数 **)
DayInMonth[year_, month_] := If[month > 0,
   If[month != 2, {31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31}~Part~month,
    If[Mod[year, 4] == 0 && (Mod[year, 400] == 0 || Mod[year, 100] != 0), 29, 28]],
   DayInMonth[year - 1, month + 12]];

(** 建立星期表，已知1900年1月1日为星期一 **)
days = Table[{i, j, DayInMonth[i, j - 1]}, {i, 1900, 2000}, {j, 1, 12}];
days[[1, 1, 3]] = 1;

(** 获取相当于1900年1月1日的天数来计算星期几 **)
days[[All, All, 3]] =
  days[[All, All, 3]] // Flatten // FoldList[Plus] // Mod[#, 7] & //
   Partition[#, 12] &;
days[[2 ;; Length@days, All, 3]] // Flatten //
  Select[#, # == 0 &] & // Length

开始手写了判断星期的矩阵，后来发现有内置的 DayName 函数可以用：

(** 直接使用内置 DayName 函数 **)
DayName /@ Flatten[#, 1]& @ Table[{i, j, 1}, {i, 1901, 2000}, {j, 1, 12}] //
  Select[#, #==Sunday&]& // Length

P20: Factorial digit sum

阶乘的数字和．

IntegerDigits@Factorial[100] // Total

P21: Amicable numbers

相亲数．

d[n_] := DivisorSigma[1, n] - n;

(** 注意排除掉 d(1)=0, d(0) 不存在，以及 d(n)=n 的情况 **)
AmicableQ[n_] := If[d[n] != 0 && d[n] != n && d[d[n]] == n, True, False];

AmicableQ /@ Range[10000 - 1] // Position[True] // Flatten // Total

P22: Names scores

姓名得分．

names = Import["attachments/p022_names.txt", "CSV"] // Flatten // Sort;
score[name_] := LetterNumber /@ StringSplit[name, ""] // Total
Sum[score[names[[i]]]*i, {i, Length@names}] // Timing

P23: Non-abundant sums

不能写成两个盈数之和的数．

AbundantQ[n_] := DivisorSigma[1, n] > n + n;
With[{te = AbundantQ /@ Range[28123] // Position[True]},
  28123*28124/2 -
    (Table[te[[i]] + te[[j]], {i, Length@te}, {j, i}] //
      Flatten // Select[#, #=28123&]& // DeleteDuplicates // Total)]

用时1分钟．

P24: Lexicographic permutations

字典序排列．

Permutations@{0, 1, 2, 3, 4, 5, 6, 7, 8, 9}~Part~1000000 // FromDigits

P25: 1000-digit Fibonacci number

一千位的 Fibonacci 数．

Module[{n=1},
  While[IntegerLength@Fibonacci@n < 1000, n++];
  n]

NestWhile[(# + 1) &, 1, IntegerLength@Fibonacci@# < 1000 &]

两者用时相同．