Problem 12

Highly divisible triangular number

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

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?

有很多约数的三角形数

三角形数是通过累加自然数所得到的数。例如，第 7 个三角形数是 1+2+3+4+5+6+7=28。前十个三角形数分别是：

列举出前七个三角形数的所有约数：

 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

可以看出，28 是第一个约数数量超过五的三角形数。

第一个约数数量超过五百的三角形数是多少？

首先我们需要三角形数的计算公式：

(defun eu12-tri (n)
  (/ (* n (1+ n)) 2))

接着我们可以考虑使用最直接的方法求出因数：

(defun eu12-factors (n)
  (let ((bound (/ n 2)))
    (named-let f ((i 2) (sum 2))
      (cond
       ((> i bound) sum)
       (t (if (zerop (% n i)) (f (+ i 1) (1+ sum))
	    (f (+ i 1) sum)))))))

(my-time
 (named-let f ((i 1))
   (let* ((n (eu12-tri i))
	  (divs (eu12-factors n)))
     (if (> divs 500) n
       (f (1+ i)))))
 )
;;0.952765941619873s
;;0.13864803314208984s
;;0.08629703521728516s
;;0.061781883239746094s
;;0.006681203842163086s
;;0.006932973861694336s
;;0.0012028217315673828s
;;0.00011301040649414062s
;;0.00012087821960449219s
;;2.09808349609375e-05s

上面的时间是我将除数个数从 5 增加至 50 分别用的时间，可见到 50 的时候就已经需要接近一秒了，使用上面的方式获取除数个数大于 500 的三角形数是不现实的。我们需要换一种思路。一个小优化是根据除数得到它的对偶，比如 28 = 4*7，我们可由 4 得到 7 而不用继续算下去：

(my-time
 (named-let f ((i 0) (num 0))
   (let* ((i1 (1+ i))
	  (num1 (+ num i1))
	  (bound (floor (sqrt num1)))
	  (sum 0))
     (cl-loop for n from 1 to bound
	      if (zerop (% num1 n))
	      do (setq sum (+ sum 2)))
     (if (> sum 500) num1
       (f i1 num1))))
 )
;;7.196847915649414s
;; 76576500

如果我们将数字进行质因数分解，然后再通过这些质因数组合获取因数的话能快上不少，这是因为质因数分解比求数字因数要快很多。在获取质因数后，我们可以将不同质因数的重数加一相乘，这样就得到了因数个数，比如下面的例子：

之所以能够这样做，是因为一个数的不同质因数相乘的结果必然是唯一的，一个数的相同质因数不同个数相乘得到的结果也是唯一的，这样就成了一个组合问题：

之所以要添加 1 是因为某一组合可能不包含某个质因数和它们组合相乘的结果。

(defun eu12-factors-2 (n)
  (let ((numof2 1))
    (while (zerop (% n 2))
      (cl-incf numof2)
      (setq n (/ n 2)))
    (named-let f ((i 3) (sum 1))
      (cond
       ((> i n) (* sum numof2))
       (t
	(if (zerop (% n i))
	    (let ((cnt (1+ (cl-loop
			    while (zerop (% n i))
			    do (setq n (/ n i))
			    sum 1))))
	      (f (+ i 2) (* sum cnt)))
	  (f (+ i 2) sum)))))))

(my-time
 (named-let f ((i 1))
   (let* ((n (eu12-tri i))
	  (divs (eu12-factors-2 n)))
     (if (> divs 500) n
       (f (1+ i)))))
 )
;;3.166064977645874s
;; 7657600

虽然这种方法能够得出结果了，但是显然它还不够快。一个比较有意思的现象是 gcd(n,n+1)=1 ，也就是说两相邻自然数是互素的。而由于 n 与 n+1 互素，那么如果 n 为偶数，n/2 与 n+1 互素，如果 n 为奇数，n 与 (n+1)/2 互素。由于互素的两个数字没有共同的质因数，我们可以做如下分解：

因为 n 或 n+1 的规模是远小于 n*(n+1)/2 的，下面的实现应该会更快：（似乎并没有）

(my-time
 (named-let f ((i 2))
   (let* ((flg (zerop (% i 2)))
	  (n (if flg (/ i 2) (/ (1+ i) 2)))
	  (m (if flg (1+ i) i)))
     (let ((D1 (eu12-factors-2 n))
	   (D2 (eu12-factors-2 m)))
       (if (> (* D1 D2) 500) (* m n)
	 (f (+ i 1))))))
 )
;;3.5967118740081787s

这个实现还能继续优化，比如我们可以把已有的 D 保存下来供后面使用。不过我有点怀疑在规模比较小时，朴素方法反而会更快…

(defun eu12-factors-3 (n)
  (let ((bound (floor (sqrt n))))
    (cl-loop for i from 1 to bound
	     if (zerop (% n i))
	     sum 2)))

(my-time
 (named-let f ((i 1))
   (let* ((n (eu12-tri i))
	  (divs (eu12-factors-3 n)))
     (if (> divs 500) n
       (f (1+ i)))))
 )
;;7.149927854537964s

(my-time
 (named-let f ((i 2))
   (let* ((flg (zerop (% i 2)))
	  (n (if flg (/ i 2) (/ (1+ i) 2)))
	  (m (if flg (1+ i) i)))
     (let ((D1 (eu12-factors-3 n))
	   (D2 (eu12-factors-3 m)))
       (if (> (* D1 D2) 500) (* m n)
	 (f (+ i 1))))))
 )
;;0.2852509021759033s

果然如此…

既然我们已经知道了数字的范围，我们可以将素数表用于 eu12-factors-2 ，这应该能快上不少。不过也不好说，毕竟数字的规模还不是很大，这里我就不继续了。