; pythagorean quadruples
 
(define-syntax fold-of
  (syntax-rules (range in is)
    ((_ "z" f b e) (set! b (f b e)))
    ((_ "z" f b e (v range fst pst stp) c ...)
      (let* ((x fst) (p pst) (s stp)
             (le? (if (positive? s) <= >=)))
        (do ((v x (+ v s))) ((le? p v) b)
          (fold-of "z" f b e c ...))))
    ((_ "z" f b e (v range fst pst) c ...)
      (let* ((x fst) (p pst) (s (if (< x p) 1 -1)))
        (fold-of "z" f b e (v range x p s) c ...)))
    ((_ "z" f b e (v range pst) c ...)
      (fold-of "z" f b e (v range 0 pst) c ...))
    ((_ "z" f b e (x in xs) c ...)
      (do ((t xs (cdr t))) ((null? t) b)
        (let ((x (car t)))
          (fold-of "z" f b e c ...))))
    ((_ "z" f b e (x is y) c ...)
      (let ((x y)) (fold-of "z" f b e c ...)))
    ((_ "z" f b e p? c ...)
      (if p? (fold-of "z" f b e c ...)))
    ((_ f i e c ...)
      (let ((b i)) (fold-of "z" f b e c ...)))))
 
(define-syntax list-of (syntax-rules ()
  ((_ arg ...) (reverse (fold-of
    (lambda (d a) (cons a d)) '() arg ...)))))
 
(define (isqrt n)
  (if (not (and (positive? n) (integer? n)))
      (error 'isqrt "must be positive integer")
      (let loop ((x n))
        (let ((y (quotient (+ x (quotient n x)) 2)))
          (if (< y x) (loop y) x)))))
 
(define (pyquad n)
  (list-of (list a b c d)
    (a range 1 (+ n 1))
    (b range a (+ n 1))
    (c range b (+ n 1))
    (s is (+ (* a a) (* b b) (* c c)))
    (d is (isqrt s))
    (= (* d d) s)))
 
(display (pyquad 25))

OyBweXRoYWdvcmVhbiBxdWFkcnVwbGVzCgooZGVmaW5lLXN5bnRheCBmb2xkLW9mCiAgKHN5bnRheC1ydWxlcyAocmFuZ2UgaW4gaXMpCiAgICAoKF8gInoiIGYgYiBlKSAoc2V0ISBiIChmIGIgZSkpKQogICAgKChfICJ6IiBmIGIgZSAodiByYW5nZSBmc3QgcHN0IHN0cCkgYyAuLi4pCiAgICAgIChsZXQqICgoeCBmc3QpIChwIHBzdCkgKHMgc3RwKQogICAgICAgICAgICAgKGxlPyAoaWYgKHBvc2l0aXZlPyBzKSA8PSA+PSkpKQogICAgICAgIChkbyAoKHYgeCAoKyB2IHMpKSkgKChsZT8gcCB2KSBiKQogICAgICAgICAgKGZvbGQtb2YgInoiIGYgYiBlIGMgLi4uKSkpKQogICAgKChfICJ6IiBmIGIgZSAodiByYW5nZSBmc3QgcHN0KSBjIC4uLikKICAgICAgKGxldCogKCh4IGZzdCkgKHAgcHN0KSAocyAoaWYgKDwgeCBwKSAxIC0xKSkpCiAgICAgICAgKGZvbGQtb2YgInoiIGYgYiBlICh2IHJhbmdlIHggcCBzKSBjIC4uLikpKQogICAgKChfICJ6IiBmIGIgZSAodiByYW5nZSBwc3QpIGMgLi4uKQogICAgICAoZm9sZC1vZiAieiIgZiBiIGUgKHYgcmFuZ2UgMCBwc3QpIGMgLi4uKSkKICAgICgoXyAieiIgZiBiIGUgKHggaW4geHMpIGMgLi4uKQogICAgICAoZG8gKCh0IHhzIChjZHIgdCkpKSAoKG51bGw/IHQpIGIpCiAgICAgICAgKGxldCAoKHggKGNhciB0KSkpCiAgICAgICAgICAoZm9sZC1vZiAieiIgZiBiIGUgYyAuLi4pKSkpCiAgICAoKF8gInoiIGYgYiBlICh4IGlzIHkpIGMgLi4uKQogICAgICAobGV0ICgoeCB5KSkgKGZvbGQtb2YgInoiIGYgYiBlIGMgLi4uKSkpCiAgICAoKF8gInoiIGYgYiBlIHA/IGMgLi4uKQogICAgICAoaWYgcD8gKGZvbGQtb2YgInoiIGYgYiBlIGMgLi4uKSkpCiAgICAoKF8gZiBpIGUgYyAuLi4pCiAgICAgIChsZXQgKChiIGkpKSAoZm9sZC1vZiAieiIgZiBiIGUgYyAuLi4pKSkpKQoKKGRlZmluZS1zeW50YXggbGlzdC1vZiAoc3ludGF4LXJ1bGVzICgpCiAgKChfIGFyZyAuLi4pIChyZXZlcnNlIChmb2xkLW9mCiAgICAobGFtYmRhIChkIGEpIChjb25zIGEgZCkpICcoKSBhcmcgLi4uKSkpKSkKICAgCihkZWZpbmUgKGlzcXJ0IG4pCiAgKGlmIChub3QgKGFuZCAocG9zaXRpdmU/IG4pIChpbnRlZ2VyPyBuKSkpCiAgICAgIChlcnJvciAnaXNxcnQgIm11c3QgYmUgcG9zaXRpdmUgaW50ZWdlciIpCiAgICAgIChsZXQgbG9vcCAoKHggbikpCiAgICAgICAgKGxldCAoKHkgKHF1b3RpZW50ICgrIHggKHF1b3RpZW50IG4geCkpIDIpKSkKICAgICAgICAgIChpZiAoPCB5IHgpIChsb29wIHkpIHgpKSkpKQoKKGRlZmluZSAocHlxdWFkIG4pCiAgKGxpc3Qtb2YgKGxpc3QgYSBiIGMgZCkKICAgIChhIHJhbmdlIDEgKCsgbiAxKSkKICAgIChiIHJhbmdlIGEgKCsgbiAxKSkKICAgIChjIHJhbmdlIGIgKCsgbiAxKSkKICAgIChzIGlzICgrICgqIGEgYSkgKCogYiBiKSAoKiBjIGMpKSkKICAgIChkIGlzIChpc3FydCBzKSkKICAgICg9ICgqIGQgZCkgcykpKQoKKGRpc3BsYXkgKHB5cXVhZCAyNSkp