prog

frontpage - thread list - new thread - preferences - ?

e and the Stern-Brocot tree

3 2021-02-09 06:56

Irrational numbers don't appear in the Stern-Brocot tree, but all the rational numbers that are ``close'' to them do. We can find the infinite representation of an irrational number α with this procedure:

if α < 1 then (output(L); α := α / (1 - α))
         else (output(R); α := α - 1)

In Scheme

(define (brocot-string num steps)
  (letrec ((helper (lambda (num acc steps)
                     (cond ((equal? steps 0) (reverse acc))
                           ((< num 1) (helper (/ num (- 1 num)) (cons 'L acc) (- steps 1)))
                           (else (helper (- num 1) (cons 'R acc) (- steps 1)))))))
    (helper num '() steps)))

Now we can look for a regular pattern

> (brocot-string 2.7182818284590452353599104282424311632484 32)
'(R R L R R L R L L L L R L R R R R R R L R L L L L L L L L R L R)

In the Stern-Brocot system e = R²LR²LRL⁴RLR⁶LRL⁸RLR¹⁰LRL¹²....

We'll use infinite streams to represent this string. We start with the sequence of bases RLRLRL...

(require srfi/41)
(define bases (stream-cons R (stream-cons L bases)))

The sequence of exponents is (2,1,1,2,1,1,4,1,1,6,1,1,8,1,1,10...) You may recognize Euler's continued fraction https://oeis.org/A003417

(define exponents
  (stream-map
   (lambda (x)
     (cond ((equal? x 1) 2)
           ((equal? (modulo x 3) 0) (- x (/ x 3)))
           (else 1)))
   (stream-from 1)))

Now we just have to zip the two streams

(define exp-pairs (stream-zip bases exponents))

and we build the infinite sequence of rational approximations of e

(define-stream (exponentiate a b)
  (if (equal? b 0)
      empty-stream
      (stream-cons a (exponentiate a (- b 1)))))

(define-stream (make-e-seq pairs)
  (let ((a (cadr (stream-car pairs)))
        (b (car (stream-car pairs))))
    (if (equal? a 1)
      (stream-cons b (make-e-seq (stream-cdr pairs)))
      (stream-append (exponentiate b a) (make-e-seq (stream-cdr pairs))))))
      
(define e-seq (make-e-seq exp-pairs))

We define a function to calculate the nth rational approximation of e in our Stern-Brocot representation like this

(define (e ref)
  (mat->rat (stream-ref (stream-scan mat-mul I e-seq) ref)))

now let's test our code

> (require math/bigfloat)
> (bf-precision 100000)
> (time (bf (eval (e 1000000))))
cpu time: 3258 real time: 3260 gc time: 1711
(time (bf (eval (e 1000000))))
(bf #e2.7182818284590452353...)

Yields a result correct to 6331 digits. While far from being the most efficient way to approximate e, that was fun!

54