An Unsuccessful Attempt for Solving Euler Project #401

Description of the problem

Euler Project, Problem 401 is a number theory question. It asks the reader to calculate the sum \[ \sum_{m=1}^n\sum_{q\|m} q^2 \] modulo \(10^9\) for \(n=10^{15}\).

Mathematical Reduction

The function we would like to compute is \[ \Upsilon(n,s) = \sum_{m=1}^n \sum_{q\|m} q^s \] for \(s=2\). I will define the following auxiliary function \[ \nu(q,n) = \begin{cases} 1 & \text{ if } q\|n\\ 0 & \text{ otherwise} \end{cases} \] Then I can rewrite \(\Upsilon(n,s)\) as follows: \[ \Upsilon(n,s) = \sum_{m=1}^n\sum_{q=1}^m \nu(q,m) q^s \] Now we can switch the summations \[ \Upsilon(n,s) = \sum_{q=1}^n \sum_{m=a}^n \nu(q,m) q^s = \sum_{q=1}^n q^s \left(\sum_{m=a}^n \nu(q,m)\right) \] Notice that, the inner sum terms are \(m=q,q+1,\ldots,n\) but the function we sum is non-zero only if \(m\) is a multiple of \(q\). This means the sum \(\sum_{m=q}^n \nu(q,m) = \left\lfloor\frac{n}{a}\right\rfloor\). Therefore, \[ \Upsilon(n,s) = \sum_{q=1}^n q^s \left\lfloor\frac{n}{q}\right\rfloor = \sum_{q=1}^n q^{s-1} n - \sum_{q=1}^n q^{s-1} (n \text{ mod } q) \] For \(s=2\) this can be written as \[ \Upsilon(n,2) = \frac{n^2(n+1)}{2} - \sum_{q=1}^n q (n \text{ mod } q) \]

Implementation without the reduction

First, I need a function which returns the list of divisors of a number

(defun divisors (n)
   (let* ((a (remove-if 'null 
                        (loop for i from 1 to (floor (sqrt n))
                              collect (if (zerop (mod n i)) i))))
          (b (mapcar (lambda (x) (/ n x)) a)))
      (sort (union a b) '<)))

DIVISORS

Then I need the small-case sigma (which I called alpha as lisp is not case sensitive) sigma function:

(defun alpha (n s) (reduce '+ (mapcar (lambda (x) (expt x s)) (divisors n))))

ALPHA

Now, the summotary function of the function I defined above:

(defun beta (n s) 
   (do ((i 1 (1+ i))
        (res 0 (incf res (alpha i s))))
       ((> i n) res)))

BETA

Let me test:

(loop for i from 1 to 6 collect (beta i 2))

(1 6 16 37 63 113)

Implementation with the reduction

Implementation of the function is very straight-forward:

(defun upsilon (n)
   (let ((res 0))
      (do ((q 1 (1+ q)))
          ((> q n) (- (* n n (1+ n) (/ 1 2)) res))
        (incf res (* q (mod n q))))))

UPSILON

Let me test

(loop for i from 1 to 6 collect (upsilon i))

(1 6 16 37 63 113)

and finally I will evaluate the function for \(n=10^9\)

(time (upsilon (expt 10 9)))

400685635084923815073475174

Evaluation took:
  177.784 seconds of real time
  177.403000 seconds of total run time (176.967000 user, 0.436000 system)
  [ Run times consist of 7.116 seconds GC time, and 170.287 seconds non-GC time. ]
  99.79% CPU
  496,545,276,339 processor cycles
  71,378,930,608 bytes consed

It took about three minutes to calculate the result. By the way, the same calculation implemented in Julia took about 6 seconds :)

 function euler401(n)
   res:: Uint128 = n^2*(n+1)/2
   for(q in 1:n)
      res = res - q*(n%q)
   end
   return(res)
 end

 tic(); euler401(10^9); toc()

 elapsed time: 5.98228752 seconds

Discussion and restrictions