;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;                                                                      ;;
;;  http://www.megabaud.fi/~karttu/matikka/Schemuli/intfuns1.scm        ;;
;;  - Often needed integer functions.                                   ;;
;;                                                                      ;;
;;  (Schemuli stands for the "Useful Library for Scheme" and "matikka"  ;;
;;   is a certain kind of fish, believe it or not).                     ;;
;;                                                                      ;;
;;  Coded by Antti Karttunen (my_firstname.my_surname@iki.fi), 2002     ;;
;;                                                                      ;;
;;  This Scheme-code is in Public Domain and runs (at least)            ;;
;;  in MIT Scheme Release 7.6.0, for which one can find documentation   ;;
;;  and the pre-compiled binaries (for various OS's running in          ;;
;;  Intel x86 architecture) under the URL:                              ;;
;;  http://www.swiss.ai.mit.edu/projects/scheme/                        ;;
;;                                                                      ;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

;;
;; Compile as:
;;
;; (cf "c:\\matikka\\Schemuli\\definech" "c:\\matikka\\Schemuli\\")
;; 
;; (load  "c:\\matikka\\Schemuli\\definech")
;; ;Loading "matikka\\schemuli\\definech.com" -- done
;; ;Value: definec
;; 
;; (fluid-let ((sf/default-syntax-table user-initial-environment))
;;   (cf "c:\\matikka\\Schemuli\\intfuns1" "c:\\matikka\\Schemuli\\")
;; )
;; 

;; From http://www.swiss.ai.mit.edu/projects/scheme/documentation/user_5.html
;;
;; If you redefine some global name in your code, for example, car, cdr
;; and cons, you should indicate it in the declaration as:
;; (declare (usual-integrations car cdr cons))

;; (Beware of using an argument name like list in the function definitions?)

(declare (usual-integrations))

(definec (fibo nakki) ;; I.e., A000045
   (if (< nakki 2)
       nakki
       (+ (fibo (-1+ nakki)) (fibo (- nakki 2)))
   )
)


;; Factorial with syntax modified from the post-fix to pre-fix !:
(definec (! n) (if (zero? n) 1 (* n (! (-1+ n))))) ;; A000142

(define (!cut i j) ;; Compute the product i*(i+1)*(i+2)*...*j
   (cond ((eq? i j) i)
         ((> i j) 1)
         (else (* i (!cut (1+ i) j)))
   )
)

;; (C r n) stands for the binomial {r choose n}
(define (C r n) (/ (!cut (1+ n) r) (! (- r n))))




(definec (binomial_n_2 n) (/ (* (-1+ n) n) 2))


(define (factbaseR->n rex) ;; Convert the reversed factorial expansion list to an integer.
   (let loop ((rex rex) (n 0) (i 1))
      (cond ((null? rex) n)
            (else (loop (cdr rex) (+ n (* (! i) (car rex))) (1+ i)))
      )
   )
)


(define (n->factbase n) ;; Convert an integer to a factorial expansion list.
   (let loop ((n n) (fex (if (zero? n) (list 0) (list))) (i 2))
      (cond ((zero? n) fex)
            (else (loop (floor->exact (/ n i)) (cons (modulo n i) fex) (1+ i)))
      )
   )
)

(define (baselist->n base bex) ;; Convert base n expansion list to an integer.
   (let loop ((bex bex) (n 0))
      (cond ((null? bex) n)
            (else (loop (cdr bex) (+ (* n base) (car bex))))
      )
   )
)

(define (baselist-as-binary  lista) (baselist->n  2 lista))
(define (baselist-as-ternary lista) (baselist->n  3 lista))
(define (baselist-as-decimal lista) (baselist->n 10 lista))

(define (halve n) (/ n 2))

(define (lsb n) (modulo n 2))

(define (floor-log-2 n) (floor->exact (/ (log n) (log 2))))

