;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;                                                                      ;;
;;  http://www.iki.fi/~kartturi/matikka/Nekomorphisms/gatorank.scm      ;;
;;  - Functions for ranking & unranking objects in Catalan families,    ;;
;;    in the standard lexicographical order (A014486)                   ;;
;;                                                                      ;;
;;  This Scheme-code is Copyright (C) 2002 by Antti Karttunen           ;;
;;  (E-mail: my_firstname.my_surname@iki.fi) and is placed under        ;;
;;  the GPL (Gnu Public License), so you are free to copy it.           ;;
;;                                                                      ;;
;;  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/                        ;;
;;                                                                      ;;
;;  Aubrey Jaffer's SLIB Scheme library is available at:                ;;
;;  http://www.swiss.ai.mit.edu/~jaffer/SLIB.html                       ;;
;;                                                                      ;;
;;  The main pointer for this code collection is:                       ;;
;;  http://www.iki.fi/~kartturi/matikka/Nekomorphisms/gatomorf.htm      ;;
;;                                                                      ;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;


;;
;; 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\\Nekomorphisms\\gatorank" "c:\\matikka\\Nekomorphisms\\")
;; )
;; 

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


;; (map A014486 (cons 0 (iota 23)))
;; --> (0 2 10 12 42 44 50 52 56 170 172 178 180 184 202 204 210 212 216 226 228 232 240 682)

(definec (A014486 n) ;; was definec
   (let ((w/2 (A072643 n)))
      (CatalanUnrank w/2 (if (zero? n) 0 (- n (A014137 (-1+ w/2)))))
   )
)


;; (map safe_w/2 (cons 0 (iota 42)))
;; --> (0 2 3 3 3 4 4 4 4 4 4 4 4 4 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 6)
(define (safe_w/2 n) (first_pos_with_funs_val_gte A000108 (1+ n)))


(define (CatalanRankLocal a)
   (if (zero? a)
       0
       (+ (CatalanRank (/ (binwidth a) 2) a))
   )
)

(define (CatalanRankGlobal a)
   (if (zero? a)
       0
       (let ((w/2 (/ (binwidth a) 2)))
          (+ (A014137 (-1+ w/2))
             (CatalanRank w/2 a)
          )
       )
   )
)

;; This should produce same as (cons 0 (iota 6919)):
;;
;; (map CatalanRankGlobal
;;      (map parenthesization->binexp
;;           (map binexp->parenthesization
;;                (map A014486 (cons 0 (iota 6919))))))
;; 

;; See http://www.iki.fi/~kartturi/matikka/tab9766.htm


(define (CatalanRank w/2 a)
    (let loop ((a a) ;; The totally balanced binary expansion
               (r 0)
               (lo 0)
               (y -1)
              )
      (if (zero? a)
          (- (/ (C (* 2 w/2) w/2) (1+ w/2))
             (1+ lo)
          )
          (if (zero? (modulo a 2))
              (loop ;; Down to the valley
                    (floor->exact (/ a 2)) ;; Was: (fix:lsh a -1) ;; a >>= 1
                    (1+ r)
                    (+ lo (CatTriangle (1+ r) y))
                    y
              )
              (loop ;; Upto the mountain high.
                    (floor->exact (/ a 2)) ;; Was: (fix:lsh a -1)
                    r lo (1+ y)
              )
          )
      )
    )
)

(define (CatalanUnrank w/2 orank)
  (let ((rank (- (/ (C (* 2 w/2) w/2) (1+ w/2))
                 (1+ orank)
              )
        )
       )
    (let loop ((a 0)         ;; Constructed bit-string
               (lo 0)
               (t w/2)       ;; The row on A009766
               (y (-1+ w/2)) ;; The position on row t of A009766
               (m (CatTriangle w/2 (-1+ w/2)))
              )
      (if (zero? t) a
          (if (> (+ lo m) rank)
              (loop (1+ (* 2 a)) ;; Up the mountain high
                    lo
                    t
                    (-1+ y)
                    (CatTriangle t (-1+ y))
              )
              (loop (* 2 a)      ;; Down to the valley low
                    (+ lo m)
                    (-1+ t)
                    y
                    (CatTriangle (-1+ t) y)
              )
          )
      )
    )
  )
)


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;                                                                    ;;
;; Function parenthesization->binexp and two versions of              ;;
;;  binexp->parenthesization, first the                               ;;
;;   straightforward version converted from Maple code,               ;;
;;    and then the more enlightened "Forth"-inspired version.         ;;
;;                                                                    ;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;


