; two integrals
(define (exp-integral x)
(let* ((gamma 0.5772156649015328606065))
(let loop ((k 1) (fact 1) (z x) (sum (+ gamma (log x))))
(let ((term (/ z k fact)))
(if (< term 1e-17) sum
(loop (+ k 1) (* fact (+ k 1)) (* z x) (+ sum term)))))))
(define (log-integral x)
(let* ((gamma 0.5772156649015328606065)
(logx (log x)) (sum (+ gamma (log logx))))
(let loop ((k 1) (fact 1) (logx^k logx) (sum sum))
(let* ((term (/ logx^k fact k)) (sum (+ sum term))
(newk (+ k 1)) (newfact (* fact newk)))
(if (< term 1e-17) sum
(loop newk newfact (* logx^k logx) sum))))))
(define (offset-log-integral x)
(- (log-integral x) 1.04516378011749278))
(display (exp-integral (log 1000000))) (newline)
(display (log-integral 1000000)) (newline)
(display (inexact->exact (round (offset-log-integral 1e21)))) (newline)