codepad
[
create a new paste
]
login
|
about
Language:
C
C++
D
Haskell
Lua
OCaml
PHP
Perl
Plain Text
Python
Ruby
Scheme
Tcl
; 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)
Private
[
?
]
Run code