;; Note for n > 0, (floor-log-2 n) = (- (binwidth n) 1)

(define (binwidth n) ;; A029837 (with a(0)=0 instead of 1)
  (let loop ((n n) (i 0)) ;;  or A036377 (with offset=0 instead of 1)
     (if (zero? n)
         i
         (loop (floor->exact (/ n 2)) (1+ i))
     )
  )
)

(definec (binrev nn) ;; A030101
  (let loop ((z 0) (n nn))
    (if (zero? n)
        z
        (loop (+ (* 2 z) (modulo n 2))
              (fix:lsh n -1) ;; n >>= 1
        )
    )
  )
)

(define (A054429 n) ;; -> (0),1,3,2,7,6,5,4,15,14,13,12,11,10,9,8,31,30,...
  (if (zero? n)
      n
      (-1+ (- (* 3 (expt 2 (floor-log-2 n))) n))
  )
)


;; ;; Moser-de Bruijn sequence: 0,1,4,5,16,17,20,21,64,65,68,69,80,...
;; (definec (A000695 n) ;; Expand bits, 0->00, 1->01, i.e. from base-2 to base-4.
;;   (if (zero? n)
;;       n
;;       (+ (modulo n 2) (fix:lsh (A000695 (fix:lsh n -1)) 2))
;;   )
;; )
;; 
;; (definec (A059905 n) ;; Take the even-positioned bits and contract them.
;;   (if (zero? n)
;;       n
;;       (+ (modulo n 2) (fix:lsh (A059905 (fix:lsh n -2)) 1))
;;   )
;; )
;; 

;; Moser-de Bruijn sequence: 0,1,4,5,16,17,20,21,64,65,68,69,80,...
(define (A000695 n) ;; Expand bits, 0->00, 1->01, i.e. from base-2 to base-4.
  (if (zero? n)
      n
      (+ (modulo n 2) (* 4 (A000695 (floor->exact (/ n 2)))))
  )
)

(define (A059905 n) ;; Take the even-positioned bits and contract them.
  (if (zero? n)
      n
      (+ (modulo n 2) (* 2 (A059905 (floor->exact (/ n 4)))))
  )
)


(definec (A000695c n) (A000695 n)) ;; The cached variant. Be careful with this!
(definec (A059905c n) (A055905 n)) ;; The cached variant. Be careful with this!

(define (A059906 n) (A059905 (floor->exact (/ n 2)))) ;; Take the odd bits and contract
(define (A059906*2 n) (* 2 (A059906 n)))



(define (A006519 n) ;; Highest power of 2 dividing n: 1,2,1,4,1,2,1,8,1,2,1,4,1,2,1,16
   (cond ((zero? n) 0)
         (else (let loop ((n n) (i 1))
                 (cond ((odd? n) i)
                       (else (loop (fix:lsh n -1) (fix:lsh i 1)))
                 )
               )
         )
   )
)


(define (A007814 n) ;; Exponent of the A006519.
   (cond ((zero? n) 0)
         (else (let loop ((n n) (i 0))
                 (cond ((odd? n) i)
                       (else (loop (fix:lsh n -1) (1+ i)))
                 )
               )
         )
   )
)


(definec (A036987 n) (if (eq? (1+ n) (A006519 (1+ n))) 1 0))

(define (A073265bi n k)
   (cond ((zero? n) n)
         ((zero? k) k)
         ((> k n) 0)
         ((eq? 1 k) (if (eq? n (A006519 n)) 1 0))
         (else
           (let sumloop ((i (floor-log-2 (-1+ n))) (s 0))
                 (cond ((negative? i) s)
                       (else
                         (sumloop (-1+ i)
                                  (+ s (A073265bi (- n (fix:lsh 1 i)) (-1+ k)))
                         )
                       )
                 )
           )
         )
   )
)


;; A073265 - A073270 reserved for us.

