; ; 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://oeis.org/A066425 [Update URL, Mar 06 2015] ; ; 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)