;; (parenthesization->binexp '())      --> 0
;; (parenthesization->binexp '(()))    --> 2
;; (parenthesization->binexp '(() ())) --> 10
;; (parenthesization->binexp '((())))  --> 12

(define (parenthesization->binexp p)
   (let loop ((s 0) (p p))
      (if (null? p)
          s
          (let* ((x (parenthesization->binexp (car p)))
                 (w (binwidth x))
                )
             (loop
;; Do not use!     (+ (fix:lsh s (+ w 2)) (fix:lsh 1 (1+ w)) (* 2 x))
                   (+ (* s (expt 2 (+ w 2))) (expt 2 (1+ w)) (* 2 x))
                   (cdr p)
             )
          )
      )
   )
)


;; PeelNextBalSubSeq and RestBalSubSeq expect their
;; integer argument nn to contain the binary expansion
;; of the (sub-)parenthesization in reverse order,
;; with the least significant bit being always 1.


(define (PeelNextBalSubSeq nn) ;; We assume that given nn is odd.
  (let loop ((z 0) (level -1) (n (fix:lsh nn -1)))
      (cond ((zero? level) (/ z 2)) ;; n on prev. iteration must has been even.
            (else (loop (+ (fix:lsh z 1) (modulo n 2)) ;; z <<= 1, z += n % 2
                        (+ level (expt -1 (modulo n 2)))
                        (fix:lsh n -1) ;; n >>= 1
                  )
            )
      )
  )
)


(define (RestBalSubSeq nn) ;; We assume that given nn is odd.
  (let loop ((level -1) (n (fix:lsh nn -1)))
      (cond ((zero? level) (PeelNextBalSubSeq (1+ (* 2 n))))
            (else (loop (+ level (expt -1 (modulo n 2)))
                        (fix:lsh n -1) ;; n >>= 1
                  )
            )
      )
  )
)


(define (ConsTBBS a b) ;; "cons" two totally balanced binary sequences
   (let ((aw (binwidth a))
         (bw (binwidth b))
        )
     (+ (expt 2 (+ 1 aw bw)) (* a (expt 2 (1+ bw))) b)
   )
)


(define (reversed_binexp->parenthesization n)
   (cons (binexp->parenthesization_in_dumb_way (PeelNextBalSubSeq n))
         (binexp->parenthesization_in_dumb_way (RestBalSubSeq n))
   )
)

(define (binexp->parenthesization_in_dumb_way n)
   (if (zero? n)
       (list)
       (reversed_binexp->parenthesization (binrev n))
   )
)


;; Now, if we remember that "Lisp" spelled backwards is "Forth",
;; and the parenthesizations have another form as Dyck paths,
;; it's much easier to implement this by scanning the totally
;; balanced binary string from the end (the rightmost = the least
;; significant bit) to the beginning (to the leftmost = the most
;; significant bit). Note how we don't need double-forked
;; recursion anymore, but just simple tail-recursion is enough.
;;


;; Sep 06 2002:  (fix:lsh n -1) changed to (floor->exact (/ n 2))

(define (binexp->parenthesization n)
   (let loop ((n n) (stack (list (list))))
       (cond ((zero? n) (car stack))
             ((zero? (modulo n 2))
                (loop (floor->exact (/ n 2)) (cons (list) stack))
             )
             (else
                (loop (floor->exact (/ n 2)) (cons2top! stack))
             )
       )
   )
)



;; Our own copy of p->Lw, copied from gatocout.scm:

(define (p->Lw p)
  (reverse! (cdr (reverse
    (let recurse ((p p))
           (cond ((null? p) (list 0))
                 (else ;; it is a list.
                    (append! (list (length p))
                             (apply append! (map recurse p))
                    )
                 )
           )
    )
  )))
)

;; The Lw->parenthesization and p->Lw (transferred to gatocout.scm)
;; are inverses of each other:
;; (first-dislocated (apply_upto_n 2056 (compose-funlist (list Lw->parenthesization p->Lw))))
;; --> ()

(define (Lw->parenthesization L)
   (let loop ((L (reverse L)) (stack (list (list)))) ;; The last leaf is implicit
       (cond ((null? L) (car stack))
             (else
                (loop (cdr L) (list-n-from-top (car L) stack))
             )
       )
   )
)



;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;                                                                    ;;
;;           Transferred from gatomain.scm                            ;;
;;                                                                    ;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;



;; (BinTree2Tree '(() . (() . ()))) --> (() (() ()))
;; (BinTree2Tree '((() . ()) . ())) --> ((() ()) ())


(define (BinTree2Tree bt) ;; Not a bijection. (A057123)
  (cond ((not (pair? bt)) bt)
        (else (list (BinTree2Tree (car bt))
                    (BinTree2Tree (cdr bt)))
        )
  )
)


(define (Tree2BinTree_if_possible gt)
  (call-with-current-continuation
    (lambda (e)
     (let recurse ((gt gt))
        (cond ((not (pair? gt)) gt)
              ((eq? 2 (length gt))
                    (cons (recurse (car gt))
                          (recurse (cadr gt))
                    )
              )
              (else (e '()))
        )
     ) ;; let
    )
  )
)

;; Only the (define A057548 (catfun1 list)) will then erroneously
;; give 0, instead of 1, if we initialize the cache like this:

(define (new-cat-cache)
   (vector 0 #f #f #f #f #f #f #f #f #f #f #f #f #f #f #f #f #f #f #f #f #f)
)


(define (catfun0 morphism) ;; only for telescoping (self-embeddable) automorphisms,
;; like df->bf or deeprevoncarside...
     (letrec ((cache (new-cat-cache))
              (fun (lambda (y)
                 (cond ((null? y) cache) ; For debugging, reveal our cache for nothing.
                       (else
                           (if (fix:>= y (vector-length cache))
                               (set! cache (vector-grow cache
                                             (max (1+ y)
                                                  (* 2 (vector-length cache))
                                             )
                                           )
                               )
                           )
                           (or (vector-ref cache y)
                               (vector-set-and-return-value! cache y
                                     (CatalanRankLocal
                                       (parenthesization->binexp
                                         (morphism
                                           (binexp->parenthesization
                                             (CatalanUnrank (safe_w/2 y) y)
                                           )
                                         )
                                       )
                                     )
                               )
                           )
                       )
                 )
                )
              )
             )
         fun
     )
)


(define (catfun1 morphism)
     (letrec ((cache (new-cat-cache))
              (fun (lambda (y)
                 (cond ((null? y) cache) ; For debugging, reveal our cache for nothing.
                       (else
                           (if (fix:>= y (vector-length cache))
                               (set! cache (vector-grow cache
                                             (max (1+ y)
                                                  (* 2 (vector-length cache))
                                             )
                                           )
                               )
                           )
                           (or (vector-ref cache y)
                               (vector-set-and-return-value! cache y
                                     (CatalanRankGlobal
                                       (parenthesization->binexp
                                         (morphism
                                           (binexp->parenthesization
                                             (A014486 y)))))
                               )
                           )
                       )
                 )
                )
              )
             )
         fun
     )
)

(define (catfun1restricted_to_bt_subset morphism)
     (letrec ((morphism_restricted_to_bt_subset
                (lambda (s)
                    (Tree2BinTree_if_possible (morphism (BinTree2Tree s)))
                )
              )
              (cache (new-cat-cache))
              (fun (lambda (y)
                 (cond ((null? y) cache) ; For debugging, reveal our cache for nothing.
                       (else
                           (if (fix:>= y (vector-length cache))
                               (set! cache (vector-grow cache
                                             (max (1+ y)
                                                  (* 2 (vector-length cache))
                                             )
                                           )
                               )
                           )
                           (or (vector-ref cache y)
                               (vector-set-and-return-value! cache y
                                     (CatalanRankGlobal
                                       (parenthesization->binexp
                                         (morphism_restricted_to_bt_subset
                                           (binexp->parenthesization
                                             (A014486 y)))))
                               )
                           )
                       )
                 )
                )
              )
             )
         fun
     )
)

;; Like previous, but there's an extra stem beneath the binary tree,
;; i.e. it's a proper trivalent-plane-tree:
(define (catfun1restricted_to_tpt_subset morphism)
     (letrec ((morphism_restricted_to_tpt_subset
                (lambda (s)
                  (Tree2BinTree_if_possible (car (morphism (list (BinTree2Tree s)))))
                )
              )
              (cache (new-cat-cache))
              (fun (lambda (y)
                 (cond ((null? y) cache) ; For debugging, reveal our cache for nothing.
                       (else
                           (if (fix:>= y (vector-length cache))
                               (set! cache (vector-grow cache
                                             (max (1+ y)
                                                  (* 2 (vector-length cache))
                                             )
                                           )
                               )
                           )
                           (or (vector-ref cache y)
                               (vector-set-and-return-value! cache y
                                     (CatalanRankGlobal
                                       (parenthesization->binexp
                                         (morphism_restricted_to_tpt_subset
                                           (binexp->parenthesization
                                             (A014486 y)))))
                               )
                           )
                       )
                 )
                )
              )
             )
         fun
     )
)

(define (catfun2 par2int-fun)
     (letrec ((cache (new-cat-cache))
              (fun (lambda (y)
                 (cond ((null? y) cache) ; For debugging, reveal our cache for nothing.
                       (else
                           (if (fix:>= y (vector-length cache))
                               (set! cache (vector-grow cache
                                             (max (1+ y)
                                                  (* 2 (vector-length cache))
                                             )
                                           )
                               )
                           )
                           (or (vector-ref cache y)
                               (vector-set-and-return-value! cache y
                                         (par2int-fun
                                           (binexp->parenthesization
                                             (A014486 y)
                                           )
                                         )
                               )
                           )
                       )
                 )
                )
              )
             )
         fun
     )
)

(define (apply_upto_n upto_n morphism) ;; For testing.
   (map CatalanRankGlobal
       (map parenthesization->binexp
            (map morphism
                 (map binexp->parenthesization
                      (map A014486 (iota0 upto_n))))))
)


(define (binseqs_of_size size)
    (map (lambda (r) (CatalanUnrank size r))
         (iota0 (-1+ (A000108 size)))
    )
)

(define (partition_by_gatomorphism size gatomorphism)
   (let ((src_set (map binexp->parenthesization (binseqs_of_size size)))
        )
     (let loop ((cur (car src_set)) (src src_set) (res (list (list))))
        (cond ((null? src) (reverse! (map reverse! res)))
              ((member cur src)
                 (loop (gatomorphism cur)
                       (delete! cur src)
                       (cons (cons cur (car res)) (cdr res))
                 )
              )
              (else
;; Completed one whole cycle, let's begin the next one with the first
;; parenthesization we have left:
                 (loop (car src) src (cons (list) res))
              )
        ) ; cond
     ) ; let loop
   )
)


;; Return a Cycle-Count function for a particular gatomorphism.
;; (This one is quite ineffective. See the module gatosiga.scm for a better one.)

(define (cc-fun gatomorphism)
  (lambda (n)
    (length (partition_by_gatomorphism n gatomorphism))
  )
)


(define (fixed-by-gatomorphism size gatomorphism)
   (keep-matching-items
        (map binexp->parenthesization (binseqs_of_size size))
        (lambda (p) (equal? (gatomorphism p) p))
   )
)


(define (indices-of-nth-forest n)
   (map (lambda (x) (+ (A014137 (-1+ n)) x)) (iota0 (-1+ (A000108 n))))
)

(define (partition-by-gatoAfun size gatomorphism)
   (let ((src_set (indices-of-nth-forest size)))
     (let loop ((cur (car src_set)) (src src_set) (res (list (list))))
        (cond ((null? src) (reverse! (map reverse! res)))
              ((member cur src)
                 (loop (gatomorphism cur)
                       (delete! cur src)
                       (cons (cons cur (car res)) (cdr res))
                 )
              )
              (else
;; Completed one whole cycle, let's begin the next one with the first
;; parenthesization we have left:
                 (loop (car src) src (cons (list) res))
              )
        ) ; cond
     ) ; let loop
   )
)



(define (number-of-1-cycles cycles)
  (let ((fes 0))
    (for-each (lambda (c)
                (if (and (pair? c) (not (pair? (cdr c)))) (set! fes (1+ fes)))
              )
              cycles
    )
    fes
  )
)


(define (num-of-ones cycles)
  (let ((fes 0))
    (for-each (lambda (c) (if (= 1 c) (set! fes (1+ fes))))
              cycles
    )
    fes
  )
)

(define (fc-Afun Afun) (lambda (n) (number-of-1-cycles (partition-by-gatoAfun n Afun))))

(define (cc-Afun Afun) (lambda (n) (length (partition-by-gatoAfun n Afun))))

(define (mc-Afun Afun)
  (lambda (n) (apply max (map length (partition-by-gatoAfun n Afun))))
)


(define (lc-Afun Afun)
  (lambda (n) (apply lcm (map length (partition-by-gatoAfun n Afun))))
)





(define (compute-and-print-count-seqs Afun outfile upto-n)
   (call-with-output-file outfile
     (lambda (outport)
       (let loop ((n 0)
                  (ccs (list 1))
                  (fcs (list 1))
                  (mcs (list 1))
                  (lcs (list 1))
                 )
          (format #t "n=~A: ccs=~A fcs=~A mcs=~A lcs=~A~%"
                      n ccs fcs mcs lcs
          )
          (format outport "n=~A: ccs=~A fcs=~A mcs=~A lcs=~A~%"
                      n ccs fcs mcs lcs
          )
          (flush-output outport)
          (cond
            ((< n upto-n)
                (let ((partlengths (map length (partition-by-gatoAfun (1+ n) Afun))))
                   (loop (1+ n)
                         (append ccs (list (length partlengths)))
                         (append fcs (list (num-of-ones partlengths)))
                         (append mcs (list (fold-left max 0 partlengths)))
                         (append lcs (list (fold-left lcm 1 partlengths)))
                   )
                )
            )
          )
       )
     )
   )
)



(define (next-left-branch L)
   (ReflectBinTree (list (gmA057164 L)))
)

(define (nth-branch n) ((compose-fun-to-nth-power next-left-branch n) '()))

(define (next-branch! branches)
  (attach! (next-left-branch (car branches)) branches)
)


(define (iterate-and-print-left-branches outfile upto-n)
   (call-with-output-file outfile
     (lambda (outport)
       (let loop ((n 0)
                  (branch (list))
                  (br-reversed (list))
                 )
          (format outport "n=~A: (length branch)=~A" n (length branch))
          (cond ((or (not (pair? branch)) (null? (car br-reversed)))
                  (format outport " ends with ()")
                  (cond ((equal? br-reversed branch)
                           (format outport " and is symmetric: ~A" branch)
                        )
                  )
                )
          )
          (newline outport)
          (flush-output outport)
          (cond
            ((< n upto-n)
                (let ((next-branch (ReflectBinTree! (list br-reversed))))
                   (loop (1+ n)
                         next-branch
                         (gmA057164 next-branch)
                   )
                )
            )
          )
       )
     )
   )
)