;
; A068009S.TXT  --  This file contains the Scheme-functions
; that compute the integer sequences for the OEIS sequences
; A068009 - A068013, A068030-A068045 & A068049.
; These sequences arose from the discussion on the sci.math
; newsgroup thread "Subsets of {1,2,3,...,n}"
; that was started by Bill Pet on Jan 18 2002,
; upto Dave Rusin's comment on Feb 05 2002.
; The whole thread can be found with the following URL:
; http://groups.google.com/groups?hl=en&threadm=36abe133.0201181401.410577fe%40posting.google.com&rnum=1&prev=/groups%3Fhl%3Den%26selm%3D36abe133.0201181401.410577fe%2540posting.google.com
;
; 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=A068009
;

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

; This scheme code runs (at least) in MIT Scheme Release 7.6.0,
; for which documentation and the pre-compiled binaries
; (for various OS's running in Intel x86 architecture)
; can be found under the URL:
; http://www.swiss.ai.mit.edu/projects/scheme/


; The functions here are based on the matrix multiplication
; method given by Will Self in the 13th article of the thread,
; which can be found with the URL:
; http://groups.google.com/groups?hl=en&selm=e9c74932.0201201826.36629677%40posting.google.com
; and which is copied here in verbatim:

;
; It seems that one really needs to answer all the questions
; 
; How many subsets of {1,2,...,k} sum to j modulo m?
; 
; for j = 0, 1, 2, ... , m-1  in order to approach this with recurrence
; relations. Here is a set of recurrence equations which can be put
; neatly in vector form.
; 
; Let I denote the m by m identity matrix, and let v denote the m by m
; matrix which has 1's on the subdiagonal and in the upper-right-hand
; corner, and 0's elsewhere.  (That is, v(1, m) = 1 and v(p, p-1) = 1
; for 2 <= p <= m and v(p, q) = 0 otherwise.)  (That is, to get v,
; remove the bottom row of the identity matrix and put it on top.)  v is
; a circulant matrix, and v^m = I.
; 
; Let t(0) be the m-vector {1, 0, 0, ... , 0}.  For k >= 1, put
; 
; t(k) = (I + v^k) t(k-1).
; 
; Then the jth coordinate (count with j starting at 0, going up to m-1)
; of t(k) gives the number of subsets of {1,2,...,k} which sum to j
; modulo m.
; 
; This agrees with the results that Sophie gave.
; 
; To compute t(k) directly from t(0), you would need to multiply
; (I + v^k) ... (I + v^3) (I + v^2) (I + v)  by t(0).
; I don't know whether something more could be done with this matrix
; product.  Since v is a circulant matrix, the whole product is too.
; One might try to diagonalize.  The eigenvalues of v are the mth roots
; of unity, and the eigenvalues of the matrix product are
; (1 + a) (1 + a^2) (1 + a^3) ... (1 + a^k), where a is an eigenvalue of
; v.  So maybe some more could be done along these lines.
; 
; Will Self
;

; AK's note: By observing that the circulant matrix v used by Will
; is actually a permutation matrix of a simple m-cycle, we can
; implement this without using any matrix machinery at all,
; as long as we are not interested about a more sophisticated
; mathematical formula.
; Instead, we start from the t(0) vector {1, 0, 0, 0, ..., 0}
; and start summing it with its rotated copy, and iterate
; with progressively larger rotating steps. See the function ti
; given later.
; Notice the similarity to the definition of Shift Register Sequences
; (except that there we would use GF(2) addition, i.e. bitwise XOR
;  of the bit-vectors?)


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;     Auxiliary variables and functions for testing & output         ;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

