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}\).
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) \]
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) '<)))
DIVISORSThen 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))))
ALPHANow, 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)))
BETALet me test:
(loop for i from 1 to 6 collect (beta i 2))
(1 6 16 37 63 113)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))))))
UPSILONLet 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 consedIt 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 secondsThe original Euler Project Question 401 is for \(n=10^{15}\) and for \(s=2\), but they require the result calculated modulo \(10^9\). The reduced calculation for \(n=10^9\) took about 180 seconds for lisp and about 6 seconds for julia. Assuming the complexity is linear (though I guess it is with a large multiplier coefficient), the real problem is out of reach for this implementation. I will have to come back to this problem later.