Problem 1

Multiples of 3 or 5

If we list all the natural numbers below 10 that are multiples of 3 or 5, we get 3, 5, 6 and 9. The sum of these multiples is 23.

Find the sum of all the multiples of 3 or 5 below 1000.

3 或 5 的倍数

在小于 10 的自然数中，3 或 5 的倍数有 3、5、6 和 9，这些数的和是 23。

求小于 1000 的自然数中所有 3 或 5 的倍数之和。

最简单的思路就是从 1 遍历到 999，通过取余 (%) 判断当前数字是否为 3 或 5 的倍数。1 和 2 不是 3 或 5 的倍数，我们可以从 3 开始：

(cl-loop for i from 3 below 1000
         if (or (zerop (% i 3)) (zerop (% i 5)))
         sum i)
;;=> 233168

容易看出在 999 内的 3 的倍数组成了一个等差数列，5 的倍数同理。我们可以使用等差数列求和分别得到 3 和 5 的倍数数列之和，然后减去同时是 3 和 5 的倍数的数列求和，即 15 倍数的数列求和。

;; fn = an; Sn = n*(a + an)/2
;; An = 3n; Bn = 5n n = 1, 2, 3...
;; n for A is 999/3 = 333, for B is 995/5 = 199
;; Cn = 15n; 990/15 = 66
(let ((a3 (/ (* 333 (+ 3 (* 333 3))) 2))
      (a5 (/ (* 199 (+ 5 (* 199 5))) 2))
      (a1 (/ (* 66 (+ 15 (* 66 15))) 2)))
  (+ a3 a5 (- a1)))
;;=> 233168

当然，我们也可以用打表的方法筛出所有是 3 或 5 倍数的数，然后把它们加起来：

(let ((bits (make-bool-vector 1000 nil)))
  (cl-do ((i 0 (+ i 3))) ((> i 999)) (aset bits i t))
  (cl-do ((j 0 (+ j 5))) ((> j 999)) (aset bits j t))
  (cl-loop for i from 0 below 1000
           if (aref bits i) sum i))
;;=> 233168

如果把这个问题扩展一下，我们可以得到如下描述：求在 [a, b) (a ≥ 0, b > a) 范围内至少是 a₀, a₁, … a_N (a_n > 1, N ≥ 1, a_i ≠ a_j, i ≠ j) 这些数中的一个的倍数的所有数字之和。

对于上面的取余判定方法和打表法，都可以很容易地写出通用的方法来：

(defun eu1-% (start end nums)
  (cl-loop for i from start below end
           if (cl-loop for n in nums
                       if (zerop (% i n)) return t
                       finally return nil)
           sum i))
(eu1-% 0 1000 '(3 5))
;;=> 233168

(defun eu1-bits (start end nums)
  (let* ((len (- end start))
         (bits (make-bool-vector len nil)))
    (while-let ((n (pop nums)))
      (cl-loop for i from (% start n) below len by n
               do (aset bits i t)))
    (cl-loop for i from 0 below len
             if (aref bits i)
             sum (+ start i))))
(eu1-bits 0 1000 '(3 5))
;;=> 233168

对于基于通项公式的方法，某个范围 [a, b] 内的等差数列 a_n = cn + d 求和的值可以由以下公式给出：

\[S = c \cdot \frac{(n_{max} - n_{min} + 1)(n_{min} + n_{max})}{2} + d(n_{max} - n_{min} + 1)\]

其中：

\[n_{min} = \lceil \frac{a-d}{c}\rceil, n_{max} = \lfloor \frac{b-d}{c} \rfloor\]

稍作修改我们就可以得到在 [a, b) 内的求和值：

(defun eu1-sum (a b c d)
  (let ((n-min (ceiling (/ (- a d) (+ c 0.0))))
        (n-max (let ((r (/ (- b d) c)))
      	   (if (= (* r c) b) (1- r) r))))
    (+ (/ (* c (+ n-max 1 (- n-min)) (+ n-min n-max)) 2)
       (* d (+ n-max 1 (- n-min))))))

使用容斥定理（inclusion-exclusion principle），我们可以得到所有这些满足条件的数字取并集后的和：

\[ \left| \bigcup_{i=1}^{n} A_i \right| = \sum_{i=1}^{n} |A_i| - \sum_{1 \le i < j \le n} |A_i \cap A_j| + \sum_{1 \le i < j < k \le n} |A_i \cap A_j \cap A_k| - \cdots + (-1)^{n+1} |\bigcap_{i=1}^{n}A_i| \]

以下代码可以根据列表元素列出所有可能大小的子集（也就是求幂集），具体的思路是：

\[ \mathcal{P}(S) = \begin{cases} \{\emptyset\} & \text{if } S = \emptyset \\ \mathcal{P}(S') \cup \{ \{e\} \cup X \mid X \in \mathcal{P}(S') \} & \text{if } S \neq \emptyset \end{cases} \]

(defun eu1-comb (ls)
  (if (null ls) '(())
    (let ((rest (eu1-comb (cdr ls))))
      (append rest (mapcar (lambda (x) (cons (car ls) x)) rest)))))

(eu1-comb '(3 5 4))
;;=> (nil (4) (5) (5 4) (3) (3 4) (3 5) (3 5 4))

接着，我们可以使用以下代码求解：

(defun eu1-iep (start end nums)
  (let ((combs (eu1-comb nums)))
    (cl-loop for s in combs
             for sgn = (if (evenp (length s)) -1 1)
             for val = (apply #'cl-lcm s)
             if (not (null s))
             sum (* sgn (eu1-sum start end val 0)))))
(eu1-iep 0 1000 '(3 5))
;;=> 233168

上面的代码还有一些可以优化的点，比如去掉求和函数的最后一个参数，使集合中的数字互质等等，最大的优化是改进幂集的计算。在利用容斥原理进行计算时，当涉及的集合数量 n 较小（通常 n < 20）时，位掩码是一种最高效的方法。该方法通过循环遍历整数 1 到 2ⁿ - 1 来隐式地代表指标集 {1 … n} 的所有非空子集。整数的二进制位 i 上的 1 表示指标 i 被选中；通过计算二进制数中 1 的个数 (popcount)，可以高效地确定子集大小，进而确定容斥原理中交集项所需的正负符号 (-1)^|I|-1。这种方式避免了显式生成和存储整个幂集列表。

Problem 1

1. Problem

2. Solution

3. Extension

3.1. 使用容斥定理