(define from0to20 '(0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20))

(define output_seq
  (lambda (seq)
      (cond ((null? seq) ())
            (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))))

(define binomial_n_2 (lambda (n) (/ (* (-1+ n) n) 2)))

(define A025581; The X component (column) of square {0..inf} arrays
  (lambda (n) (- (binomial_n_2 (1+ (floor->exact (+ (/ 1 2) (sqrt (* 2 (1+ n))))))) (1+ n))))

(define A002262; The Y component (row) of square {0..inf} arrays
  (lambda (n) (- n (binomial_n_2 (floor->exact (+ (/ 1 2) (sqrt (* 2 (1+ n)))))))))

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


;(map binomial_n_2 from0to20)
;Value 14: (0 0 1 3 6 10 15 21 28 36 45 55 66 78 91 105 120 136 153 171 190)
;(map A025581 from0to20)
;Value 12: (0 1 0 2 1 0 3 2 1 0 4 3 2 1 0 5 4 3 2 1 0)
;(map A002262 from0to20)
;Value 13: (0 0 1 0 1 2 0 1 2 3 0 1 2 3 4 0 1 2 3 4 5)


; Given a binary function N(r,i), with r (the "row") ranging from
; 1 onward, and i (the index on the row) ranging from 1 to r,
; this "macro" will create a corresponding unary function
; which takes arguments from 1 onward, and gives the corresponding
; values of the triangle read row by row:


(define triangle_1_1
  (lambda (fun) (lambda (n) (fun (A002024 n) (1+ (A002262 (-1+ n))))))
)

; (define kolmio (triangle_1_1 cons))
; (map kolmio (cdr from0to20))
; ((1 . 1)
;  (2 . 1) (2 . 2)
;  (3 . 1) (3 . 2) (3 . 3)
;  (4 . 1) (4 . 2) (4 . 3) (4 . 4)
;  (5 . 1) (5 . 2) (5 . 3) (5 . 4) (5 . 5)
;  (6 . 1) (6 . 2) (6 . 3) (6 . 4) (6 . 5))
;


; Given a binary function N(r,c), with r (the "row", y) ranging from
; 1 to inf, and c (the column index on the row, x) ranging from 0 to inf,
; this "macro" will create a corresponding unary function
; which takes arguments from 0 onward, and gives the
; values of the corresponding infinite array read
; antidiagonal by antidiagonal:

(define table_1_0
  (lambda (fun) (lambda (n) (fun (1+ (A002262 n)) (A025581 n))))
)

; I.e.
; (define taulu (table_1_0 cons))
; (map taulu from0to20)
; --> ((1 . 0) (1 . 1) (2 . 0) (1 . 2) (2 . 1) (3 . 0) (1 . 3) (2 . 2) (3 . 1)
;      (4 . 0) (1 . 4) (2 . 3) (3 . 2) (4 . 1)
;      (5 . 0) (1 . 5) (2 . 4) (3 . 3) (4 . 2) (5 . 1) (6 . 0))
;


(define bin2unary_first_arg_fixed
  (lambda (fun first_arg) (lambda (n) (fun first_arg n)))
)



;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;        Functions implementing cyclic rotation of a vector          ;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

; When lista n grows to 31, starts whining that:
; The object 2147483648, passed as the second argument to general-car-cdr, is not in the correct range.

(define nthcdr_too_clever (lambda (n lista) (general-car-cdr lista (expt 2 n))))

(define nthcdr
   (lambda (n lista)
       (if (or (zero? n) (null? lista))
           lista
           (nthcdr (-1+ n) (cdr lista))
       )
   )
)


(define circularize! (lambda (lista) (append! lista lista)))

(define upto_previous
   (lambda (olp clp)
       (if (eq? olp (cdr clp))
           clp
           (upto_previous olp (cdr clp))
       )
   )
)

(define uncircularize!
   (lambda (lista)
       (set-cdr! (upto_previous lista lista) ())
       lista
   )
)

(define rol ; Rotate lista left c steps
   (lambda (c lista)
       (uncircularize! (nthcdr c (circularize! (list-copy lista))))
   )
)

; Rotate lista right c steps
(define ror (lambda (c lista) (rol (modulo (- c) (length lista)) lista)))

; Sum the "vector" v with its copy rotated right i steps:
(define sum_with_rotated_copy (lambda (v i) (map + v (ror i v))))

; The function ti computes the t(i):th m-vector, as occuring
; in Will Self's article, and where t(0) is the m-vector {1, 0, 0, ... , 0},
; and for i >= 1, t(i) = (I + v^i) t(i-1).
; I.e. as long as they fit, the successive rows of the funny table A053632
; (giving coefficients in expansion of Product_{k=1..n} (1+x^k))
; 1; 1,1; 1,1,1,1; 1,1,1,2,1,1,1; 1,1,1,2,2,2,2,2,1,1,1; 1,1,1,2,2,3,3,3,3,3,3,2,2,1,1,1;


(define ti
   (lambda (i m)
     (if (zero? i) (cons 1 (make-list (-1+ m) 0))
         (sum_with_rotated_copy (ti (-1+ i) m) i)
     )
   )
)

; Returns the number of subsets {1,...,i} that
; sum to 0 mod m, i.e. the first element of
; the corresponding vector computed with ti:

(define sim (lambda (m i) (car (ti i m))))

(define first_terms_greater_than_1 (lambda (m) (sim m (A002024 m)))) ; Gives A068049.

; Equal to A000009+1 and to A052839 (from the second term onward)
(define left_edge_of_A068049 (lambda (n) (sim (1+ (binomial_n_2 n)) n)))


; This kind of triangle not needed:
; (define A068009 (triangle_1_1 sim))
; (define A068009seq (map A068009 (iota (binomial_n_2 15))))
; (output_seq A068009seq)
; 2,1,2,1,2,4,1,1,2,4,1,1,2,4,8,1,1,2,3,6,12,1,1,1,3,5,10,20,1,1,1,2,4,8,16,32,1,1,1,2,4,8,15,30,60,1,1,1,2,4,7,13,26,52,104,1,1,1,1,3,6,12,24,47,94,188,1,1,1,1,3,6,11,22,44,86,172,344,1,1,1,1,2,5,10,20,39,79,158,316,632,1,1,1,1,2,5,10,19,37,73,147,293,586,1172,

; Use the table form instead:

(define A068009 (table_1_0 sim))

; (define A068009seq (map A068009 (cons 0 (iota (-1+ (binomial_n_2 15))))))
; (output_seq A068009seq)
; 1,2,1,4,1,1,8,2,1,1,16,4,2,1,1,32,8,4,1,1,1,64,16,6,2,1,1,1,128,32,12,4,2,1,1,1,256,64,24,8,4,2,1,1,1,512,128,44,16,8,3,1,1,1,1,1024,256,88,32,14,6,3,1,1,1,1,2048,512,176,64,26,12,5,2,1,1,1,1,4096,1024,344,128,52,22,10,4,2,1,1,1,1,8192,2048,688,256,104,44,20,8,4,2,1,1,1,1,


(define sim_row_1 (bin2unary_first_arg_fixed sim 1)) ; A000079
(define sim_row_2 (bin2unary_first_arg_fixed sim 2)) ; A011782
(define sim_row_3 (bin2unary_first_arg_fixed sim 3)) ; A068010
(define sim_row_4 (bin2unary_first_arg_fixed sim 4)) ; 1 & A011782, or 1,1 & A000079
(define sim_row_5 (bin2unary_first_arg_fixed sim 5)) ; A068011
(define sim_row_6 (bin2unary_first_arg_fixed sim 6)) ; A068012
(define sim_row_7 (bin2unary_first_arg_fixed sim 7)) ; A068013
(define sim_row_8 (bin2unary_first_arg_fixed sim 8)) ; 1,1,1 & A000079 ?
(define sim_row_9 (bin2unary_first_arg_fixed sim 9)) ; A068030
(define sim_row_10 (bin2unary_first_arg_fixed sim 10)) ; A068031
(define sim_row_11 (bin2unary_first_arg_fixed sim 11)) ; A068032
(define sim_row_12 (bin2unary_first_arg_fixed sim 12)) ; A068033
(define sim_row_13 (bin2unary_first_arg_fixed sim 13)) ; A068034
(define sim_row_14 (bin2unary_first_arg_fixed sim 14)) ; A068035
(define sim_row_15 (bin2unary_first_arg_fixed sim 15)) ; A068036
(define sim_row_16 (bin2unary_first_arg_fixed sim 16)) ; A068037 (interesting...)
(define sim_row_17 (bin2unary_first_arg_fixed sim 17)) ; A068038
(define sim_row_18 (bin2unary_first_arg_fixed sim 18)) ; A068039
(define sim_row_19 (bin2unary_first_arg_fixed sim 19)) ; A068040
(define sim_row_20 (bin2unary_first_arg_fixed sim 20)) ; A068041
(define sim_row_21 (bin2unary_first_arg_fixed sim 21)) ; A068042
(define sim_row_25 (bin2unary_first_arg_fixed sim 25)) ; A068043
(define sim_row_32 (bin2unary_first_arg_fixed sim 32)) ; A068044
(define sim_row_64 (bin2unary_first_arg_fixed sim 64)) ; A068045

(define print_rows
     (lambda (row upto_row upto_nth_term)
          (cond ((> row upto_row) ())
                (else
                   (write-string "The row ") (write row) (write-string ": ")
                   (output_seq (map (bin2unary_first_arg_fixed sim row)
                                    (cons 0 (iota upto_nth_term)))
                   )
                   (newline)
                   (print_rows (1+ row) upto_row upto_nth_term)
                )
          )
     )
)