; approximating pi
(define (logint x)
(let ((gamma 0.57721566490153286061) (log-x (log x)))
(let loop ((k 1) (fact 1) (num log-x)
(sum (+ gamma (log log-x) log-x)))
(if (< 100 k) sum
(let* ((k (+ k 1))
(fact (* fact k))
(num (* num log-x))
(sum (+ sum (/ num fact k))))
(loop k fact num sum))))))
(define (factors n) ; trial division
(let loop ((n n) (fs (list)))
(if (even? n) (loop (/ n 2) (cons 2 fs))
(if (= n 1) (if (null? fs) (list 1) fs)
(let loop ((n n) (f 3) (fs fs))
(cond ((< n (* f f)) (reverse (cons n fs)))
((zero? (modulo n f))
(loop (/ n f) f (cons f fs)))
(else (loop n (+ f 2) fs))))))))
(define (mobius-mu n)
(if (= n 1) 1
(let loop ((fs (factors n)) (prev 0) (m 1))
(cond ((null? fs) m)
((= (car fs) prev) 0)
(else (loop (cdr fs) (car fs) (- m)))))))
(define riemann-r
(let ((ms (let loop ((n 1) (k 100) (ms (list)))
(if (zero? k) (reverse ms)
(let ((m (mobius-mu n)))
(if (zero? m) (loop (+ n 1) k ms)
(loop (+ n 1) (- k 1) (cons (* m n) ms))))))))
(lambda (x)
(let loop ((ms ms) (sum 0))
(if (null? ms) sum
(let* ((m (car ms)) (m-abs (abs m)) (m-recip (/ m)))
(loop (cdr ms) (+ sum (* m-recip (logint (expt x (/ m-abs))))))))))))
(display "10^3") (display #\tab)
(display (- (inexact->exact (round (logint 1e3))) 168)) (display #\tab)
(display (- (inexact->exact (round (riemann-r 1e3))) 168)) (newline)
(display "10^6") (display #\tab)
(display (- (inexact->exact (round (logint 1e6))) 78498)) (display #\tab)
(display (- (inexact->exact (round (riemann-r 1e6))) 78498)) (newline)
(display "10^9") (display #\tab)
(display (- (inexact->exact (round (logint 1e9))) 50847534)) (display #\tab)
(display (- (inexact->exact (round (riemann-r 1e9))) 50847534)) (newline)
(display "10^12") (display #\tab)
(display (- (inexact->exact (round (logint 1e12))) 37607912018)) (display #\tab)
(display (- (inexact->exact (round (riemann-r 1e12))) 37607912018)) (newline)