循环功能耗时过长

Question

我正在尝试做一个实现 n 个立方体的总和的函数：

1^3 + 2^3 + 3^3 + ... + n^3 = 总和

如果n 不存在，我的函数应该接收sum 并返回n 或-1。

一些例子：

(find-n 9)   ; should return 2 because 1^3 + 2^3 = 9
(find-n 100) ; should return 4 because 1^3 + 2^3 + 3^3 + 4^3 = 100
(find-n 10)  ; should return -1

经过一些工作，我做了这两个功能：

; aux function
(defn exp-3 [base] (apply *' (take 3 (repeat base))))

; main function
(defn find-n [m]
  (loop [sum 0
         actual-base 0]
       (if (= sum m) 
           actual-base
           (if (> sum m)
               -1
               (recur (+' sum (exp-3 (inc actual-base))) (inc actual-base))))))

这些函数工作正常，但用BigNumbers 评估操作花费的时间太长，例如：

(def sum 1025247423603083074023000250000N)
(time (find-n sum))
; => "Elapsed time: 42655.138544 msecs"
; => 45001000

我问这个问题是为了提出一些关于如何使这个功能更快的建议。

Answer 1

这都是关于代数的，与 Clojure 或编程无关。由于本站不支持数学排版，所以我们用 Clojure 来表达。

定义

(defn sigma [coll] (reduce + coll))

和

(defn sigma-1-to-n [f n]
  (sigma (map f (rest (range (inc n))))))

（或

(defn sigma-1-to-n [f n]
  (->> n inc range rest (map f) sigma))

)

那么问题是，给定n，找到i 使得(= (sigma-1-to-n #(* % % %) i) n)。

快速执行此操作的关键是 Faulhaber's formula 用于立方体。它告诉我们，对于任何自然数i，以下是相等的：

(#(*' % %) (sigma-1-to-n identity i))

(sigma-1-to-n #(* % % %) i)

(#(*' % %) (/ (*' i (inc i)) 2))

所以，作为立方体的总和，数字

• 必须是完美的正方形
• 其平方根是前这么多数字的总和。

为了确定一个整数是否是一个完美的平方，我们取它的近似浮点平方根，看看对最接近的整数求平方是否能恢复我们的整数：

(defn perfect-square-root [n]
  (let [candidate (-> n double Math/sqrt Math/round)]
    (when (= (*' candidate candidate) n)
      candidate)))

如果参数不是完美正方形，则返回 nil。

现在我们有了平方根，我们必须确定它是否是一系列自然数的和：在普通代数中，是(j (j + 1)) / 2，对于一些自然数j。

我们可以使用类似的技巧直接回答这个问题。

j (j + 1) = (j + 1/2)^2 + 1/4

因此，如果有一个，则以下函数返回与参数相加的连续数字的数量：

(defn perfect-sum-of [n]
  (let [j (-> n (*' 2)
                (- 1/4)
                double
                Math/sqrt
                (- 0.5)
                Math/round)]
    (when (= (/ (*' j (inc j)) 2) n)
      j)))

我们可以结合这些来做你想做的事：

(defn find-n [big-i]
  {:pre [(integer? big-i) ((complement neg?) big-i)]}
  (let [sqrt (perfect-square-root big-i)]
    (and sqrt (perfect-sum-of sqrt))))

(def sum 1025247423603083074023000250000N)

(time (find-n sum))
"Elapsed time: 0.043095 msecs"
=> 45001000

（请注意，时间比以前快了大约二十倍，可能是因为 HotSpot 必须在 find-n 上工作，已被附加测试彻底锻炼）

这显然比原来的要快很多。

---

警告

由于浮点的有限精度，我担心上述过程可能会产生假阴性（它永远不会产生假阳性）。然而，测试表明该过程对于问题使用的数字类型是牢不可破的。

---

Java double 有 52 位精度，大约小数点后 15.6 位。令人担忧的是，如果数字远大于此，该过程可能会错过精确的整数解，因为舍入只能与它开始的浮点数一样准确。

但是，该过程正确地解决了 31 位整数的示例。用许多（一千万！）类似的数字进行测试不会导致失败。

---

为了测试解决方案，我们生成了一个（惰性）[limit cube-sum] 对序列：

(defn generator [limit cube-sum]
  (iterate
    (fn [[l cs]]
      (let [l (inc l)
            cs (+' cs (*' l l l))]
        [limit cs]))
    [limit cube-sum]))

例如，

(take 10 (generator 0 0))
=> ([0 0] [1 1] [2 9] [3 36] [4 100] [5 225] [6 441] [7 784] [8 1296] [9 2025])

现在我们

• 从给定的例子开始，
• 尝试接下来的一千万个案例，然后
• 删除有效的。

所以

(remove (fn [[l cs]] (= (find-n cs) l)) (take 10000000 (generator 45001000 1025247423603083074023000250000N)))
=> () 

它们都有效。没有失败。只是为了确保我们的测试有效：

(remove (fn [[l cs]] (= (find-n cs) l)) (take 10 (generator 45001001 1025247423603083074023000250000N)))
=>
([45001001 1025247423603083074023000250000N]
 [45001002 1025247514734170359564546262008N]
 [45001003 1025247605865263720376770289035N]
 [45001004 1025247696996363156459942337099N]
 [45001005 1025247788127468667814332412224N]
 [45001006 1025247879258580254440210520440N]
 [45001007 1025247970389697916337846667783N]
 [45001008 1025248061520821653507510860295N]
 [45001009 1025248152651951465949473104024N]
 [45001010 1025248243783087353664003405024N])

一切都应该失败，而他们确实如此。

Answer 2

仅仅避免apply（在 CLJ 中并不是那么快）就可以让您获得 4 倍的加速：

(defn exp-3 [base]
  (*' base base base))

还有 10%：

(defn find-n [m]
  (loop [sum 0
         actual-base 0]
    (if (>= sum m)
      (if (= sum m) actual-base -1)
      (let [nb (inc actual-base)]
        (recur (+' sum (*' nb nb nb)) nb)))))

Answer 3

以下基于算法的方法依赖于一个简单的公式，即前 N 个自然数的立方和为：(N*(N+1)/2)^2

(defn sum-of-cube
  "(n*(n+1)/2)^2"
  [n]
  (let [n' (/ (*' n (inc n)) 2)]
    (*' n' n')))

(defn find-nth-cube
  [n]
  ((fn [start end prev]
     (let [avg (bigint (/ (+' start end) 2))
           cube (sum-of-cube avg)]
       (cond (== cube n) avg
             (== cube prev) -1
             (> cube n) (recur start avg cube)
             (< cube n) (recur avg end cube))))
    1 n -1))

(time (find-nth-cube 1025247423603083074023000250000N))
"Elapsed time: 0.355177 msecs"
=> 45001000N

我们想要找到数字 N，使得 1..N 立方的总和是某个数字 X。要找到这样的数字是否存在，我们可以通过应用上述公式在某个范围内对其执行二进制搜索看看公式的结果是否等于 X。这种方法有效，因为顶部的函数正在增加，因此任何值 f(n) 太大意味着我们必须寻找较小的数字 n，而任何值 @ 987654325@太小意味着我们必须寻找更大的数字n。

我们选择 0 到 X 的（大于必要的，但简单且安全的）范围。如果我们的公式应用于给定的候选数产生 X，我们将知道该数存在。如果不存在，我们继续二进制搜索直到我们找到该数字，或者直到我们尝试了两次相同的数字，这表明该数字不存在。

上限为logN，计算 1E100 (1 googol) 仅需 1 毫秒，因此对于算法方法来说非常有效。

Answer 4

您可能想要使用一些数学技巧。

(a-k)^3 + (a+k)^3 = 2a^3+(6k^2)a

所以，这样的总和：

(a-4)^3+(a-3)^3+(a-2)^3+(a-1)^3+a^3+(a+1)^3+(a+2)^3+(a+3)^3+(a+4)^3 
= 9a^3+180a

（请确认计算的正确性）。

使用这个等式，您可以跳 9（或任意 2k+1），而不是每次增加 1。每当您遇到比 n 更大的数字时，您都可以检查确切的数字。

另一种改进方法是通过一次进行一批计算，然后在函数 find-n 中使用该表，创建一个包含 n 和 sum 的表.