; fibonacci primes

(define (greplin2 n)
  (let loop ((f-2 0) (f-1 1) (f 1))
    (if (and (< n f) (prime? f))
        (sum (factors (+ f 1)))
        (loop f-1 f (+ f-1 f)))))

(define (factors n)
  (let loop ((n n) (x 2) (fs '()))
    (cond ((< n (* x x)) (reverse (cons n fs)))
          ((zero? (modulo n x)) (loop (/ n x) x (cons x fs)))
          (else (loop n (+ x 1) fs)))))

(define (witness? a n)
  (let loop ((r 0) (s (- n 1)))
    (if (even? s) (loop (+ r 1) (/ s 2))
      (if (= (expm a s n) 1) #t
        (let loop ((j 0) (s s))
          (cond ((= j r) #f)
                ((= (expm a s n) (- n 1)) #t)
                (else (loop (+ j 1) (* s 2)))))))))

(define (expm b e m)
  (define (m* x y) (modulo (* x y) m))
  (cond ((zero? e) 1)
        ((even? e) (expm (m* b b) (/ e 2) m))
        (else (m* b (expm (m* b b) (/ (- e 1) 2) m)))))

(define (prime? n)
  (cond ((< n 2) #f) ((= n 2) #t) ((even? n) #f)
        (else (let loop ((k 50))
                (cond ((zero? k) #t)
                      ((not (witness? (randint 1 n) n)) #f)
                      (else (loop (- k 1))))))))

(define rand
  (let ((a 16807) (m 2147483647) (seed 1043618065))
    (lambda args
      (if (pair? args) (set! seed (car args)))
      (set! seed (modulo (* a seed) m))
      (/ seed m))))

(define (randint first past)
  (+ first (inexact->exact (floor (* (rand) (- past first))))))

(define (sum xs) (apply + xs))

(display (greplin2 227000))

352