loop function is taking too long

Question

I'm trying to do a function who implements a sum of n cubes:

1^3 + 2^3 + 3^3 + ... + n^3 = sum

My function should receive a sum and return a n or -1 if n doesn't exists.

Some examples:

(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

After some work I made these two functions:

; 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))))))

These functions are working properly but is taking too long to evaluate operations with BigNumbers, as example:

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

I'm asking this question to raise some advices of how can I make this function faster.

Answer 1

This is all about algebra, and little to do with Clojure or programming. Since this site does not support mathematical typography, let's express it in Clojure.

Define

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

and

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

(or

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

)

Then the question is, given n, to find i such that (= (sigma-1-to-n #(* % % %) i) n).

The key to doing this quickly is Faulhaber's formula for cubes. It tells us that the following are equal, for any natural number i:

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

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

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

So, to be the sum of cubes, the number

• must be a perfect square
• whose square root is the sum of the first so many numbers.

To find out whether a whole number is a perfect square, we take its approximate floating-point square root, and see whether squaring the nearest integer recovers our whole number:

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

This returns nil if the argument is not a perfect square.

Now that we have the square root, we have to determine whether it is the sum of a range of natural numbers: in ordinary algebra, is it (j (j + 1)) / 2, for some natural number j.

We can use a similar trick to answer this question directly.

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

So the following function returns the number of successive numbers that add up to the argument, if there is one:

(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)))

We can combine these to do what you want:

(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

(Notice that the time is about twenty times faster than before, probably because HotSpot has got to work on find-n, which has been thoroughly exercised by the appended testing)

This is obviously a lot faster than the original.

---

Caveat

I was concerned that the above procedure might produce false negatives (it will never produce a false positive) on account of the finite precision of floating point. However, testing suggests that the procedure is unbreakable for the sort of number the question uses.

---

A Java double has 52 bits of precision, roughly 15.6 decimal places. The concern is that with numbers much bigger than this, the procedure may miss the exact integer solution, as the rounding can only be as accurate as the floating point number that it starts with.

However, the procedure solves the example of a 31 digit integer correctly. And testing with many (ten million!) similar numbers produces not one failure.

---

To test the solution, we generate a (lazy) sequence of [limit cube-sum] pairs:

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

For example,

(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])

Now we

• start with the given example,
• try the next ten million cases and
• remove the ones that work.

So

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

They all work. No failures. Just to make sure our test is valid:

(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])

All ought to fail, and they do.