(define (A073265 n) (A073265bi (1+ (A025581 n)) (1+ (A002262 n))))
(define (A073266 n) (A073265bi (1+ (A003056 n)) (1+ (A002262 n))))
(define (A073267 n) (A073265bi n 2)) ;; Occurs also as the FIX-counts of form 105. 6o6
;; (A073265bi n 1) ;; gives the characteristic function of A000079,
;; i.e. A036987 with offset 1 instead of 0, and a(0) = 0.

;; A073265(6,1) = 0
;; A073265(6,2) = 2 4+2, 2+4 -> 3+(1)+1, 1+(1)+3
;; A073265(6,3) = 4 2+2+2, 4+1+1, 1+4+1, 1+1+4 -> 1+(1)+1+(1)+1, 3+(1)+0+(1)+0, etc.
;; A073265(6,4) = 6 2+2+1+1, 2+1+2+1, 2+1+1+2, 1+2+1+2, 1+2+2+1, 1+1+2+2
;; A073265(6,5) = 5 2+1+1+1+1 * 5 --> 1+(1)+0+(1)+0+(1)+0+(1)+0
;; A073265(6,6) = 1 1+1+1+1+1+1 --> 0+(1)+0+(1)+0+(1)+0+(1)+0+(1)+0

(define (A073268 n) ;; The fix sequence for form 41, SwapBinTree o SwapDownCar (0 o 6)
  (if (zero? n)
      1
      (let sumloop ((i (floor->exact (/ (log n) (log 2)))) (s 0))
        (cond ((negative? i) s)
              (else (sumloop (-1+ i) (+ s (A000108 (- n (expt 2 i))))))
        )
      )
  )
)


(define (A073345bi n k)
  (cond ((zero? n) (if (zero? k) 1 0))
        ((zero? k) k)
        ((> k n) 0)
        (else
          (let ((half-n (fix:lsh (-1+ n) -1))
                (k-1 (-1+ k))
               )
           (+ (* 2 (add (lambda (i)
                          (* (A073345bi (- (-1+ n) i) k-1)
                             (add (lambda (j) (A073345bi i j))
                                  0 k-1
                             )
                          )
                        )
                        0 half-n
                   )
              )
              (* 2 (add (lambda (i)
                          (* (A073345bi (- (-1+ n) i) k-1)
                             (add (lambda (j) (A073345bi i j))
                                  0 (- k 2)
                             )
                          )
                        )
                        (1+ half-n) (-1+ n)
                   )
              )
              (if (odd? n) (* -1 (expt (A073345bi half-n k-1) 2)) 0)
           ) ;; +
          ) ;; let
        ) ;; else
  ) ;; cond
)


(define (A073346bi n k)
  (cond ((zero? k) (A036987 n))
        ((zero? n) 0)
        ((> k n) 0)
        (else
          (let ((half-n (fix:lsh (-1+ n) -1))
                (k-1 (-1+ k))
               )
           (+ (* 2 (add (lambda (i)
                          (* (A073346bi (- (-1+ n) i) k-1)
                             (add (lambda (j) (A073346bi i j))
                                  0 k-1
                             )
                          )
                        )
                        0 half-n
                   )
              )
              (* 2 (add (lambda (i)
                          (* (A073346bi (- (-1+ n) i) k-1)
                             (add (lambda (j) (A073346bi i j))
                                  0 (- k 2)
                             )
                          )
                        )
                        (1+ half-n) (-1+ n)
                   )
              )
              (if (odd? n) (* -1 (expt (A073346bi half-n k-1) 2)) 0)
              (if (= 1 k) (* -1 (A036987 n)) 0)
           ) ;; +
          ) ;; let
        ) ;; else
  ) ;; cond
)


(define (A073345 n) (A073345bi (A025581 n) (A002262 n)))
(define (A073346 n) (A073346bi (A025581 n) (A002262 n)))
(define (A073429 n) (A073345bi (A003056 n) (A002262 n)))
(define (A073430 n) (A073346bi (A003056 n) (A002262 n)))

