; once in a blue moon
(define (julian year month day)
(let* ((a (quotient (- 14 month) 12))
(y (+ year 4800 (- a)))
(m (+ month (* 12 a) -3)))
(+ day
(quotient (+ (* 153 m) 2) 5)
(* 365 y)
(quotient y 4)
(- (quotient y 100))
(quotient y 400)
(- 32045))))
(define (gregorian julian)
(let* ((j (+ julian 32044))
(g (quotient j 146097))
(dg (modulo j 146097))
(c (quotient (* (+ (quotient dg 36524) 1) 3) 4))
(dc (- dg (* c 36524)))
(b (quotient dc 1461))
(db (modulo dc 1461))
(a (quotient (* (+ (quotient db 365) 1) 3) 4))
(da (- db (* a 365)))
(y (+ (* g 400) (* c 100) (* b 4) a))
(m (- (quotient (+ (* da 5) 308) 153) 2))
(d (+ da (- (quotient (* (+ m 4) 153) 5)) 122))
(year (+ y (- 4800) (quotient (+ m 2) 12)))
(month (+ (modulo (+ m 2) 12) 1))
(day (+ d 1)))
(values year month day)))
(define (blue-moons from to)
(define (int x) (inexact->exact (floor x)))
(let* ((delta 29.530588853) (offset 23.9704840005)
(moon (+ offset (* delta (floor (/ from delta))))))
(let loop ((moon moon) (prev 0) (moons (list)))
(if (< to moon) (reverse moons)
(let-values (((y m d) (gregorian (int moon))))
(loop (+ moon delta) m
(if (= m prev)
(cons (list m y) moons)
moons)))))))
(display (blue-moons (julian 2001 1 1) (julian 2100 12 31)))