数学的乐趣:Scala和Project Euler

缘起

刚用Scala做了一个像样的项目,对Scala热情满满,无事时用它做了几题Project Euler, 发现Scala非常适合这个,写一点心得,希望抛砖引玉,能有更多的人喜欢上Scala这门语言。

欧拉工程

从第12题说起:

Problem 12 :Highly divisible triangular number

The sequence of triangle numbers is generated by adding the natural numbers. So the 7th triangle number would be 1 + 2 + 3 + 4 + 5 + 6 + 7 = 28. The first ten terms would be:

1, 3, 6, 10, 15, 21, 28, 36, 45, 55, …

Let us list the factors of the first seven triangle numbers:

 1: 1

 3: 1,3

 6: 1,2,3,6

10: 1,2,5,10

15: 1,3,5,15

21: 1,3,7,21

28: 1,2,4,7,14,28

We can see that 28 is the first triangle number to have over five divisors.

What is the value of the first triangle number to have over five hundred divisors?

快速实现

根据题目描述,所谓triangle数,是指第n个数为1到n相加得出的数。需要求出第一个有500个因子的triangle数。所以解题需要两步:

  1. 找出所有的triangle数。
  2. 计算每个triangle的因子,找出第一个因子数大于500的triangle数。

简单!

triangular数的定义

triangular数按定义最直接的实现。

def triangulars(n: Int): Int = n match {
 case 1 => 1
 case num => triangulars(num-1) + num
}

如果是第一个数返回1,否则返回前一个triangular数和n之和,递归实现。

检查因子数:

def divisisorsNum(n: Int) = {
  (1 to n).filter( n % _ == 0).size 
}

n和所有小于n的数取模,并统计能被整除(取模结果为0)的数的个数,即为n的因子数。

获取最小有n个因子的triangle 数:

def euler12(n: Int) = {
  Stream.from(1).map(triangulars).find(divisisorsNum(_) > n)
}

Stream.from(1)构建一个从1开始的整数列表,像整数的定义一样,这个列表是无穷大的,但其中的值在被引用之前并不存在于内存中,而是通过call by name实现惰性求值。可以理解为该列表是由一个head元素+动态计算出后续元素的函数组成。

通过对整数列表执行triangulars函数的map,生成一个triangular数组成的列表,同样,这个列表也是无穷大和惰性求值的。

最后查找其中divisisorsNum也就是因子数大于n的第一个triangular数。

运行:

scala> euler12(5)
res4: Int = 28

scala> euler12(50)
res5: Int = 25200

scala> euler12(500)

因子数较小的测试都没有问题,但在求本题的答案时(500个因子)挂死,原因是我们的程序效率太低,CPU算不过来了。

改进

尾递归

首先triangulars函数中,我们用到了递归,但不是尾递归,如若递归层数太深,会发生栈溢出。

scala> triangulars(Int.MaxValue)
java.lang.StackOverflowError
  at .triangulars(<console>:12)
  at .triangulars(<console>:12)
  at .triangulars(<console>:12)

改成尾递归

@tailrec
def triangulars(n: Long, acc: Long = 1): Long = n match {
  case 1 => acc
  case num => triangulars(n-1, acc + n)
}
scala> triangulars(Int.MaxValue)
res2: Long = 2305843008139952128

但改完后求本题答案还是挂死(其实跑几个小时还是能算出答案的),仔细观察,euler12获取triangular数列表的效率非常低,反复的对每个整数计算对应的triangular数,而这些计算有很大一部分是重复的,如triangulars(100),其实就是triangulars(99)+100,不需要从1开始重新计算。

所以需要

更高效的构建triangular数列表

稍加推导,除第一个元素,第n个triangular数是前一个triangular数和n之和,直接用这个定义构成triangular数列表。

def triangulars(a: Int=1, n: Int=2): Stream[Int] = a #:: triangulars(a+n, n+1)

新的triangulars函数直接生成一个Stream,第一个参数作为递归的初始值,第二个参数标记index,每次递归时加1。

scala> triangulars().take(10).force
res6: scala.collection.immutable.Stream[Int] = Stream(1, 3, 6, 10, 15, 21, 28, 36, 45, 55)

作为测试,从这个列表中取出10个数,并强制求值,查看结果。

更高效的检查因子数

前面的divisisorsNum函数检查因子数的方法是把n和所有小于n的整数都取模,过滤出其中能被除尽的数,这样效率显然太低了,尤其是取模运算很耗CPU。

更高效的方法是:只检查小于根号n的数,并将结果乘以2(每个能整除的因子必然还有一个大于根号n小于n的同伴)。

def divisisorNumbers(n: Int) = {
  Stream.from(1).takeWhile(i => i*i <= n)
    .foldLeft(0)((acc, i) => if (n % i == 0) acc + 2 else acc)
}

这里用到了foldLeft,其作用是把Stream中所有的元素过一遍指定的函数,并将结果累加返回。

foldLeft(0)((acc, i) => if (n % i == 0) acc + 2 else acc) 第一个参数0是累加器的初始值,acc是累加器,i为Stream中的各个元素。

改进后的版本:

def euler12(n: Int) = {
  triangulars().find(divisisorNumbers(_) > n)
}
scala> euler12(500)
res21: Option[Int] = Some(76576500)

统计一下时间,花了4秒多得出结果。

scala> time {euler12(500)}
4 seconds 433 microseconds 155 milliseconds
res24: Option[Int] = Some(76576500)

已经到了可用级别了。但是:

能否更快

数论的知识,任何一个整数都能唯一的表示为几个素数的乘积。

其中n为任意整数,p1-ps为素数。

如,4200可以表示成:

所以4200有 (3+1)(1+1)(2+1)(1+1) = 48个因子。

而triangle数可以表示为:

所以triangles(n)的因子数为:n的因子数乘以n+1的因子数,其中n或者n+1中偶数那个因子需要减掉一个因子(因为要除了2),而n和n+1中必然有一个偶数。

代码实现

/*
* 计算n的因子个数,排除一次2因子
* */
def divisisorNumbers(n: Int, sum: Int, primes: Stream[Int])
  : Int = primes match {
  case head #:: tail if head > n => sum
  case prime #:: tail =>
    var number = n
    var count = 0
    while (number % prime == 0) {
      number = number / prime
      count = count + 1
    }
    if (prime == 2 && count > 0) {
      count = count - 1
    }
    divisisorNumbers(number, sum * (count + 1), tail)
}
def euler12(n: Int) = {
  val primeStreamInt  = Stream.from(1)
    .filter(BigInt(_).isProbablePrime(15))
  Stream.from(1)
    .find(i => divisisorNumbers(i,1, primeStreamInt) *
          divisisorNumbers(i + 1, 1,primeStreamInt) > n)
    .map(i => i * (i + 1) / 2)
}

代码中用到的primeStreamInt是利用JDK的isProbablePrime函数过滤出的素数列表,用于查找因子,算是一个小cheat。

统计时间:性能提高到了毫秒级,但是我还是更喜欢前面那个4秒多的实现:)

scala> time(euler12(500))
0 seconds 156 microseconds 568 milliseconds
res6: Option[Int] = Some(76576500)

结语

Scala是一门很有趣的语言,同时具有函数式语言的灵敏和Java生态系统的严谨。如果你还有学习的动力,把它作为你的下一门语言吧。