;
; a066425s.txt  --  This file contains the Scheme-functions
; that compute the integer sequences for the OEIS sequences
; A066425, A068054 - A068059, A068221 - A068224.
;
; The associated OEIS entry can be found with the URL:
; http://www.research.att.com/cgi-bin/access.cgi/as/njas/sequences/eisA.cgi?Anum=A066425
;
; The code here is still exponential, although significantly faster in
; practice than Maple code given at A066425.
;
; To compute even higher values, use either a different implementation
; of Scheme, implement the same algorithm in C, using arbitrary
; precision integer library like Gnu MP (GMP)
; (see http://www.gnu.org/software/gmp/gmp.html
; or http://www.swox.com/gmp/ for that)
; or use some insight/heuristic to reduce the memory
; usage significantly from the current A068739(n) bits.
;
;

; Coded by Antti Karttunen (E-mail: my_firstname.my_surname(AT)iki.fi)
; in February 2002. These functions are in public domain.


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;           Auxiliary functions for testing & output                 ;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

(define output_seq
  (lambda (seq)
      (cond ((null? seq) (newline))
            (else (write (car seq))
                  (write-string ",")
                  (output_seq (cdr seq))
            )
      )
  )
)

(define reversed_iota
  (lambda (n)
      (if (zero? n) ()
          (cons n (reversed_iota (-1+ n)))
      )
  )
)

(define iota (lambda (n) (reverse! (reversed_iota n))))


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;    Macro defunac for defining cached unary (integer) functions     ;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;


; See http://www.swiss.ai.mit.edu/projects/scheme/mit/macros.txt
; why this is needed in MIT Scheme, releases before 7.7.0:

(define-macro (define-macro-both pattern . body)
  `(begin
     (define-macro ,pattern ,@body)
     (syntax-table-define user-initial-syntax-table ',(car pattern)
       (macro ,(cdr pattern)
	 ,@body))))

(syntax-table-define user-initial-syntax-table 'define-macro-both
  (macro (pattern . body)
    `(begin
       (define-macro ,pattern ,@body)
       (syntax-table-define system-global-syntax-table ',(car pattern)
	 (macro ,(cdr pattern)
	   ,@body)))))



(define vector-set-and-return-value!
        (lambda (vec ind val) (vector-set! vec ind val) val)
)

; define unary cached functions. Syntax is like
; defun of Lisp (without any keywords):
; (defunac funcname (one_arg) body...)
; Note that one should never use argument named
; _cache_ here!
; 