(define (A073431 n)
   (cond
     ((zero? n) 1)
     (else
        (/ (add (lambda (i) (add (lambda (j) (A073346bi n j)) 0 (A007814 i)))
                1 (expt 2 (-1+ n))
           )
           (expt 2 (-1+ n))
        )
     )
   )
)


;; For some reason this seems to produce the same answers:
;; (Some peculiar interaction of A007814 & A073346 ?)

(define (A073431v2 n)
   (cond
     ((zero? n) 1)
     (else
       (-1+ (/ (add (lambda (i) (A073346bi n (A007814 i)))
                    1 (expt 2 (-1+ n))
               )
               (expt 2 (- n 2))
            )
       )
     )
   )
)



(definec (A048678 n) ;; Rewrite to Zeckendorf-expansion: 0->0, 1->01
  (cond ((zero? n) n)
        ((even? n) (* 2 (A048678 (floor->exact (/ n 2)))))
        (else      (1+ (* 4 (A048678 (floor->exact (/ n 2))))))
  )
)


(define (A022290 n) ;; I.e. interpret_as_zeckendorf_expansion
  (let loop ((n n) (s 0) (i 2))
     (cond ((zero? n) s)
           (else (loop (floor->exact (/ n 2))
                       (+ s (* (modulo n 2) (fibo i)))
                       (1+ i)
                 )
           )
     )
  )
)


(define *Sqrt5* (sqrt 5))
(define *Phi* (/ (1+ *Sqrt5*) 2))
(define *LogPhi* (log *Phi*))

(define (A072648 n) ;; An approximate "inverse" of A000045 (of the fibonacci numbers)
   (cond ((zero? n) n)
         (else (floor->exact (/ (log (* n *Sqrt5*)) *LogPhi*)))
   )
)

;; 1,2,3,4,5,6,7,8,9,A,B,C,D,E,F,         ,21
;; 1,2,3,3,4,4,4,5,5,5,5,5,6,6,6,6,6,6,6,6,7,...
(define (A072649 n) ;; n occurs fibo(n) times.
  (let ((b (A072648 n)))
     (+ -1 b (floor->exact (/ n (fibo (1+ b)))))
  )
)


(define (A072650 n) ;; rewrite_0to0_x1to1
  (let loop ((n n) (s 0) (i 0))
     (cond ((zero? n) s)
           ((even? n) (loop (floor->exact (/ n 2)) s (1+ i)))
           (else (loop (floor->exact (/ n 4)) (+ s (expt 2 i)) (1+ i)))
     )
  )
)



(definec (A003714c n) ;; The cached variant.
  (cond ((< n 3) n)
        (else (+ (expt 2 (-1+ (A072649 n)))
                 (A003714c (- n (fibo (1+ (A072649 n)))))
              )
        )
  )
)


(define (A003714 n) ;; Iterative (tail-recursive) variant, not cached.
   (let loop ((n n) (s 0))
     (cond ((< n 3) (+ s n))
           (else (loop (- n (fibo (1+ (A072649 n))))
                       (+ s (expt 2 (-1+ (A072649 n))))
                 )
           )
     )
  )
)


(define (A048679 n) (A072650 (A003714 n)))
(define (A048680 n) (A022290 (A048678 n)))


;; (map A025581 (cons 0 (iota 20))) --> (0 1 0 2 1 0 3 2 1 0 4 3 2 1 0 5 4 3 2 1 0)

;; (definec (A025581 n) ;; The X component (column) of square {0..inf} arrays
;;   (- (binomial_n_2 (1+ (floor->exact (+ (/ 1 2) (sqrt (* 2 (1+ n))))))) (1+ n))
;; )
;; 
;; ;; (map A002262 (cons 0 (iota 20))) --> (0 0 1 0 1 2 0 1 2 3 0 1 2 3 4 0 1 2 3 4 5)
;; (definec (A002262 n) ;; The Y component (row) of square {0..inf} arrays
;;   (- n (binomial_n_2 (floor->exact (+ (/ 1 2) (sqrt (* 2 (1+ n)))))))
;; )

