; 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)))

((10 2001) (7 2004) (5 2007) (12 2009) (8 2012) (7 2015) (1 2018) (3 2018) (10 2020) (8 2023) (5 2026) (12 2028) (9 2031) (6 2034) (1 2037) (3 2037) (10 2039) (7 2042) (4 2045) (12 2047) (9 2050) (6 2053) (1 2056) (3 2056) (10 2058) (7 2061) (4 2064) (12 2066) (8 2069) (5 2072) (1 2075) (3 2075) (10 2077) (7 2080) (5 2083) (12 2085) (8 2088) (6 2091) (1 2094) (3 2094) (10 2096) (8 2099))