5

I've been slowly working my way through the list of Project Euler problems, and I've come to one that I know how to solve, but it seems like I can't (given the way my solution was written).

I am using Common Lisp to do this with and my script has been running for over 24 hours (well over their one minute goal).

For the sake of conciseness, here's my solution (it's a spoiler, but only if you have one hell of a fast processor):

(defun square? (num)
  (if (integerp (sqrt num)) T))

(defun factors (num)
  (let ((l '()))
    (do ((current 1 (1+ current)))
        ((> current (/ num current)))
      (if (= 0 (mod num current))
          (if (= current (/ num current))
              (setf l (append l (list current)))
              (setf l (append l (list current (/ num current)))))))
    (sort l #'< )))

(defun o_2 (n)
  (reduce #'+ (mapcar (lambda (x) (* x x)) (factors n))))

(defun sum-divisor-squares (limit)
  (loop for i from 1 to limit when (square? (o_2 i)) summing i))

(defun euler-211 ()
 (sum-divisor-squares 64000000))

The time required to solve the problem using smaller, more friendly, test arguments seems to grow larger than exponentialy... which is a real problem.

It took:

  • 0.007 seconds to solve for 100
  • 0.107 seconds to solve for 1000
  • 2.020 seconds to solve for 10000
  • 56.61 seconds to solve for 100000
  • 1835.385 seconds to solve for 1000000
  • 24+ hours to solve for 64000000

I'm really trying to figure out which part(s) of the script is causing it to take so long. I've put some thought into memoizing the factors function, but I'm at a loss as to how to actually implement that.

For those that want to take a look at the problem itself, here it be.

Any ideas on how to make this thing go faster would be greatly appreciated.

**sorry if this is a spoiler to anyone, it's not meant to be.... but if you have the computing power to run this in a decent amount of time, more power to you.

2
  • Your timings is increasing roughly geometrically with the number. An approximation is 0.00001 * n^(1.36). This means 64M will take about 4.6 days. Dec 3, 2008 at 1:37
  • the square? is wrong. Common Lisp does not specify that (sqrt 4) is the integer 2 . It could also be the float 2.0 . Also note that (if (pred-p argument) t) is the same as just (pred- argument) for most practical purposes. Mar 11, 2009 at 8:05

7 Answers 7

14

Here's a solution, keeping in mind the spirit of [Project] Euler. [Warning: spoiler. I've tried to keep the hints slow, so that you can read only part of the answer and think on your own if you want. :)]

When you are confronted with a problem having to do with numbers, one good strategy (as you probably already know from solving 210 Project Euler problems) is to look at small examples, find a pattern, and prove it. [The last part may be optional depending on your attitude to mathematics ;-)]

In this problem, though, looking at small examples -- for n=1,2,3,4,... will probably not give you any hint. But there is another sense of "small examples" when dealing with number-theoretic problems, which you also probably know by now -- primes are the building blocks of the natural numbers, so start with the primes.

For a prime number p, its only divisors are 1 and p, so the sum of the squares of its divisors is 1+p2.
For a prime power pk, its only divisors are 1, p, p2, … pk, so the sum of the squares of its divisors is 1+p+p2+…+pk=(pk+1-1)/(p-1).
That was the simplest case: you've solved the problem for all numbers with only one prime factor.

So far nothing special. Now suppose you have a number n that has two prime factors, say n=pq. Then its factors are 1, p, q, and pq, so the sum of the squares of its divisors is 1+p2+q2+p2q2=(1+p2)(1+q2).
What about n=paqb? What is the sum of the squares of its factors?

[............................Dangerous to read below this line...................]

It is ∑0≤c≤a, 0≤d≤b(pcqd)2 = ((pa+1-1)/(p-1))((qb+1-1)/(q-1)).

That should give you the hint, both on what the answer is and how to prove it: the sum of the divisors of n is simply the product of the (answer) for each of the prime powers in its factorization, so all you need to do is to factorize 64000000 (which is very easy to do even in one's head :-)) and multiply the answer for each (=both, because the only primes are 2 and 5) of its prime powers.

That solves the Project Euler problem; now the moral to take away from it.

The more general fact here is about multiplicative functions -- functions on the natural numbers such that f(mn) = f(m)f(n) whenever gcd(m,n)=1, i.e. m and n have no prime factors in common. If you have such a function, the value of the function at a particular number is completely determined by its values at prime powers (can you prove this?)

The slightly harder fact, which you can try to prove[it's not that hard], is this: if you have a multiplicative function f [here, f(n)=n2] and you define the function F as F(n) = ∑d divides nf(d), (as the problem did here) then F(n) is also a multiplicative function.

[In fact something very beautiful is true, but don't look at it just yet, and you'll probably never need it. :-)]

4
  • To be honest, I've not done 210, I've merely done 67. This one struck my attention, so I tried it out. But thanks for the insight on how to more effectively do this one. Dec 3, 2008 at 3:22
  • The Möbius inversion formula is very, very beautiful; unfortunately, it is also very, very difficult to intuit. It was covered in my introductory number theory class but less than half of the students were able to apply it properly even when given hints. Some days, it feels like dark magic to me, too. :)
    – ephemient
    May 22, 2009 at 22:10
  • "For a prime power pk, its only divisors are 1, p, p2, … pk, so the sum of the squares of its divisors is 1+p+p2+…+pk=(pk+1-1)/(p-1)." There is a typo in above statement, should be "...so the sum of the squares of its divisors is 1+p^2+p^2^2+p^3^2+...+pk^2=(pk^(2+1)-1)/(pk-1)".
    – dz902
    Dec 15, 2016 at 14:31
  • Sorry, should be "1+p^2+p^2^2+p^3^2+...+pk^2=(p(k+1)^2-1)/(p^2-1)". But I really wonder how this is deduced...Cheers.
    – dz902
    Dec 15, 2016 at 14:56
3

I think that your algorithm is not the most efficient possible. Hint: you may be starting from the wrong side.

edit: I'd like to add that choosing 64000000 as the upper limit is likely the problem poster's way of telling you to think of something better.

edit: A few efficiency hints:

  • instead of
(setf l (append l (...)))

you can use

(push (...) l)

which destructively modifies your list by consing a new cell with your value as car and the former l as cdr, then points l to this cell. This is much faster than appending which has to traverse the list once each. If you need the list in the other order, you can nreverse it after it is complete (but that is not needed here).

  • why do you sort l?

  • you can make (> current (/ num current)) more efficient by comparing with the square root of num instead (which only needs to be computed once per num).

  • is it perhaps possible to find the factors of a number more efficiently?

And a style hint: You can put the scope of l into the do declaration:

(do ((l ())
     (current 1 (+ current 1)))
    ((> current (/ num current))
     l)
  ...)
0
2

I would attack this by doing the prime factorization of the number (for example: 300 = 2^2 * 3^1 * 5^2), which is relatively fast, especially if you generate this by sieve. From this, it's relatively simple to generate the factors by iterating i=0..2; j=0..1; k=0..2, and doing 2^i * 3^j * 5^k.

5 3 2
-----
0 0 0 = 1
0 0 1 = 2
0 0 2 = 4
0 1 0 = 3
0 1 1 = 6
0 1 2 = 12
1 0 0 = 5
1 0 1 = 10
1 0 2 = 20
1 1 0 = 15
1 1 1 = 30
1 1 2 = 60
2 0 0 = 25
2 0 1 = 50
2 0 2 = 100
2 1 0 = 75
2 1 1 = 150
2 1 2 = 300

This might not be fast enough

1
  • It's definitely fast enough, and this is a generic solution that works for many problems... this should be the accepted answer instead ;-) May 12, 2009 at 16:53
2

The clever trick you are missing is that you don't need to factor the numbers at all How many numbers from 1..N are multiples of 1? N How many numbers from 1..N are multiples of 2? N/2

The trick is to sum each number's factors in a list. For 1, add 1^2 to every number in the list. For 2, add 2^2 to every other number. For 3, add 3^2 to every 3rd number.

Don't check for divisibility at all. At the end, you do have to check whether the sum is a perfect square, and that's it. In C++, this worked in 58 seconds for me.

1

Sorry, I don't understand LISP well enough to read your answer. But my first impression is that the time cost of the brute force solution should be:

open bracket

sqrt(k) to find the divisors of k (by trial division), square each one (constant time per factor), and sum them (constant time per factor). This is σ2(k), which I will call x.

plus

not sure what the complexity of a good integer square root algorithm is, but certainly no worse than sqrt(x) (dumb trial multiplication). x might well be big-O larger than k, so I reserve judgement here, but x is obviously bounded above by k^3, because k has at most k divisors, each itself no bigger than k and hence its square no bigger than k^2. It's been so long since my maths degree that I have no idea how fast Newton-Raphson converges, but I suspect it's faster than sqrt(x), and if all else fails a binary chop is log(x).

close bracket

multiplied by n (as k ranges 1 .. n).

So if your algorithm is worse than O(n * sqrt(n^3)) = O(n ^ (5/2)), in the dumb-sqrt case, or O(n * (sqrt(n) + log(n^3)) = O(n ^ 3/2) in the clever-sqrt case, I think something has gone wrong which should be identifiable in the algorithm. At this point I'm stuck because I can't debug your LISP.

Oh, I've assumed that arithmetic is constant-time for the numbers in use. It darn well should be for numbers as small as 64 million, and the cube of that fits in a 64bit unsigned integer, barely. But even if your LISP implementation is making arithmetic worse than O(1), it shouldn't be worse than O(log n), so it won't have much affect on the complexity. Certainly won't make it super-polynomial.

This is where someone comes along and tells me just how wrong I am.

Oops, I just looked at your actual timing figures. They aren't worse than exponential. Ignoring the first and last values (because small times aren't accurately measurable and you haven't finished, respectively), multiplying n by 10 multiplies time by no more than 30-ish. 30 is about 10^1.5, which is about right for brute force as described above.

0

I think you can attack this problem with something like a prime sieve. That's only my first impression though.

1
  • Given the source, it's more likely there's some smart mathematical relation between sigma-2(n) and sigma-2 of the prime factors of n, which makes the whole thing much, much faster. Ofc you're right that calculating a prime sieve speeds up finding factors significantly. Dec 3, 2008 at 1:44
0

I've reworked the program with some notes taken from the comments here. The 'factors' function is now ever so slightly more efficient and I also had to modify the σ_(2)(n) function to accept the new output.

'factors' went from having an output like:

$ (factors 10) => (1 2 5 10)

to having one like

$ (factors 10) => ((2 5) (1 10))

Revised function looks like this:

(defun o_2 (n)
"sum of squares of divisors"
  (reduce #'+ (mapcar (lambda (x) (* x x)) (reduce #'append (factors n)))))

After the modest re-writes I did, I only saved about 7 seconds in the calculation for 100,000.

Looks like I'm going to have to get off of my ass and write a more direct approach.

1
  • You don't have to; you can do as FryGuy said: first write a prime-factors function that does (prime-factors 64000000) => (2 2 2 2 2 2 2 2 2 2 2 2 5 5 5 5 5 5), then generate the list of factors using that prime-factors function. This should help you for other Project Euler problems as well. Dec 3, 2008 at 14:51

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.