;; At some point these will produce incorrect values, because of the
;; limited precision of IEEE 64-bit floating point numbers.
;; What is that point, and how to recode these with strictly fixnum-only
;; way? (I need a fixnum-only square root algorithm...)
(definec (A025581 n) ;; The X component (column) of square {0..inf} arrays
  (- (binomial_n_2 (1+ (floor->exact (flo:+ 0.5 (flo:sqrt (exact->inexact (* 2 (1+ n)))))))) (1+ n))
)

;; (map A002262 (cons 0 (iota 20))) --> (0 0 1 0 1 2 0 1 2 3 0 1 2 3 4 0 1 2 3 4 5)
(definec (A002262 n) ;; The Y component (row) of square {0..inf} arrays
  (- n (binomial_n_2 (floor->exact (flo:+ 0.5 (flo:sqrt (exact->inexact (* 2 (1+ n))))))))
)


(define (A002024 n) ;; repeat n n times, starting from n = 1.
  (floor->exact (+ (/ 1 2) (sqrt (* 2 n))))
)

(define (A003056 n) ;; repeat n n+1 times, starting from n = 0.
  (floor->exact (- (sqrt (* 2 (1+ n))) (/ 1 2)))
)

(define (A001477 n) n)

(define (packA001477 x y) (/ (+ (expt (+ x y) 2) x (* 3 y)) 2))
(define (packA061579 x y) (/ (+ (expt (+ x y) 2) (* 3 x) y) 2))

;; I.e. (define (id n) (packA001477 (A025581 n) (A002262 n)))
;; and

(define (A061579 n) (packA061579 (A025581 n) (A002262 n)))

(define (A038722 n) (if (zero? n) n (1+ (A061579 (-1+ n)))))

;; Gives id (A001477): (+ (A000695 (A059905 n)) (* 2 (A000695 (A059906 n))))
(define (A057300 n) (+ (A000695 (A059906 n)) (* 2 (A000695 (A059905 n)))))
(define (A054238 n) (+ (A000695 (A025581 n)) (* 2 (A000695 (A002262 n)))))
(define (A054239 n) (packA001477 (A059905 n) (A059906 n)))

(define (A072661 n) (A059905 (A048679 n)))
(define (A072662 n) (A059906 (A048679 n)))


(define (packA054238 x y) (+ (A000695 x) (* 2 (A000695 y))))
(define (packA054238tr x y) (+ (A000695 y) (* 2 (A000695 x))))

(define (packA048680oA054238 x y) (A048680 (packA054238 x y)))
(define (packA048680oA054238tr x y) (A048680 (packA054238tr x y)))

(define (A072793 n) (A048680 (A054238 n))) ;; EIS # reserved for packA048680oA054238
(define (A072794 n) (A054239 (A048679 n))) ;; inverse for above.

;; A972732 - A072741 reserved for us.

(define (A072732 n) (packA072732 (A025581 n) (A002262 n)))

(define (packA072732 x y)
   (let ((x-y (- x y)))
     (cond ((<= x-y 0) ;; i.e. (x <= y)
               (packA001477
                  (+ (* 2 x) (modulo x-y 2))
                  (+ (* 2 x) (floor->exact (/ (1+ (- x-y)) 2)))
               )
           )
;;         ((< x-y 2) ;; either x=y, or x=y+1.
;;             (packA001477 (+ (* 2 y) x-y) (* 2 y))
;;         )
           (else  ;; x > y. (floor->exact (/ -1 2)) has to be -1 !!!
               (packA001477
                  (+ (* 2 (1+ y)) (floor->exact (/ (- x-y 2) 2)))
                  (+ (* 2 y) (modulo (1+ x-y) 2))
               )
           )
     )
   )
)

