[ create a new paste ] login | about

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

programmingpraxis - Scheme, pasted on Sep 9:
; tribonacci numbers

(define (iterate n f . bs)
  (let loop ((n n) (b (car bs)) (bs (cdr bs)) (xs '()))
    (if (zero? n) (reverse xs)
      (let ((new-bs (append bs (list (apply f b bs)))))
        (loop (- n 1) (car new-bs) (cdr new-bs) (cons b xs))))))

(define (tribs n)
  (iterate n + 1 1 2))

(display (tribs 25)) (newline)

(define (make-matrix rows columns . value)
  (do ((m (make-vector rows)) (i 0 (+ i 1)))
      ((= i rows) m)
    (if (null? value)
        (vector-set! m i (make-vector columns))
        (vector-set! m i (make-vector columns (car value))))))

(define (matrix-rows x) (vector-length x))

(define (matrix-cols x) (vector-length (vector-ref x 0)))

(define (matrix-ref m i j) (vector-ref (vector-ref m i) j))

(define (matrix-set! m i j x) (vector-set! (vector-ref m i) j x))

(define-syntax for
  (syntax-rules ()
    ((for (var first past step) body ...)
      (let ((ge? (if (< first past) >= <=)))
        (do ((var first (+ var step)))
            ((ge? var past))
          body ...)))
    ((for (var first past) body ...)
      (let* ((f first) (p past) (s (if (< first past) 1 -1)))
        (for (var f p s) body ...)))
    ((for (var past) body ...)
      (let* ((p past)) (for (var 0 p) body ...)))))

(define (matrix-multiply a b)
  (let ((ar (matrix-rows a)) (ac (matrix-cols a))
        (br (matrix-rows b)) (bc (matrix-cols b)))
    (if (not (= ac br))
        (error 'matrix-multiply "incompatible matrices")
        (let ((c (make-matrix ar bc 0)))
          (for (i ar)
            (for (j bc)
              (for (k ac)
                (matrix-set! c i j
                  (+ (matrix-ref c i j)
                     (* (matrix-ref a i k)
                        (matrix-ref b k j)))))))
          c))))

(define (matrix-power m n)
  (cond ((= n 1) m)
        ((even? n)
          (matrix-power
            (matrix-multiply m m)
            (/ n 2)))
        (else (matrix-multiply m
                (matrix-power
                  (matrix-multiply m m)
                  (/ (- n 1) 2))))))

(define (trib n)
  (matrix-ref
    (matrix-power
      #( #(1 1 0)
         #(1 0 1)
         #(1 0 0))
      n) 2 0))

(display (trib 1000)) (newline)

(display (exact->inexact (/ (trib 1000) (trib 999))))


Output:
1
2
3
(1 1 2 4 7 13 24 44 81 149 274 504 927 1705 3136 5768 10609 19513 35890 66012 121415 223317 410744 755476 1389537)
1499952522327196729941271196334368245775697491582778125787566254148069690528296568742385996324542810615783529390195412125034236407070760756549390960727215226685972723347839892057807887049540341540394345570010550821354375819311674972209464069786275283520364029575324
1.839286755214161


Create a new paste based on this one


Comments: