; the divisor function sigma

(define sort #f)
(define merge #f)
(let ()
  (define dosort
    (lambda (pred? ls n)
      (if (= n 1)
          (list (car ls))
          (let ((i (quotient n 2)))
            (domerge pred?
                     (dosort pred? ls i)
                     (dosort pred? (list-tail ls i) (- n i)))))))
  (define domerge
    (lambda (pred? l1 l2)
      (cond
        ((null? l1) l2)
        ((null? l2) l1)
        ((pred? (car l2) (car l1))
         (cons (car l2) (domerge pred? l1 (cdr l2))))
        (else (cons (car l1) (domerge pred? (cdr l1) l2))))))
  (set! sort
    (lambda (pred? l)
      (if (null? l) l (dosort pred? l (length l)))))
  (set! merge
    (lambda (pred? l1 l2)
      (domerge pred? l1 l2))))

(define (uniq-c eql? xs)
  (if (null? xs) xs
    (let loop ((xs (cdr xs)) (prev (car xs)) (k 1) (result '()))
      (cond ((null? xs) (reverse (cons (cons prev k) result)))
            ((eql? (car xs) prev) (loop (cdr xs) prev (+ k 1) result))
            (else (loop (cdr xs) (car xs) 1 (cons (cons prev k) result)))))))

(define (factors n)
  (let twos ((n n) (fs '()))
    (if (even? n) (twos (/ n 2) (cons 2 fs))
      (if (= n 1) fs
        (let odds ((n n) (d 3) (fs fs))
          (cond ((< n (* d d))
                  (reverse (cons n fs)))
                ((zero? (modulo n d))
                  (odds (/ n d) d (cons d fs)))
                (else (odds n (+ d 2) fs))))))))

(define (divisor x n)
  (define (prod xs) (apply * xs))
  (if (= n 1) 1
    (let ((fs (uniq-c = (sort < (factors n)))))
      (if (zero? x)
          (prod (map add1 (map cdr fs)))
          (prod (map (lambda (p a)
                       (/ (- (expt p (* (+ a 1) x)) 1)
                          (- (expt p x) 1)))
                     (map car fs) (map cdr fs)))))))

(define (summatory-divisor x n)
  (do ((d 1 (+ d 1))
       (s 0 (+ s (* (expt d x) (quotient n d)))))
      ((< n d) s)))

(display (divisor 0 60)) (newline)
(display (divisor 1 60)) (newline)
(display (divisor 2 60)) (newline)

(display (summatory-divisor 0 60)) (newline)
(display (summatory-divisor 1 60)) (newline)
(display (summatory-divisor 2 60)) (newline)

(define (euler-401 n)
  (define (p x) (quotient (* x (+ x 1) (+ x x 1)) 6))
  (let loop ((i 1) (s 0))
    (if (< n i) s
      (let* ((m (quotient n i)) (j (+ (quotient n m) 1)))
        (loop j (+ s (* m (- (p (- j 1)) (p (- i 1))))))))))

(time (display (summatory-divisor 2 1000000)) (newline))
(time (display (euler-401 1000000)) (newline))

12
168
5460
261
3014
89286
400686363385965077
cpu time: 510 real time: 4006 gc time: 210
400686363385965077
cpu time: 0 real time: 33 gc time: 0