(define (A072733 n) (packA072733 (A025581 n) (A002262 n)))

(define (packA072733 x y)
     (cond ((<= x y) ;; i.e. (x <= y)
             (let ((half-x (floor->exact (/ x 2))))
               (packA001477
                  half-x
                  (+ half-x
                     (* 2 (- y (* 2 half-x) (modulo x 2)))
                     (modulo x 2)
                  )
               )
             )
           )
           (else  ;; x > y
             (let ((half-y (floor->exact (/ y 2))))
               (packA001477
                  (+ 1 half-y
                     (* 2 (- (-1+ x) (* 2 half-y) (modulo y 2)))
                     (modulo y 2)
                  )
                  half-y
               )
             )
           )
     )
)


(define (A072734 n) (packA072734 (A025581 n) (A002262 n)))

;; And code this as an inverse of A072735 (the next one)...

;; Using this as a NxN->N packing bijection for a global arithmetic
;; unranking function for Catalan structures gives a better result
;; than A072642:
;; (define weird1tr (max-n-fun-with-arithrank-scheme packA072734))
;; (map binwidth (map weird1tr (iota0 10)))
;; (0 1 2 4 6  8 13 23 43 82 160)
;; And using it's "4th power" gives even better results in a certain range:
;; (0 1 2 8 9 11 15 23 39 71 134)
;;

(define (packA072734 x y)
   (let ((x-y (- x y)))
     (cond
           ((negative? x-y) ;; i.e. (y > x)
               (packA001477
                  (+ (* 2 x) (modulo (1+ x-y) 2))
                  (+ (* 2 x) (floor->exact (/ (+ (- x-y) (modulo x-y 2)) 2)))
               )
           )
           ((< x-y 3) ;; i.e. (x >= y and x <= y+2)
               (packA001477 (+ (* 2 y) x-y) (* 2 y))
           )
           (else  ;; x >= y+3
               (packA001477
                  (+ (* 2 y) (floor->exact (/ (1+ x-y) 2)) (modulo (1+ x-y) 2))
                  (+ (* 2 y) (modulo x-y 2))
               )
           )
     )
   )
)


(define (A072735 n) (packA072735 (A025581 n) (A002262 n)))

(define (packA072735 x y)
     (cond ((<= x y) ;; i.e. (x <= y)
             (let ((half-x (floor->exact (/ x 2))))
               (packA001477
                  half-x
                  (+ half-x
                     (* 2 (- y (* 2 half-x)))
                     (modulo x 2)
                     (if (and (eq? x y) (even? x)) 0 -1)
                  )
               )
             )
           )
           (else  ;; x > y
             (let ((half-y (floor->exact (/ y 2))))
               (packA001477
                  (+ half-y
                     (* 2 (- (-1+ x) (* 2 half-y)))
                     (modulo y 2)
                     (if (and (eq? x (1+ y)) (even? y)) 1 0)
                  )
                  half-y
               )
             )
           )
     )
)


(define (A072736 n) (A025581 (A072733 n))) ;; X-projection of A072732
(define (A072737 n) (A002262 (A072733 n))) ;; Y-projection of A072732

(define (A072738 n) (A025581 (A072732 n))) ;; X-projection of A072733
(define (A072739 n) (A002262 (A072732 n))) ;; Y-projection of A072733

(define (A072740 n) (A025581 (A072735 n))) ;; X-projection of A072734
(define (A072741 n) (A002262 (A072735 n))) ;; Y-projection of A072734

;; Also A072781 - A072800 reserved for us.

(define (A072781 n) (A025581 (A072734 n))) ;; X-projection of A072735
(define (A072782 n) (A002262 (A072734 n))) ;; Y-projection of A072735

