[ create a new paste ] login | about

Project: programmingpraxis
Link: http://programmingpraxis.codepad.org/Fvwfmjc5    [ raw code | output | fork ]

programmingpraxis - Scheme, pasted on Jul 28:
; 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)


Output:
1
2
3
4
10^3	10	0
10^6	130	29
10^9	1701	-79
10^12	38263	-1476


Create a new paste based on this one


Comments: