[plt-scheme] Numerical precision (sqrt 2)

From: Andreas Rottmann (a.rottmann at gmx.at)
Date: Sat Jun 6 19:23:38 EDT 2009

rw <robert.winslow at gmail.com> writes:

> I am seeking a (hopefully built-in) way to decide the precision of
> irrational numbers. By default (PLT 4.1.5) gives:
>
> (sqrt 2) -> 1.4142135623730951
>
> I want to get an answer in (say) 50 or 100 digits precision. I am
> considering building a generating function for each digit of the
> irrational number, but I expect this to be computationally slow.
> Ideally, it would operate like this:
>
> (precision-take (sqrt 2) 50) -> 2^0.5 with 50 significant digits
>
> Is functionality like this already out there?
>
I've done something similiar for the "golden number". The code works by
using an approximation function that works with rationals of growing
precision, periodically checking if the approximation is good
enough. You could probably transform the algorithm produce a digit
stream to get an interface like you described. Perhaps you find the code
useful (hereby placed in the public domain):

#lang scheme

(define (prev-power num base)
  (let ((l (/ (log num) (log base))))
    (max 1 (inexact->exact (expt base (inexact->exact (floor l)))))))

(define *decimal-digits* '#(#\0 #\1 #\2 #\3 #\4 #\5 #\6 #\7 #\8 #\9))

(define (number->string/prec num base prec)
  (let loop ((digits '())
             (n num)
             (prec prec)
             (place-val (prev-power num base)))
    (if (= prec 0)
        (list->string (reverse digits))
        (let* ((digit (inexact->exact (truncate (/ n place-val))))
               (digit-char (vector-ref *decimal-digits* digit)))
          (loop (if (= place-val 1)
                    (cons #\. (cons digit-char digits))
                    (cons digit-char digits))
                (- n (* digit place-val))
                (- prec 1)
                (/ place-val base))))))

(define (golden-ratio epsilon)
  (let loop ((f1 1) (f2 1) (i 0))
    (let ((f3 (+ f1 f2)))
      (if (= i 1000)
          (let ((r1 (/ f2 f1))
                (r2 (/ f3 f2)))
            (if (< (abs (- r2 r1)) epsilon)
                r2
                (loop f2 f3 0)))
          (loop f2 f3 (+ i 1))))))

(define (golden n)
  (display (number->string/prec (golden-ratio (expt 1/10 n)) 10 n))
  (newline))

(golden 100)


Posted on the users mailing list.