(define (A072783 n) (- (A072740 n) (A072736 n)))
(define (A072784 n) (- (A072741 n) (A072737 n)))
(define (A072785 n) (- (A072781 n) (A072738 n)))
(define (A072786 n) (- (A072782 n) (A072739 n)))


;; (map binexp->runcount1list (iota0 16))
;; --> (() (1) (1 1) (2) (1 2) (1 1 1) (2 1) (3) (1 3)
;;      (1 2 1) (1 1 1 1) (1 1 2) (2 2) (2 1 1) (3 1) (4) (1 4))

(define (binexp->runcount1list n) ;; (length (binexp->runcount1list n)) gives A005811
   (if (zero? n)
       (list)
       (let loop ((n n) (rc (list)) (count 0) (prev-bit (modulo n 2)))
          (if (zero? n)
              (cons count rc)
              (if (eq? (modulo n 2) prev-bit) ;; If the lsb has not changed
                  (loop (floor->exact (/ n 2))
                        rc
                        (1+ count)
                        (modulo n 2)
                  )
                  (loop (floor->exact (/ n 2)) ;; If it has changed
                        (cons count rc)
                        1
                        (modulo n 2)
                  )
              )
          )
       )
  )
)


;; Replace in the binary expansion of n each 1-bit
;; with the distance to the next 1-bit:
;;
;; (binexp->siteswap-list 5)  --> (2 0 1)
;; (binexp->siteswap-list 6)  --> (1 2 0)
;; (binexp->siteswap-list 7)  --> (1 1 1)
;; (binexp->siteswap-list 8)  --> (4 0 0 0)
;; (binexp->siteswap-list 9)  --> (3 0 0 1)
;; (binexp->siteswap-list 10) --> (2 0 2 0)

(define (binexp->siteswap-list n)
   (if (zero? n)
       (list)
       (let loop ((n n) (rc (list)) (count 1))
          (if (zero? n)
              rc
              (if (zero? (modulo n 2))
                  (loop (/ n 2)
                        (cons 0 rc)
                        (1+ count)
                  )
                  (loop (floor->exact (/ n 2))
                        (cons count rc)
                        1
                  )
              )
          )
       )
  )
)


;; CatTrianglDirect := (r,m) -> `if`((m < 0),0,((r-m+1)*(r+m)!)/(r! * m! * (r+1)));
;; A009766 gives also:
;; a(n,m) = C(n+m,n)*(n-m+1)/(n+1), n >= m >= 0.

(define (CatTriangle r m)
   (if (< m 0)
       0
       (/ (* (1+ (- r m)) (! (+ r m)))
          (* (! r) (! m) (1+ r))
       )
   )
)


(define (self-convolve fun upto_n)
  (+ (if (zero? (modulo upto_n 2)) (expt (fun (/ upto_n 2)) 2) 0)
     (* 2
        (let loop ((i 0) (j upto_n))
          (if (>= i j)
              0
              (+ (* (fun i) (fun j))
                 (loop (1+ i) (-1+ j))
              )
          )
        )
     )
  )
)


(definec (A000108 n) ;; Catalan's numbers.
   (if (zero? n)
       1
       (self-convolve A000108 (-1+ n))
   )
)

(definec (A014137 n) ;; Partial sums of A000108: 1,2,4,9,23,65,197,...
   (if (zero? n)
       1
       (+ (A014137 (-1+ n)) (A000108 n))
   )
)

(define (first_pos_with_funs_val_gte fun n)
   (let loop ((i 0))
      (if (>= (fun i) n) i
          (loop (1+ i))
      )
   )
)

;; (map A072643 (iota0 23)) ;; n occurs (A000108 n) times.
;; --> (0 1 2 2 3 3 3 3 3 4 4 4 4 4 4 4 4 4 4 4 4 4 4 5)

;; Was defined with definec, but without reason:

(define (A072643 n) (first_pos_with_funs_val_gte A014137 (1+ n))) ;; Was ranks_w/2