(define-macro-both (defunac funam . a_et_b)
  `(define ,funam
     (letrec
        ((_cache_ (make-vector 16))
         (,funam
          (lambda ,(car a_et_b)
            (cond ((null? ,(caar a_et_b)) _cache_)
                  (else
                      (if (fix:>= ,(caar a_et_b) (vector-length _cache_))
                          (set! _cache_
                                (vector-grow _cache_
                                             (max (1+ ,(caar a_et_b))
                                                  (* 2 (vector-length _cache_))
                                             )
                                )
                          )
                      )
                      (or (vector-ref _cache_ ,(caar a_et_b))
                          (vector-set-and-return-value! _cache_ ,(caar a_et_b)
                                                                ,@(cdr a_et_b))
                      )
                  )
            ) ; cond
          )
         )
        ) ; letrec-definitions
       ,funam
     ) ; letrec
   )
)


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;       Simple transformations: DIFF, PSUM and DISTSUMS              ;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

(define PSUMaux
   (lambda (a b sum)
        (if (null? a) b
            (PSUMaux (cdr a)
                     (cons (+ (car a) sum) b)
                     (+ (car a) sum)
            )
        )
   )
)

(define PSUM (lambda (a) (reverse! (PSUMaux a () 0))))

(define DIFF (lambda (a) (map - (cdr a) (except-last-pair a))))
; (define DIFF (lambda (a) (map - (cdr a) (list-head a (-1+ (length a))))))

(define sum_elems_by_binexp
   (lambda (v c s)
      (cond ((zero? c) s)
            ((zero? (fix:and c 1))
               (sum_elems_by_binexp (cdr v) (fix:lsh c -1) s))
            (else
               (sum_elems_by_binexp (cdr v) (fix:lsh c -1) (+ s (car v))))
      )
   )
)


(define DistSumsTransformAux
   (lambda (v z i lim)
      (if (> i lim)
          z
          (DistSumsTransformAux v
                                (cons (sum_elems_by_binexp v i 0) z)
                                (1+ i)
                                lim
          )
      )
   )
)

(define DISTSUMS
   (lambda (v)
     (reverse! (DistSumsTransformAux v () 0 (-1+ (expt 2 (length v)))))))


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;      Computing A066425 via mask sequences A068221 & A068222        ;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;


(define (binwidth_does_not_work_with_large_n n)
  (if (zero? n) 0 (1+ (floor->exact (/ (log n) (log 2))))))

(define (binwidth_aux n i)
  (if (zero? n) i (binwidth_aux (floor->exact (/ n 2)) (1+ i)))
)

(define (binwidth n) (binwidth_aux n 0))


(define (pos_of_first_zero_bit_aux n w)
   (bit-substring-find-next-set-bit (bit-string-not (unsigned-integer->bit-string w n)) 0 w)
)


(define (pos_of_first_zero_bit n) (pos_of_first_zero_bit_aux n (1+ (binwidth n))))

(define (or_with_right_shifted_copy m i)
  (let* ((w (binwidth m))
         (bs1 (unsigned-integer->bit-string w m))
         (bs2 (unsigned-integer->bit-string w (floor->exact (/ m (expt 2 i)))))
        )
    (bit-string-or! bs1 bs2)
    (bit-string->unsigned-integer bs1)
  )
)

(define (or_a_with_right_shifted_b a b i)
  (let* ((w (binwidth a))
         (bs1 (unsigned-integer->bit-string w a))
         (bs2 (unsigned-integer->bit-string w (floor->exact (/ b (expt 2 i)))))
        )
    (bit-string-or! bs1 bs2)
    (bit-string->unsigned-integer bs1)
  )
)


(defunac A068058 (n) (pos_of_first_zero_bit (A068222 n)))

(defunac A068059 (n) (if (< n 2) n (+ (A068059 (-1+ n)) (A068058 (-1+ n)))))

(defunac A068221 (n)
  (if (< n 2) n
      (* (1+ (expt 2 (+ (binwidth (A068221 (-1+ n)))
                        (-1+ (A068059 (-1+ n)))
                     )
             )
         )
         (A068221 (-1+ n))
      )
  )
)


(defunac A068222 (n)
 (if (< n 2) n
   (or_a_with_right_shifted_b (A068221 n) (A068222 (-1+ n)) (A068058 (-1+ n)))
 )
)


(defunac A066425 (n) (if (< n 2) n (+ (* 2 (A066425 (-1+ n))) (A068058 (-1+ n)))))


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; Compute A066425 with procedures A066425compute_way1 & A066425compute_way2 ;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;




(define bit-string-set-multiple!
  (lambda (bv indices) (for-each (lambda (i) (bit-string-set! bv i)) indices) bv)
)


; Occupy bits: 0,1 (+0) & 1,2, (+1) the next free: 3, the next distinct sums: 3,4
; Thus occupy bits 3,4 & 3+1,4+1 & 3+3,4+3, next free: 8, the next distinct sums: 8,9,11,12
; Thus occupy bits (+0) 8,9,11,12 & (+1) 9,10,12,13 & (+3) 11,12,14,15 & (+8) 16,17,19,20,
; next free is 18: the next subset of distinct sums (= nesudisu)
; is 18,19,21,22,26,27,29,30 (+0), which are occupied
; &  19,20,22,23,27,28,30,31 (+1)
; &  21,22,24,25,29,30,32,33 (+3)
; &  26,27,29,30,34,35,37,38 (+8)
; &  36,37,39,40,44,47,49,48 (+18)
; so the next free is 41.


(define A066425compute_way1
   (lambda (upto_n)
     (A066425aux_way1 0 upto_n (list 0) (list 0)
         ((lambda (bv) (bit-string-set! bv 0) bv)
             (make-bit-string (expt 2 (+ upto_n 3)) #f))
     )
   )
)


(define A066425compute_way2
   (lambda (upto_n)
     (A066425aux_way2 0 upto_n (list 0) (list 0)
         ((lambda (bv) (bit-string-set! bv 0) bv)
             (make-bit-string (expt 2 (+ upto_n 3)) #f))
     )
   )
)



; A straightforward, naive approach:

(define bit-string-mark-next-positions-in-loop!
     (lambda (bv next results nesudisu)
        (for-each
          (lambda (r)
             (for-each (lambda (d) (bit-string-set! bv (fix:+ r d))) nesudisu)
          )
          (cons next results)
        )
     )
)


;
; This writes most of the bits en bloc:
;

(define bit-string-mark-next-positions!
     (lambda (bv next results nesudisu)
        (bit-substring-move-right! bv 0  (1+ (+ (car (last-pair nesudisu)) next)) bv next)
        (for-each (lambda (d) (bit-string-set! bv (fix:+ next d))) nesudisu)
     )
)

; (load-option 'format) To use format, do this.


(define A066425aux_way1
   (lambda (i upto_n results distsums bv)
     (cond ((eq? i upto_n) (cdr (reverse! results)))
           (else
          (write results) (newline)
        (let* ((next (bit-substring-find-next-set-bit
                        (bit-string-not bv) (car results) (bit-string-length bv)))
;              (nesudisu (map + (make-list (expt 2 i) next) distsums))
               (nesudisu (map (lambda (n) (fix:+ n next)) distsums))
              )
;          (format #t "bv=~S~%next=~S, results=~S, distsums=~S~%nesudisu=~S~%"
;                       bv next results distsums nesudisu)
          (bit-string-mark-next-positions-in-loop! bv next results nesudisu)

          (A066425aux_way1 (1+ i) upto_n (cons next results)
                             (append! distsums nesudisu) bv
          )
        )
     ))
   )
)

(define A066425aux_way2
   (lambda (i upto_n results distsums bv)
     (cond ((eq? i upto_n) (cdr (reverse! results)))
           (else
          (write results) (newline)
        (let* ((next (bit-substring-find-next-set-bit
                        (bit-string-not bv) (car results) (bit-string-length bv)))
;              (nesudisu (map + (make-list (expt 2 i) next) distsums))
               (nesudisu (map (lambda (n) (fix:+ n next)) distsums))
              )
;          (format #t "bv=~S~%next=~S, results=~S, distsums=~S~%nesudisu=~S~%"
;                       bv next results distsums nesudisu)
          (bit-string-mark-next-positions! bv next results nesudisu)

          (A066425aux_way2 (1+ i) upto_n (cons next results)
                             (append! distsums nesudisu) bv
          )
        )
     ))
   )
)

; Do:
; set MITSCHEME_LARGE_HEAP=5000
; before starting Scheme, if you want to compute this
; upto the 21st term with (A066425compute_way2 21)

(define A066425seq '(1 3 8 18 41 84 181 364 751 1512 3037 6107 12216 24547 49117 98236 196544 393178 786407 1573201 3146426))
(newline) (write-string "A066425=") (output_seq A066425seq)

(define A068054seq (DIFF A066425seq))
(write-string "A068054=") (output_seq A068054seq)

(define A068055seq (PSUM A066425seq))
(write-string "A068055=") (output_seq A068055seq)

(define A068056seq (DISTSUMS (list-head A066425seq 7)))
(write-string "A068056=") (output_seq A068056seq)

(define A068057seq (DIFF A068056seq))
(write-string "A068057=") (output_seq A068057seq)

(define double (lambda (n) (* 2 n)))

; (define A068058seq (map A068058 (iota 14)))
(define A068058seq (map -  (cdr A066425seq) (map double (except-last-pair A066425seq))))
(write-string "A068058=") (output_seq A068058seq)

; (define A068059seq (map A068059 (iota 14)))
(define A068059seq (cons 1 (map -  (cdr A066425seq) (except-last-pair A068055seq))))
(write-string "A068059=") (output_seq A068059seq)

(define A068221seq (map A068221 (iota 10)))
(write-string "A068221=") (output_seq A068221seq)

(define A068222seq (map A068222 (iota 10)))
(write-string "A068222=") (output_seq A068222seq)



; A066425=1,3,8,18,41,84,181,364,751,1512,3037,6107,12216,24547,49117,98236,196544,393178,786407,1573201,3146426
; A068054=2,5,10,23,43,97,183,387,761,1525,3070,6109,12331,24570,49119,98308,196634,393229,786794,1573225
; A068055=1,4,12,30,71,155,336,700,1451,2963,6000,12107,24323,48870,97987,196223,392767,785945,1572352,3145553,6291979
; A068056=0,1,3,4,8,9,11,12,18,19,21,22,26,27,29,30,41,42,44,45,49,50,52,53,59,60,62,63,67,68,70,71,84,85,87,88,92,93,95,96,102,103,105,106,110,111,113,114,125,126,128,129,133,134,136,137,143,144,146,147,151,152,154,155,181,182,184,185,189,190,192,193,199,200,202,203,207,208,210,211,222,223,225,226,230,231,233,234,240,241,243,244,248,249,251,252,265,266,268,269,273,274,276,277,283,284,286,287,291,292,294,295,306,307,309,310,314,315,317,318,324,325,327,328,332,333,335,336
; A068057=1,2,1,4,1,2,1,6,1,2,1,4,1,2,1,11,1,2,1,4,1,2,1,6,1,2,1,4,1,2,1,13,1,2,1,4,1,2,1,6,1,2,1,4,1,2,1,11,1,2,1,4,1,2,1,6,1,2,1,4,1,2,1,26,1,2,1,4,1,2,1,6,1,2,1,4,1,2,1,11,1,2,1,4,1,2,1,6,1,2,1,4,1,2,1,13,1,2,1,4,1,2,1,6,1,2,1,4,1,2,1,11,1,2,1,4,1,2,1,6,1,2,1,4,1,2,1
; A068058=1,2,2,5,2,13,2,23,10,13,33,2,115,23,2,72,90,51,387,24
; A068059=1,2,4,6,11,13,26,28,51,61,74,107,109,224,247,249,321,411,462,849,873
; A068221=1,3,27,6939,1819024155,4000076419257644882715,77372730618755190082837640023742562619343248155,237146729315719909074685533605289611245100501015339264129426202250460094330698830235303562543675218715,8911187074616127159630138357075376547105533627449036419992725829026420755813478002181762031684518586027031663576818914678933298957191317061964541961930290219206812943157985636166985187949005447441861847007632155,105550988443128699847561088957135066295771053152870985894051303513623470170918963309279940118635455494292568456350580230255605462968748315430207607329505017449556550512448292542860165311056301669951619199191289929709640691578515598784709391754840066514585236587762271315791884648864972535118380359104009269860190472736832178919823530064213646657225566244885486463766625834181687216546803046519573357361745907412726523803113172963085261595
; A068222=1,3,27,6943,1819024347,4000076419257964896255,77372730618755190082837640024229164347468807163,237146729315719909074685533605289611245100501015339264143038071525981840734161007646995989118658281471,8911187074616127159630138357075376547105533627449036419992725829026420755813478002181762031684518586027031663576818939820031318888851486918027744750598302906858440892925063402349886288923546289902149182064417791,105550988443128699847561088957135066295771053152870985894051303513623470170918963309279940118635455494292568456350580230255605462968748315430207607329505017449556550512448292542860165311056301669951619199191289929709640691578515604884743906920666741158615889926739229770783146862432879234580138971689969096225678794080782492012099258205618338044437946244495467114728249912695756080000435289942855794230024102257204681848384210463855992831

; A068223 (= A068221 in binary) & A068224 (= A068222 in binary)