The following program to compute the squareroot of a number x (I teach
as some of you know) works of course not very well:
(setf (ext:long-float-digits) 20000)
(setq epsilon 0.0000000000000000000000001L0)
(defun heron-sqrt (x)
(let ((tmp 1.0L0))
(loop
(setq tmp (* 0.5L0 (+ tmp (/ x tmp))))
(if (< (abs (- (* tmp tmp) x)) epsilon)
(return tmp)))))
(print (heron-sqrt 2.0L0))
The problem is that my epsilon has been introduced artificially. How do I
perform this computation (after setting long-float-digits) as precisely
aas possible? I do not think I can use long-float-epsilon here?
Additionally, I do not understand either, where I have to put the "L0" at
the end of a nuber and when this is not important.
TIA,
jb
From: Steven E. Harris
Subject: Re: Question about a numerical computation (possibly OT as Clisp question)
Date:
Message-ID: <836571el9f.fsf@torus.sehlabs.com>
jblazi <······@hotmail.com> writes:
> The problem is that my epsilon has been introduced artificially. How do I
> perform this computation (after setting long-float-digits) as precisely
> aas possible? I do not think I can use long-float-epsilon here?
Not taking your user-defined float precision into account, but using
scaled-epsilon�:
(defun heron-sqrt (x)
(loop with epsilon = (scaled-epsilon x '-)
and r = 1.0
do (let ((next (/ (+ r (/ x r)) 2)))
(cond ((< (abs (- x (* next next))) epsilon)
(return next))
((< (abs (- r next))
(scaled-epsilon r '-))
(error "Convergence failure."))
(t (setf r next))))))
For (heron-sqrt 2.0) with CLISP, I get a convergence failure. The
function converges on a "perfect" r, but squaring that r does not get
it close enough to x. That is, |x - r^2| >= epsilon. This is
disappointing, but it may not be too surprising, though, for:
CL-USER> (expt (sqrt 2) 2)
1.9999999
CL-USER> (>= (abs (- 2 (expt (sqrt 2) 2)))
(scaled-epsilon 2.0 '-))
T
Footnotes:
� http://groups.google.com/groups?selm=3283568052702575KL2065E%40naggum.no
http://groups.google.com/groups?selm=jk41xj455cn.fsf%40W003275.na.alarismed.com
--
Steven E. Harris
jblazi <······@hotmail.com> writes:
> The following program to compute the squareroot of a number x (I teach
> as some of you know) works of course not very well:
>
>
> (setf (ext:long-float-digits) 20000)
> (setq epsilon 0.0000000000000000000000001L0)
>
> (defun heron-sqrt (x)
> (let ((tmp 1.0L0))
> (loop
> (setq tmp (* 0.5L0 (+ tmp (/ x tmp))))
> (if (< (abs (- (* tmp tmp) x)) epsilon)
> (return tmp)))))
>
> (print (heron-sqrt 2.0L0))
>
> The problem is that my epsilon has been introduced artificially. How do I
> perform this computation (after setting long-float-digits) as precisely
> aas possible? I do not think I can use long-float-epsilon here?
cf. http://www.lispworks.com/reference/HyperSpec/Body/c_number.htm
For example: DOUBLE-FLOAT-EPSILON
But of course, for long-floats on clisp, there is no epsilon, since
you speficy yourself the precision you want:
(defparameter epsilon (expt 0.1 (1- (ext:long-float-digits))))
[Personally, I prefer:
(defun average (&rest args) (/ (apply (function +) args) (length args)))
(defun square (x) (* x x))
(defun double (x) (+ x x))
(defun fixed-point (fun start &optional epsilon trace)
(unless epsilon (setq epsilon 0.001))
(do* ((prev start curr)
(curr (funcall fun prev)
(funcall fun prev)))
((< (abs (- curr prev)) epsilon) curr)
(when trace (format *trace-output* "fixed-point step f(~A)=~A~%" prev curr))
));;fixed-point
CL-USER> (let ((x 9.0))
(fixed-point (lambda (y) (/ (+ y (/ x y)) 2)) (/ x 2) 0.00001 t))
fixed-point step f(4.5)=3.25
fixed-point step f(3.25)=3.0096154
fixed-point step f(3.0096154)=3.0000153
fixed-point step f(3.0000153)=3.0
3.0
]
> Additionally, I do not understand either, where I have to put the "L0" at
> the end of a nuber and when this is not important.
--
__Pascal Bourguignon__ http://www.informatimago.com/
Our enemies are innovative and resourceful, and so are we. They never
stop thinking about new ways to harm our country and our people, and
neither do we.
From: Alain Picard
Subject: Re: Question about a numerical computation (possibly OT as Clisp question)
Date:
Message-ID: <87eklpvxuw.fsf@memetrics.com>
Pascal Bourguignon <····@mouse-potato.com> writes:
>
> [Personally, I prefer:
>
> (defun average (&rest args) (/ (apply (function +) args) (length args)))
A note for newbies who might see this later and cut and paste
this into their code: Don't.
Use something like
(defun average (args)
(assert args () "average not defined without 1 or more arguments")
(/ (reduce (function +) args)
(length args)))
This way you won't croak if you have more than CALL-ARGUMENTS-LIMIT args.
Cheers,
--ap
Alain Picard <············@memetrics.com> writes:
> Pascal Bourguignon <····@mouse-potato.com> writes:
>
> >
> > [Personally, I prefer:
> >
> > (defun average (&rest args) (/ (apply (function +) args) (length args)))
Yes, this was my newbie code...
> A note for newbies who might see this later and cut and paste
> this into their code: Don't.
>
> Use something like
>
> (defun average (args)
> (assert args () "average not defined without 1 or more arguments")
> (/ (reduce (function +) args)
> (length args)))
>
> This way you won't croak if you have more than CALL-ARGUMENTS-LIMIT args.
Note however that CALL-ARGUMENTS-LIMIT minimum is specified to be 50,
so the "newbie" code could lead to problems only on a multi-line
(average ...) sexp, or an (apply (function average) arg-list).
--
__Pascal Bourguignon__ http://www.informatimago.com/
Our enemies are innovative and resourceful, and so are we. They never
stop thinking about new ways to harm our country and our people, and
neither do we.
From: Rob Warnock
Subject: Re: Question about a numerical computation (possibly OT as Clisp question)
Date:
Message-ID: <p7SdnUqfluRWfa7cRVn-oQ@speakeasy.net>
Pascal Bourguignon <····@mouse-potato.com> wrote:
+---------------
| Alain Picard <············@memetrics.com> writes:
| > Pascal Bourguignon <····@mouse-potato.com> writes:
| > > [Personally, I prefer:
| > > (defun average (&rest args) (/ (apply (function +) args) (length args)))
| >
| > ... (/ (reduce (function +) args) (length args))
| > This way you won't croak if you have more than CALL-ARGUMENTS-LIMIT args.
|
| Note however that CALL-ARGUMENTS-LIMIT minimum is specified to be 50,
| so the "newbie" code could lead to problems only on a multi-line
| (average ...) sexp, or an (apply (function average) arg-list).
+---------------
Sounds like a great opportunity for a tutorial on compiler macros!
First, show the REDUCE version [which should always be the default],
then show a compiler macro that examines the length of the arg list
and if it's less than CALL-ARGUMENTS-LIMIT rewrites it into the APPLY
version.
-Rob
-----
Rob Warnock <····@rpw3.org>
627 26th Avenue <URL:http://rpw3.org/>
San Mateo, CA 94403 (650)572-2607
From: Paul Foley
Subject: Re: Question about a numerical computation (possibly OT as Clisp question)
Date:
Message-ID: <m2r7pnx57s.fsf@mycroft.actrix.gen.nz>
On Mon, 30 Aug 2004 21:22:03 -0500, Rob Warnock wrote:
> Sounds like a great opportunity for a tutorial on compiler macros!
> First, show the REDUCE version [which should always be the default],
> then show a compiler macro that examines the length of the arg list
> and if it's less than CALL-ARGUMENTS-LIMIT rewrites it into the APPLY
> version.
If the list is a literal, the compiler macro can produce the result
directly. If not, it can't check the length.
[Unless you're thinking of using &rest and looking at the number of
arguments, but then it must be less than CALL-ARGUMENTS-LIMIT, so you
can always use APPLY!]
--
The State is the great fiction by which everyone seeks to live at the
expense of everyone else. -- Fr�d�ric Bastiat
(setq reply-to
(concatenate 'string "Paul Foley " "<mycroft" '(··@) "actrix.gen.nz>"))
From: Gareth McCaughan
Subject: Re: Question about a numerical computation (possibly OT as Clisp question)
Date:
Message-ID: <87u0ulbev5.fsf@g.mccaughan.ntlworld.com>
jblazi <······@hotmail.com> writes:
> The following program to compute the squareroot of a number x (I teach
> as some of you know) works of course not very well:
>
>
> (setf (ext:long-float-digits) 20000)
> (setq epsilon 0.0000000000000000000000001L0)
>
> (defun heron-sqrt (x)
> (let ((tmp 1.0L0))
> (loop
> (setq tmp (* 0.5L0 (+ tmp (/ x tmp))))
> (if (< (abs (- (* tmp tmp) x)) epsilon)
> (return tmp)))))
>
> (print (heron-sqrt 2.0L0))
>
> The problem is that my epsilon has been introduced artificially. How do I
> perform this computation (after setting long-float-digits) as precisely
> as possible? I do not think I can use long-float-epsilon here?
Roughly, the relative error in y^2 is twice the relative
error in y. If it's the case that the algorithm can be
expected to get within long-float-epsilon of the correct
answer (relatively) but no nearer, then you would want
approximately (* 2 epsilon x) in the termination criterion;
but checking the details is a non-trivial bit of numerical
analysis and it's 2am local time so I shan't attempt it
right now. (The key questions are: how close to the ideal
answer can the iteration actually get, and what's the
effect of the extra roundoff error incurred when calculating
the quantity you look at to test for convergence.)
> Additionally, I do not understand either, where I have to put the "L0" at
> the end of a nuber and when this is not important.
You have to put it at the end of a number when you want
that number to be a long-float rather than whatever kind
of float is the reader's default. (Which will be
SINGLE-FLOAT unless it's told otherwise. Alas.)
When would you want *that*? To answer that question,
you need to understand the rules of float contagion;
the HyperSpec explains these. (I think you are using
CLISP, which at one time had default float contagion
rules incompatible with ANSI's; IIRC, this has been
changed, but you might want to check.)
Oh, there's another reason that doesn't have to do
with float contagion. If the first thing that's going
to be done with a number is to convert it to a long-float
(so that some arithmetic operation can be done between
it and another long-float, perhaps) then you might as
well make it a long-float to begin with; that will be
faster.
--
Gareth McCaughan
.sig under construc
From: Paul Foley
Subject: Re: Question about a numerical computation (possibly OT as Clisp question)
Date:
Message-ID: <m26571z02v.fsf@mycroft.actrix.gen.nz>
On Sun, 29 Aug 2004 13:40:09 +0200, jblazi wrote:
> The following program to compute the squareroot of a number x (I teach
> as some of you know) works of course not very well:
> (setf (ext:long-float-digits) 20000)
> (setq epsilon 0.0000000000000000000000001L0)
> (defun heron-sqrt (x)
> (let ((tmp 1.0L0))
> (loop
> (setq tmp (* 0.5L0 (+ tmp (/ x tmp))))
> (if (< (abs (- (* tmp tmp) x)) epsilon)
> (return tmp)))))
> (print (heron-sqrt 2.0L0))
> The problem is that my epsilon has been introduced artificially. How do I
> perform this computation (after setting long-float-digits) as precisely
> aas possible? I do not think I can use long-float-epsilon here?
Forget epsilon; just look to see whether tmp changed at all. When you
run out of bits, it won't change:
(defun heron-sqrt (x)
(let ((val 1l0))
(loop
(let ((next (* 0.5l0 (+ val (/ x val)))))
(if (= next val)
(return val)
(setq val next))))))
For a more portable version (the adjustable-precision long-float only
exists in CLISP), calculate with rationals, instead, and use either
(RATIONAL X) or (RATIONALIZE X) to avoid float contagion -- of course
then you have to go back to using some "epsilon".
> Additionally, I do not understand either, where I have to put the "L0" at
> the end of a nuber and when this is not important.
If you want the number to be a long-float, and the value of
*READ-DEFAULT-FLOAT-FORMAT* is not CL:LONG-FLOAT, you need it there.
If *READ-DEFAULT-FLOAT-FORMAT* is CL:LONG-FLOAT, you don't need it.
If you don't care what subtype of float you get, you don't need it.
[I'm not sure if that's all you're asking]
--
The State is the great fiction by which everyone seeks to live at the
expense of everyone else. -- Fr�d�ric Bastiat
(setq reply-to
(concatenate 'string "Paul Foley " "<mycroft" '(··@) "actrix.gen.nz>"))
Paul Foley <···@below.invalid> writes:
> [...]
> Forget epsilon; just look to see whether tmp changed at all. When you
> run out of bits, it won't change:
>
> (defun heron-sqrt (x)
> (let ((val 1l0))
> (loop
> (let ((next (* 0.5l0 (+ val (/ x val)))))
> (if (= next val)
> (return val)
> (setq val next))))))
Well, how do you know it won't alternate between two values that
only differ in the least significant bit?
---Vassil.
--
Vassil Nikolov <········@poboxes.com>
Hollerith's Law of Docstrings: Everything can be summarized in 72 bytes.
From: Paul Foley
Subject: Re: Question about a numerical computation (possibly OT as Clisp question)
Date:
Message-ID: <m2r7ppxd06.fsf@mycroft.actrix.gen.nz>
On Mon, 30 Aug 2004 01:15:54 -0400, Vassil Nikolov wrote:
> Paul Foley <···@below.invalid> writes:
>> [...]
>> Forget epsilon; just look to see whether tmp changed at all. When you
>> run out of bits, it won't change:
>>
>> (defun heron-sqrt (x)
>> (let ((val 1l0))
>> (loop
>> (let ((next (* 0.5l0 (+ val (/ x val)))))
>> (if (= next val)
>> (return val)
>> (setq val next))))))
> Well, how do you know it won't alternate between two values that
> only differ in the least significant bit?
In order for that to happen, (+ VAL (/ X VAL)) would have to be
identical to VAL except in that last bit (multiplying by 0.5 only
affects the exponent, not the mantissa), but that result can't be less
than VAL (assuming X is positive), so if the least significant bit of
VAL is set, it can't be clear in NEXT unless they differ in other bits
as well...so it must terminate, right?
--
The State is the great fiction by which everyone seeks to live at the
expense of everyone else. -- Fr�d�ric Bastiat
(setq reply-to
(concatenate 'string "Paul Foley " "<mycroft" '(··@) "actrix.gen.nz>"))
Paul Foley <···@below.invalid> writes:
> On Mon, 30 Aug 2004 01:15:54 -0400, Vassil Nikolov wrote:
>
>> Paul Foley <···@below.invalid> writes:
>>> [...]
>>> Forget epsilon; just look to see whether tmp changed at all. When you
>>> run out of bits, it won't change:
>>>
>>> (defun heron-sqrt (x)
>>> (let ((val 1l0))
>>> (loop
>>> (let ((next (* 0.5l0 (+ val (/ x val)))))
>>> (if (= next val)
>>> (return val)
>>> (setq val next))))))
>
>> Well, how do you know it won't alternate between two values that
>> only differ in the least significant bit?
>
> In order for that to happen, (+ VAL (/ X VAL)) would have to be
> identical to VAL except in that last bit (multiplying by 0.5 only
> affects the exponent, not the mantissa), but that result can't be less
> than VAL (assuming X is positive), so if the least significant bit of
> VAL is set, it can't be clear in NEXT unless they differ in other bits
> as well...so it must terminate, right?
Yes, I should have figured it out myself. Thanks for setting me
straight.
---Vassil.
--
Vassil Nikolov <········@poboxes.com>
Hollerith's Law of Docstrings: Everything can be summarized in 72 bytes.