Home > Archive > Scheme > December 2005 > using collatz to compute euler's constant
You are viewing an archived Text-only version of the thread.
To view this thread in it's original format and/or if you want to reply to
this thread please [click here]
| Author |
using collatz to compute euler's constant
|
|
| ronaldo 2005-11-29, 7:03 pm |
| For those who don't know it, the collatz function is defined
(define (collatz n) ; n > 0
(if (even? n)
(/ n 2)
(+ 1 (* 3 n))))
The mathematician Lothar Collatz conjectured (back in 1937) that
(define (prove n) ; n > 1
(if (=3D n 1)
#t
(prove (collatz n))))
always terminates.
Much research has been done into this prove function, and one
conclusion gets drawn very easily: there is absolutely no relation
between the input number and the number of recursive calls in the
function (i.e. the proof length).
Here is yet another program to investigate the collatz problem.
But, where most research has gone into proof lengths for single
numbers,
I started investigating proof lengths for intervals [1..n].
Running the program shows (this is not a proof but numerical evidence
only) that
prooflength(1..n)
lim(n -> oo) ----------------------- =3D e
n
which is a result that AFAIK has not been found before.
There are three ways to operate the program:
(prove-collatz n 'debug) displays a full trace of the proof for numbers
1 to n
(prove-collatz n 'stats) only displays some state variables at each
proof step
(prove-collatz n 'special) shows the convergence to e
Here is the program (best viewed in 120 col text mode)
;;; collatzr.scm -*-R5RS Scheme-*- Ronald Schr=F6der, 29th of november,
2005
(define (prove-collatz limit mode) (letrec ;;; limit > 0
((collatz (lambda (n) (if (even? n) (/ n 2) (+ 1 (* 3 n)))))
(dump (lambda (it) (begin (display it) (newline))))
(start (lambda (f) (continue f 1 (f 1) '() '() 0 0 0 0 0 0 '())))
(special (lambda (n s) (if (> n 0) (/ (exact->inexact (/ s n)) (exp
1)) 1.0)))
(continue (lambda (f n fn p q r s t u v w y)
(begin
(case mode
((debug) (dump (list (list n fn) (list (+ r s t u v w) r s t u v
w) y (length p) p (length q) q)))
((stats) (dump (list (- n 1) (+ r s t u v w) (list r s t u v w)
(length p) (length q))))
((special) (dump (list (- n 1) (+ r s t u v w) (special (- n 1) (+
r s t u v w))))))
(cond
((> n limit) #t)
((in n p) (continue f (+ n 1) (f (+ n 1)) (rmv n (mrg p q)) '() (+
r 1) s t u v w y))
((in fn p) (continue f (+ n 1) (f (+ n 1)) (mrg p q) '() r
(+ s 1) t u v w y))
((< fn n) (continue f (+ n 1) (f (+ n 1)) (mrg p q) '() r
s (+ t 1) u v w y))
((=3D fn n) (continue f (+ n 1) (f (+ n 1)) (mrg p q) '() r
s t (+ u 1) v w (add n y))) ;;; CYCLE
((> fn n)
(cond
((in fn q) (continue f (+ n 1) (f (+ n 1)) (mrg p q) '() r
s t u (+ v 1) w (add fn y))) ;;; CYCLE
(else (continue f n (f fn) p (add fn q) r
s t u v (+ w 1) y)))))))))
(start collatz)))
(define (in x s)
(if (null? s)
#f
(if (=3D x (car s))
#t
(if (< x (car s))
#f
(in x (cdr s))))))
(define (add x s)
(if (null? s)
(cons x s)
(if (=3D x (car s))
s
(if (< x (car s))
(cons x s)
(cons (car s) (add x (cdr s)))))))
(define (rmv x s)
(if (null? s)
s
(if (=3D x (car s))
(cdr s)
(if (< x (car s))
s
(cons (car s) (rmv x (cdr s)))))))
(define (mrg s1 s2)
(if (null? s1)
s2
(if (null? s2)
s1
(if (=3D (car s1) (car s2))
(mrg (cdr s1) s2)
(if (< (car s1) (car s2))
(cons (car s1) (mrg (cdr s1) s2))
(cons (car s2) (mrg (cdr s2) s1)))))))
| |
| Bradley J Lucier 2005-11-29, 7:03 pm |
| In article <1133274120.307507.305540@f14g2000cwb.googlegroups.com>,
ronaldo <ronald.schroder@hotmail.com> wrote:
<program for verifying Collatz's cojecture>
At first I was thrown off by the subject line, to me Euler's constant is
$$
\gamma=\lim_{n\to\infty}(\log n - \sum_1^n {1\over n});
$$
see the very nice book "Gamma : Exploring Euler's Constant" by
Julian Havil.
Can you say in a few words what your proof is doing? Perhaps I can
ask some questions in between comments.
I see that it's counting the number of times each clause in the two
nested cond's is called.
p is the set of integers that have been proved to either finish or have
a cycle (of course, since there are no cycles up to 3*2^{53}, we won't
find a cycle here in a reasonable amount of computer time). So if you want
to just count the number of steps for the "proof" you could eliminate the
cycle checks, right?
q is the set of integers we visit during the computation for the current n.
You remove n from the set p, but then add an extra check for (< fn n),
so it would be the same if you *didn't* remove n from p. I suppose checking
(< fn n) is faster than checking for n in p.
It appears that the only time n will be in p will be when q is empty, so
perhaps you don't need to call the mrg? Also, you call (in n p) each step
in verifying the conjecture for n, and since p doesn't change in the iteration
for fixed p, it will be false each time it's called (if it's false the first
time) so this seems rather artificial.
I did a quick revision to just return the total number of steps in the proof,
got rid of the reporting, and got:
> (time (prove-collatz 10000 #t))
(time (prove-collatz 10000 #t))
3286 ms real time
2902 ms cpu time (2733 user, 169 system)
5 collections accounting for 17 ms real time (16 user, 1 system)
523294256 bytes allocated
no minor faults
no major faults
27338
> (time (prove-collatz 100000 #t))
(time (prove-collatz 100000 #t))
480681 ms real time
384399 ms cpu time (380801 user, 3598 system)
1109 collections accounting for 12997 ms real time (10617 user, 116 system)
108690291088 bytes allocated
no minor faults
no major faults
273946
(This was compiled with no declarations in Gambit-C 4.0b.)
Do you have an estimate of the big-O complexity of your function? It appears
to be O(N^2).
Thanks for the interesting function/conjecture. There's a lot of stuff
on Collatz's conjecture on the web and I didn't find your conjecture after
a quick search.
Brad
| |
| Bradley J Lucier 2005-11-29, 7:03 pm |
| In article <dmiikb$d1r@arthur.cs.purdue.edu>,
Bradley J Lucier <bjl@cs.purdue.edu> wrote:
>In article <1133274120.307507.305540@f14g2000cwb.googlegroups.com>,
>ronaldo <ronald.schroder@hotmail.com> wrote:
><program for verifying Collatz's cojecture>
....
>It appears that the only time n will be in p will be when q is empty, so
>perhaps you don't need to call the mrg? Also, you call (in n p) each step
>in verifying the conjecture for n, and since p doesn't change in the iteration
>for fixed p, it will be false each time it's called (if it's false the first
>time) so this seems rather artificial.
Sorry, I see that you increment one variable each time through the loop,
so it doesn't matter that there are some extra tests in the loop.
Brad
| |
| ronaldo 2005-11-29, 7:03 pm |
| >Can you say in a few words what your proof is doing? Perhaps I can
>ask some questions in between comments.
>I see that it's counting the number of times each clause in the two
>nested cond's is called.
Basically, I count the number of times that a recursive call to collatz
is being done (which is exactly one time in each call to continue) -
and I count the whys of it.
>You remove n from the set p, but then add an extra check for (< fn n),
>so it would be the same if you *didn't* remove n from p. I suppose checking
>(< fn n) is faster than checking for n in p.
The loop invariant is that p is the set of numbers which are known to
reduce to 1, but of course keeping numbers smaller than n in a set is
far less efficient than a simple test. So yes, this is an optimization
hack. Nevertheless, the program slows down rapidly as n increases.My
guess is O(n^2.718281828459045) ;-).
But come to speak of the big O, you are referring to time complexity.
My algoritm was designed to explore another kind of complexity so
execution timings were not my main interest.
| |
| ronaldo 2005-11-29, 7:03 pm |
| >Thanks for the interesting function/conjecture.
You're welcome.
| |
| Bradley J Lucier 2005-12-02, 9:24 pm |
| In article <1133274120.307507.305540@f14g2000cwb.googlegroups.com>,
ronaldo <ronald.schroder@hotmail.com> wrote:
>Running the program shows (this is not a proof but numerical evidence
>only) that
>
> prooflength(1..n)
> lim(n -> oo) ----------------------- = e
> n
>
>which is a result that AFAIK has not been found before.
OK, I wrote a small bitset library using some of Gambit-C's internal
bignum routines and used exact integers as bitsets. Unfortunately,
the bitsets are *very* sparse, so while this representation has O(1)
access time rather than O(N) access time for a set of $N$ elements,
I rapidly ran out of space. (For example, it's impossible in 8GB
of memory to run prove-collatz using exact-integers as bitsets.)
So I turned to the built-in "tables" in Gambit-C 4.0b to implement bitsets.
Tables turned out to have a good balance between access time and space
requirements.
So here's the result. If the limit exists, I don't think the limit is e :-(.
Each run reports, in order, the number of steps of the "proof", the maximum
value computed by the collatz function in the run, and the number of nonzero
elements in the bitsets.
> (time (prove-collatz 10000 #t))
(time (prove-collatz 10000 #t))
85 ms real time
84 ms cpu time (62 user, 22 system)
1 collection accounting for 26 ms real time (6 user, 19 system)
1451696 bytes allocated
no minor faults
no major faults
27337
27114424
21664
> (time (prove-collatz 100000 #t))
(time (prove-collatz 100000 #t))
979 ms real time
877 ms cpu time (843 user, 34 system)
5 collections accounting for 48 ms real time (41 user, 7 system)
16409160 bytes allocated
no minor faults
no major faults
273945
1570824736
217212
> (time (prove-collatz 1000000 #t))
(time (prove-collatz 1000000 #t))
13869 ms real time
11353 ms cpu time (10956 user, 397 system)
10 collections accounting for 526 ms real time (379 user, 65 system)
185659296 bytes allocated
no minor faults
no major faults
2736247
56991483520
2168611
> (time (prove-collatz 10000000 #t))
(time (prove-collatz 10000000 #t))
464667 ms real time
390909 ms cpu time (385759 user, 5150 system)
11 collections accounting for 3914 ms real time (2814 user, 528 system)
1485295920 bytes allocated
no minor faults
no major faults
27407525
60342610919632
21730849
Here's the program:
(declare (standard-bindings)(extended-bindings)(block)(not safe))
(define-structure bitset body)
#|
;;; Bitsets as exact integers
(define (empty-bitset)
(make-bitset 0))
;;; note that all the following functions work even if (bitset-body s) is an
;;; unnormalized bignum.
(define (bitset-bit-set? n s)
(bit-set? n (bitset-body s)))
(define (bitset-size s)
(bit-count (bitset-body s)))
(define (bitset-maximum-element s)
(let ((body (bitset-body s)))
(and (not (eq? body 0))
(- (integer-length body) 1))))
(define (bitset-bit-set! n s)
(define (set-bit! n body)
(let ((index (##fixnum.quotient n ##bignum.mdigit-width))
(shift (##fixnum.remainder n ##bignum.mdigit-width)))
(##bignum.mdigit-set! body index (##fixnum.bitwise-ior (##bignum.mdigit-ref body index)
(##fixnum.arithmetic-shift-left 1 shift)))))
(let ((body (bitset-body s)))
(cond ((##fixnum? body)
(if (and (##fixnum? n)
(##fixnum.< n (##fixnum.- ##fixnum-width 1)))
(bitset-body-set! s (##fixnum.bitwise-ior body (##fixnum.arithmetic-shift-left 1 n)))
(let ((body (##bignum.make (##quotient (+ n ##bignum.adigit-width 1) ##bignum.adigit-width) (##bignum.<-fixnum body) #f)))
(bitset-body-set! s body)
(set-bit! n body))))
((< n (##fixnum.- (##fixnum.* (##bignum.mdigit-length body)
##bignum.mdigit-width)
1))
(set-bit! n body))
(else
;; we increase the capacity of the bitset by at least 50%
;; note that this means that the body is an un-normalized bignum, but that's OK for our purposes.
(let ((body (##bignum.make (##max (##quotient (+ ##bignum.adigit-width n 1) ##bignum.adigit-width)
(##quotient (* (##bignum.adigit-length body) 3) 2))
body
#f)))
(bitset-body-set! s body)
(set-bit! n body))))
s))
|#
;;; bitsets as tables
(define (empty-bitset)
(make-bitset (make-table)))
(define (bitset-bit-set? n s)
(table-ref (bitset-body s) n #f))
(define (bitset-size s)
(table-length (bitset-body s)))
(define (bitset-maximum-element s)
(let ((max-element -1))
(table-for-each (lambda (key value)
(set! max-element (max max-element value)))
(bitset-body s))
max-element))
(define (bitset-bit-set! n s)
(table-set! (bitset-body s) n n)
s)
(define (prove-collatz limit mode)
;; mode is ignored
(letrec ;;; limit > 0
((collatz (lambda (n) (if (even? n) (quotient n 2) (+ 1 (* 3 n)))))
(start (lambda (f) (continue f 1 (f 1) (bitset-bit-set! 1 (empty-bitset)) 0 0 0)))
(continue (lambda (f n fn p r s t)
(cond
((##fixnum.> n limit)
(values (+ r s t)
(bitset-maximum-element p)
(bitset-size p)))
((bitset-bit-set? n p)
(continue f (+ n 1) (f (+ n 1)) p (+ r 1) s t))
(else
(letrec ((loop (lambda (f n fn p r s t)
(if (bitset-bit-set? fn p)
(continue f (+ n 1) (f (+ n 1)) (bitset-bit-set! n p) r (+ s 1) t)
(loop f n (f fn) (bitset-bit-set! fn p) r s (+ t 1))))))
(loop f n fn p r s t)))))))
(start collatz)))
| |
| ronaldo 2005-12-02, 9:24 pm |
| Gee, thanks, Brad! These are numbers I could never reproduce myself.
I'm running the rhizome/pi compiler on a 800MhZ machine with only 256
MB of memory. And rhizome/pi I doesn't support bitsets.
Meanwhile, I found some more results.
First I removed some reduntant tests, as you suggested. Then, instead
of the collatz function I used a kind of generalized collatz function
where you add an extra division (by the old multiplier) and replace the
multiplier with the next prime number. The next function then is:
(define (c5 n)
(if (even? n)
(/ n 2)
(if (divisible n 3)
(/ n 3)
(+ 1 (* 5 n)))))
I named these functions c3, c5, c7, etc. to name the divisor (now c3 =
collatz). I let run upto c31 for n = 10000 and upto c7 for n = 100000.
I found cycles for c11: n =17, c13: n = 19, c17: n = 43, c23: n = 179
and c31: n = 67. A premature conclusion might be that cycles occur for
small primes only.
Results for n = 10000:
f calls to f |p|
c3 27339 21664
c5 22223 18273
c7 22532 19555
c11 22749 20364
c13 25794 23696
c17 22845 20974
c19 25599 23883
c23 23724 22253
c29 24303 22977
c31 26398 25125
where |p| is the size of the set of all numbers proven to reduce either
to 1 or to another attractor.
These results prompt me to stick to my conjecture - for the time being.
Showing it to be true by means of computation seems to be impossible
though.
| |
| ronaldo 2005-12-03, 7:00 pm |
| ;;; collatzr.scm -*-R5RS Scheme-*- Ronald Schr=F6der, 3rd of december,
2005
;;;
;;; EUREKA ?
;;;
;;; The program as it was just doesn't show it as well for
; some numbers as for others. You miss out on this if you turn
; off the display.
;;; I've been asking the wrong question. Instead of asking how
; many proof steps are required for n I should ask upto which n
; I can progress given a number of proof steps.
; In other words, fuel it and see how far it gets.
; Let's feed it the decimal expansion of e
; collatzr: (prove c3 1000000 272 #t)
; (79 272 221)
; collatzr: (prove c3 1000000 2718 #t)
; (967 2718 2171)
; collatzr: (prove c3 1000000 27183 #t)
; (9934 27183 21548)
; collatzr: (prove c3 1000000 271828 #t)
; (99199 271828 215543)
; Does it converge? Unfortunately, the next step is beyond
; my computer's capacity. Would you be so kind?
; The program -below- has been transformed to be
; "copy-and-paste"-able. Run it through pp.
(define (prove collatz max-n max-z mode)
(letrec
((dump (lambda (it)
(begin (display it) (newline))))
(start (lambda (f)
(continue f 0 (f 0) '() '() 0 0 0 0 0 0 #f '())))
(continue (lambda (f n fn p q r s t u v w nn y)
(let ((z (+ r s t u v w)))
(begin
(case mode
((prove)
(dump
(list
(list n fn)
(list z r s t u v w) y
(length p) p (length q) q)))
((stats)
(if nn (dump (list (- n 1) z (list r s t u v w))) #t))
(else #t))
(cond
((> n max-n)
(list (- n 1) z (+ (- n 1) (length p))))
((> z max-z)
(list n (- z 1) (+ n (length p))))
((and nn (in n p))
(continue
f (+ n 1) (f (+ n 1))
(rmv n (mrg p q)) '() (+ r 1) s t u v w #t y))
((=3D fn n)
(continue
f (+ n 1) (f (+ n 1))
(mrg p q) '() r (+ s 1) t u v w #t (add n y)))
((and (>=3D fn 0) (< fn n))
(continue
f (+ n 1) (f (+ n 1))
(mrg p q) '() r s (+ t 1) u v w #t y))
((in fn q)
(continue
f (+ n 1) (f (+ n 1))
(mrg p q) '() r s t (+ u 1) v w #t (add fn y)))
((in fn p)
(continue
f (+ n 1) (f (+ n 1))
(mrg p q) '() r s t u (+ v 1) w #t y))
(else
(continue
f n (f fn)
p (add fn q) r s t u v (+ w 1) #f y))))))))
(start collatz)))
(define (? p q) (=3D (modulo p q) 0))
(define (c3 n) (if (even? n) (/ n 2) (+ 1 (* n 3))))
(define (c5 n) (if (even? n) (/ n 2) (if (? n 3) (/ n 3)
(+ 1 (* n 5)))))
(define (c7 n) (if (even? n) (/ n 2) (if (? n 3) (/ n 3)
(if (? n 5) (/ n 5) (+ 1 (* n 7))))))
(define (c11 n) (if (even? n) (/ n 2) (if (? n 3) (/ n 3)
(if (? n 5) (/ n 5) (if (? n 7) (/ n 7) (+ 1 (* n 11)))))))
(define (in x s)
(if (null? s) #f
(if (=3D x (car s)) #t
(if (< x (car s)) #f
(in x (cdr s))))))
(define (add x s)
(if (null? s) (cons x s)
(if (=3D x (car s)) s
(if (< x (car s)) (cons x s)
(cons (car s) (add x (cdr s)))))))
(define (rmv x s)
(if (null? s) s
(if (=3D x (car s)) (cdr s)
(if (< x (car s)) s
(cons (car s) (rmv x (cdr s)))))))
(define (mrg s1 s2)
(if (null? s1) s2
(if (null? s2) s1
(if (=3D (car s1) (car s2)) (mrg (cdr s1) s2)
(if (< (car s1) (car s2)) (cons (car s1) (mrg s2 (cdr s1)))
(cons (car s2) (mrg (cdr s2) s1)))))))
| |
| ronaldo 2005-12-04, 7:56 am |
| UPDATE
collatzr: (prove d3 1000000 2718 #t)
(996 2718 2159)
collatzr: (prove d3 1000000 27183 #t)
(9981 27183 21518)
collatzr: (prove d3 1000000 271828 'stats)
(98957 271828 215654)
The function d3 is almost the same as c3, replacing the
addition with a subtraction:
(define (d3 n)
(if (even? n)
(/ n 2)
(- (* n 3) 1)))
This function seems to approach the same limit as c3.
By the way, d3 contains 4 cycles in the proven range
(n one of (0 1 5 17)).
| |
| ronaldo 2005-12-05, 3:57 am |
| MORE NEWS
I have generated graphics for n -> prooflength(1..n)/n for c3 and d3
for n upto 10000..
As usual, they're worth more than a thousand words. Since it's
impossible to include
them here I will only say that they confirm my suspicions. These are
irregular sawtoothed
functions as one would expect - but they do seem to converge on e.
My scheme has no plotting library, I used the free Euler mathematics
program instead. It took
me less than 20 minutes to install it and learn some basic operations
like reading data from
a file and generating a function graph.
Links:
http://mathsrv.ku-eichstaett.de/MGF...grothman/euler/
| |
| ronaldo 2005-12-05, 7:05 pm |
| OK. My graphics are now up to n = 100000.
I combined c3 and d3 in one picture and it's
beautiful. I'd like to show it to everybody
on the Internet but I don't have my own website.
Any suggestions?
| |
| Anton van Straaten 2005-12-05, 9:58 pm |
| ronaldo wrote:
> OK. My graphics are now up to n = 100000.
> I combined c3 and d3 in one picture and it's
> beautiful. I'd like to show it to everybody
> on the Internet but I don't have my own website.
> Any suggestions?
It's probably a bit off-topic for the Scheme Cookbook, and the Community
Scheme Wiki doesn't currently support image uploads, but you could post
it at base.google.com -- just select "Post an Item", and you can
describe it and attach GIF, JPEG, PNG or TIFF images.
Anton
| |
|
|
| ronaldo 2005-12-06, 7:09 pm |
| I'm now gathering statistics about the
individual program branches. More
surprises: the following results
appear to be independent of n!
>(bitset-bit-set? n p)
This happens about 20 iterations out of 100.
> (if (bitset-bit-set? fn p)
This happens about 15 iterations out of 100.
Approximately 4 times out of these 15 fn < n.
Room for optimization!
>(loop f n (f fn) ...
This happens about 63 iterations out of 100.
| |
|
|
|
|
|