First, the code:
;;; Simple Continued Fraction
;;; [a0, a1, a2 ... an]
(defmacro while (test &body body)
`(do ()
((not ,test))
,@body))
(defun make-stack ()
(make-array 1
:element-type 'fixnum
:adjustable t
:fill-pointer 0))
(defun push-stack (x stack)
(vector-push-extend x stack))
(defun pop-stack (stack)
(vector-pop stack))
(defun empty-stack-p (stack)
(= (fill-pointer stack) 0))
(defun convergent-float (generator convergent)
(when (/= convergent 0)
(let ((stack (make-stack))
(val 1.0s0))
(dotimes (i convergent)
(push-stack (funcall generator) stack))
(while (not (empty-stack-p stack))
(setf val (/ 1.0s0 val))
(incf val (pop-stack stack)))
val)))
;;; generates Bill Gosper's e expansion as
;;; [1 0 1 1 2 1 1 4 1 1 6 1 1 8 1 1 10 1 1 12 ... ]
(defun e-generator ()
(let ((index 0))
(lambda ()
(let ((val 0))
(cond ((eq index 0) (setf val 1))
((eq index 1) (setf val 0))
(t (progn
(setf val (- index 2))
(if (eq (rem val 3) 2)
(setf val (* 2 (+ 1 (floor (/ val 3)))))
(setf val 1)))))
(incf index)
val))))
CL-USER 1 > (convergent-float (e-generator) 20)
2.718277S0
Now my problems:
The code works, well, sort of. In fact it doesn't, because
what I want is to get to an as close as possible approximation
of the convergent, and no fixed floating point format whether single
or double is adequate.
The question is: what are my options? Ideally, I would want some
kind of "decimal bignum", ie a kind of arbitrary large float
format.
My second problem is with the inelegance/inefficiency of my solution.
I couldn't come up with an algorithm that would iterate from the
start of the sequence, so, internally, I use a stack such that the
generator can still produce its coefficients in order, yet I can
compute the result from the end.
But that's very undesirable, because under some circumstances,
I will want the generator to produce an arbitrarily large number
of coefficients, and yet only decide at runtime to stop calling
it after a given precision has been reached, which I can't know
until after the stack has completely absorbed all the generator
output (and run out of memory in the process ...)
Or any math package that I missed on Cliki that handles
convergents?
Many Thanks
--
JFB
On May 2, 7:31 am, verec wrote:
>
> The code works, well, sort of. In fact it doesn't, because
> what I want is to get to an as close as possible approximation
> of the convergent, and no fixed floating point format whether single
> or double is adequate.
>
> The question is: what are my options? Ideally, I would want some
> kind of "decimal bignum", ie a kind of arbitrary large float
> format.
I wrote a small continued fraction program in Scheme. It would
be fairly easy to port it to Lisp. See
http://jrm-code-project.googlecode.com/svn/trunk/scheme/contfrac.scm
(cf:render e 100)
2.718281828459452353628747135266249775724793699959574966967627724766303535475945713821785251664274274663919320359921817413...
Hi Joe,
> I wrote a small continued fraction program in Scheme. It would
> be fairly easy to port it to Lisp. See
> http://jrm-code-project.googlecode.com/svn/trunk/scheme/contfrac.scm
>
> (cf:render e 100)
> 2.718281828459452353628747135266249775724793699959574966967627724766303535475945713821785251664274274663919320359921817413...
>
I love this comment:
;;; ``... the problem of finding the continued
;;; fraction for a sum from the continued fractions
;;; representing the addends is exceedingly complicated,
;;; and unworkable in computational practice.''
;;;
;;; --A. YA. KHinchin, 1935
;;; Well, it is *more* complicated than simple addition, but it really
;;; isn't that hard at all.
;;;
;;; This implements simple arithmetic on continued fractions as described
;;; by Gosper in
;;;
;;; http://www.inwap.com/pdp10/hbaker/hakmem/cf.html
;;; and
;;;
;;; http://www.tweedledum.com/rwg/cfup.htm
Is it from "Continued Fractions" by Khinchin?
(A delightful book - I really enjoy his style of
exposition).
--
Jens Axel S�gaard
On May 2, 11:29 am, Jens Axel Søgaard <······@soegaard.net> wrote:
>
> I love this comment:
>
> ;;; ``... the problem of finding the continued
> ;;; fraction for a sum from the continued fractions
> ;;; representing the addends is exceedingly complicated,
> ;;; and unworkable in computational practice.''
> ;;;
> ;;; --A. YA. KHinchin, 1935
>
> Is it from "Continued Fractions" by Khinchin?
> (A delightful book - I really enjoy his style of
> exposition).
That's the one. ISBN 0486696308
On May 2, 2:06 pm, Joe Marshall <··········@gmail.com> wrote:
> I wrote a small continued fraction program in Scheme. It would
> be fairly easy to port it to Lisp. Seehttp://jrm-code-project.googlecode.com/svn/trunk/scheme/contfrac.scm
>
> (cf:render e 100)
> 2.7182818284594523536287471352662497757247936999595749669676277247663035354 75945713821785251664274274663919320359921817413...
Joe:
The comment at the beginning of the file:
;;; 7. Try computing this function (from Muller):
;;;
;;; a0 = 61/11
;;;
;;; a1 = 11/2
;;;
;;; 3000
;;; 1130 - --------
;;; an-2
;;; an = 111 - ------------------
;;; an-1
;;;
;;; (imagine that the n-1 and n-2 are subscripts)
;;;
;;; The answer converges to 6, but if you use floating point, it
will
;;; converge to 100.
is not quite right, it should be
;;; a1 = 61/11
;;;
;;; a0 = 11/2
(see the actual code at the end of the file).
I decided to test this with my computable reals package, and with the
corrected initial conditions got
> (computable->inexact (a_ 30))
5.995804952329115
which matches what your rational computation gives for a_29 (so we
disagree on the numbering). I also find
> (computable->inexact (a_ 100))
5.999999987925326
> (computable->inexact (a_ 200))
6.
Fun.
Brad
On May 2, 9:28 pm, ······@math.purdue.edu wrote:
>
> Joe:
>
> The comment at the beginning of the file:
> is not quite right, it should be
>
> ;;; a1 = 61/11
> ;;;
> ;;; a0 = 11/2
Fixed.
> > (computable->inexact (a_ 30))
>
> 5.995804952329115
>
> which matches what your rational computation gives for a_29 (so we
> disagree on the numbering).
My mistake, n was off by one. I fixed this as well.
> I also find
>
> > (computable->inexact (a_ 100))
> 5.999999987925326
> > (computable->inexact (a_ 200))
>
> 6.
I get these values as well.
Thanks for the notes!
On Wed, 2 May 2007 15:31:33 +0100, verec wrote:
>Or any math package that I missed on Cliki that handles
>convergents?
>
There is this library:
http://www.haible.de/bruno/MichaelStoll/reals.html
which probably does what you want.
--
|Don't believe this - you're not worthless ,gr---------.ru
|It's us against millions and we can't take them all... | ue il |
|But we can take them on! | @ma |
| (A Wilhelm Scream - The Rip) |______________|
From: ···············@space.at
Subject: Re: Simple Continued Fraction
Date:
Message-ID: <u647bustr.fsf@space.at>
>>>>> "verec" == verec <verec> writes:
verec> First, the code:
verec> ;;; Simple Continued Fraction
verec> ;;; [a0, a1, a2 ... an]
[snip]
verec> (defun convergent-float (generator convergent)
verec> (when (/= convergent 0)
verec> (let ((stack (make-stack))
verec> (val 1.0s0))
^^^^^
verec> (dotimes (i convergent)
verec> (push-stack (funcall generator) stack))
verec> (while (not (empty-stack-p stack))
verec> (setf val (/ 1.0s0 val))
^^^^^
verec> (incf val (pop-stack stack)))
verec> val)))
[snip]
verec> Now my problems:
verec> The code works, well, sort of. In fact it doesn't, because
verec> what I want is to get to an as close as possible
verec> approximation of the convergent, and no fixed floating
verec> point format whether single or double is adequate.
verec> The question is: what are my options? Ideally, I would want
verec> some kind of "decimal bignum", ie a kind of arbitrary large
verec> float format.
You could use rational numbers, it probably would be sufficient to
replace the 1.0s0 marked in your code by 1, but please check.
If you use CLISP, you can use the long float which has adjustable
precision. I believe that Maxima has something similar (bigfloats),
but I have no experience with it.
verec> My second problem is with the inelegance/inefficiency of my
verec> solution.
verec> I couldn't come up with an algorithm that would iterate
verec> from the start of the sequence, so, internally, I use a
verec> stack such that the generator can still produce its
verec> coefficients in order, yet I can compute the result from
verec> the end.
verec> But that's very undesirable, because under some
verec> circumstances, I will want the generator to produce an
verec> arbitrarily large number of coefficients, and yet only
verec> decide at runtime to stop calling it after a given
verec> precision has been reached, which I can't know until after
verec> the stack has completely absorbed all the generator output
verec> (and run out of memory in the process ...)
See http://mathworld.wolfram.com/ContinuedFraction.html, equations
(16) and (17).
Since you mentioned Bill Gosper, have you read his entries in HAKMEM?
http://home.pipeline.com/~hbaker1/hakmem/cf.html
verec> Or any math package that I missed on Cliki that handles
verec> convergents?
Hope this helps,
Roland
+ ···············@space.at:
| See http://mathworld.wolfram.com/ContinuedFraction.html, equations
| (16) and (17).
Indeed. Before I saw Roland's response, I played around a bit and
hacked up this code:
(defpackage :cfrac (:use :cl))
(in-package :cfrac)
(defun next-convergent (a b c d)
(lambda (next-coefficient)
(values (+ (* a next-coefficient) b)
a
(+ (* c next-coefficient) d)
c)))
(defun show-first-n-convergents (generator n)
(flet ((gen () (funcall generator)))
(let ((a (gen)) (b 1) (c 1) (d 0))
(dotimes (k n (/ a c))
(multiple-value-setq (a b c d)
(funcall (multiple-value-call #'next-convergent a b c d) (gen)))
(format t "~&(~d+~dq)/(~d+~dq) ~a=~f" a b c d
(/ a c)
(float (/ a c) 1.0d0))))))
(defun e-generator ()
"Generate the sequence 2 1 2 1 1 4 1 1 6 1 1 8 1 1 10 ..."
(let ((c nil) (n 0))
(lambda ()
(cond ((null c)
(setf c '#1=(nil nil t . #1#)) ; couldn't resist ...
2)
((null (car (setf c (cdr c))))
1)
(t (incf n 2))))))
And here is a trial run:
cfrac> (show-first-n-convergents (e-generator) 20)
(3+2q)/(1+1q) 3=3.0
(8+3q)/(3+1q) 8/3=2.6666666666666665
(11+8q)/(4+3q) 11/4=2.75
(19+11q)/(7+4q) 19/7=2.7142857142857144
(87+19q)/(32+7q) 87/32=2.71875
(106+87q)/(39+32q) 106/39=2.717948717948718
(193+106q)/(71+39q) 193/71=2.7183098591549295
(1264+193q)/(465+71q) 1264/465=2.718279569892473
(1457+1264q)/(536+465q) 1457/536=2.718283582089552
(2721+1457q)/(1001+536q) 2721/1001=2.7182817182817183
(23225+2721q)/(8544+1001q) 23225/8544=2.7182818352059925
(25946+23225q)/(9545+8544q) 25946/9545=2.7182818229439496
(49171+25946q)/(18089+9545q) 49171/18089=2.718281828735696
(517656+49171q)/(190435+18089q) 517656/190435=2.7182818284454013
(566827+517656q)/(208524+190435q) 566827/208524=2.718281828470584
(1084483+566827q)/(398959+208524q) 1084483/398959=2.7182818284585633
(13580623+1084483q)/(4996032+398959q) 13580623/4996032=2.718281828459065
(14665106+13580623q)/(5394991+4996032q) 14665106/5394991=2.718281828459028
(28245729+14665106q)/(10391023+5394991q) 28245729/10391023=2.718281828459046
(410105312+28245729q)/(150869313+10391023q) 410105312/150869313=2.718281828459045
410105312/150869313
--
* Harald Hanche-Olsen <URL:http://www.math.ntnu.no/~hanche/>
- It is undesirable to believe a proposition
when there is no ground whatsoever for supposing it is true.
-- Bertrand Russell
+ Harald Hanche-Olsen <······@math.ntnu.no>:
| (funcall (multiple-value-call #'next-convergent a b c d) (gen)))
What a brain fart. That is of course just a super complicated way of
saying (funcall (next-convergent a b c d) (gen)). Duh!
I must have had multiple values on my brain.
--
* Harald Hanche-Olsen <URL:http://www.math.ntnu.no/~hanche/>
- It is undesirable to believe a proposition
when there is no ground whatsoever for supposing it is true.
-- Bertrand Russell
verec writes:
> First, the code:
> (defun convergent-float (generator convergent)
> (when (/= convergent 0)
> (let ((stack (make-stack))
> (val 1.0s0))
> (dotimes (i convergent)
> (push-stack (funcall generator) stack))
> (while (not (empty-stack-p stack))
> (setf val (/ 1.0s0 val))
> (incf val (pop-stack stack)))
> val)))
Change convergent-float to
(defun convergent-float (generator convergent)
(when (/= convergent 0)
(let ((stack (make-stack))
(val 1))
(dotimes (i convergent)
(push-stack (funcall generator) stack))
(while (not (empty-stack-p stack))
(setf val (/ 1 val))
(incf val (pop-stack stack)))
val)))
> Now my problems:
>
> The code works, well, sort of. In fact it doesn't, because
> what I want is to get to an as close as possible approximation
> of the convergent, and no fixed floating point format whether single
> or double is adequate.
>
> The question is: what are my options? Ideally, I would want some
> kind of "decimal bignum", ie a kind of arbitrary large float
> format.
>
Then you have arbitrarily large decimal bignums.
CL-USER> (convergent-float (e-generator) 20)
14665106/5394991
CL-USER> (float *)
2.7182817
CL-USER> (float ** 1d0)
2.718281828459028D0
CL-USER>
Wade
On 2007-05-02 15:31:33 +0100, verec said:
> First, the code:
[...]
> Many Thanks
Thank you all for all the pointers and examples!
I see I've got much homework to do! :-)
Thanks again
--
JFB
verec writes:
> (defun convergent-float (generator convergent)
> (when (/= convergent 0)
> (let ((stack (make-stack))
> (val 1.0s0))
> (dotimes (i convergent)
> (push-stack (funcall generator) stack))
> (while (not (empty-stack-p stack))
> (setf val (/ 1.0s0 val))
> (incf val (pop-stack stack)))
> val)))
It seems to me you may possibly have used the short-float syntax by
mistake; if not, I apologize. Are you aware that 1.0s0 is a
short-float, not a single-float, and that short-floats may have very
little precision?
http://clhs.lisp.se/Body/02_cbb.htm
http://clhs.lisp.se/Body/t_short_.htm