Hi,
I hope this is the right group to post in; if it isn't, redirections are
very welcome.
My problem is this: I've tried computing the value of pi using the
bisection method (for math class, as you might have guessed), and tried to
stop after computing a certain number of significant digits. The code to
check this looks something like this (where 'f' is #'sin):
(if (< (- (funcall f b) (funcall f a)) 1L-7)
m
(if (< (funcall f m) 0)
(bisect m b f)
(bisect a m f)))
That is, if my boundaries a and b differ by a slight enough margin, the
function should terminate. This works nicely for values larger than
2.3841858E-7, but fails for smaller ones. I am guessing that this is some
kind of precision limitation for a type I am using, but I am not able to
see where this might be the case. I'd be grateful if somebody could
enlighten me in this respect... thanks in advance :)
--
Bye: Andi S.
Andi Scharfstein <····@synchron.org> writes:
> (if (< (- (funcall f b) (funcall f a)) 1L-7)
> m
> (if (< (funcall f m) 0)
> (bisect m b f)
> (bisect a m f)))
>
> That is, if my boundaries a and b differ by a slight enough margin, the
> function should terminate. This works nicely for values larger than
> 2.3841858E-7, but fails for smaller ones. I am guessing that this is some
> kind of precision limitation for a type I am using, but I am not able to
> see where this might be the case. I'd be grateful if somebody could
> enlighten me in this respect... thanks in advance :)
To clarify what others might be trying to say, you are using 10^-7 as
your "stop computing" difference. Of course, if you start off with
numbers that are close to 10^-7, you'll stop right away.
Also, why do you think you're guaranteed that (funcall f b) is always
greater than (funcall f a)? (There might be constraints on f you didn't
mention, but you should be cognizant of the problem even then.)
--
Rahul Jain
·····@nyct.net
Professional Software Developer, Amateur Quantum Mechanicist
One thing that bothers me about this approach is that
the you are not measuring the difference between
a and b. Might it be the case that b becomes
very close to a, but in a region where the function f
is very steep? If f is sufficiently steep you might
not be able to get (- (funcall f b) (funcall f a))
less than epsilon.
|b - a| might become much less than epsilon before
|f(b) - f(a)| becomes less than epsilon.
> Andi Scharfstein <····@synchron.org> writes:
>
>
>> (if (< (- (funcall f b) (funcall f a)) 1L-7)
>> m
>> (if (< (funcall f m) 0)
>> (bisect m b f)
>> (bisect a m f)))
>>
>>That is, if my boundaries a and b differ by a slight enough margin, the
>>function should terminate. This works nicely for values larger than
>>2.3841858E-7, but fails for smaller ones. I am guessing that this is some
>>kind of precision limitation for a type I am using, but I am not able to
>>see where this might be the case. I'd be grateful if somebody could
>>enlighten me in this respect... thanks in advance :)
>
>
> To clarify what others might be trying to say, you are using 10^-7 as
> your "stop computing" difference. Of course, if you start off with
> numbers that are close to 10^-7, you'll stop right away.
>
> Also, why do you think you're guaranteed that (funcall f b) is always
> greater than (funcall f a)? (There might be constraints on f you didn't
> mention, but you should be cognizant of the problem even then.)
>
On Mon, 17 Jan 2005 07:47:07 +0100, <·····@rdrop.com> wrote:
> One thing that bothers me about this approach is that
> the you are not measuring the difference between
> a and b. Might it be the case that b becomes
> very close to a, but in a region where the function f
> is very steep? If f is sufficiently steep you might
> not be able to get (- (funcall f b) (funcall f a))
> less than epsilon.
This comment reminded me of when I had to do such problems. There are
other cases, i.e. where f is not flat between a and b, but f(a)
happens to equal f(b). Several more checks are needed in this
approach.
>
>|b - a| might become much less than epsilon before
>|f(b) - f(a)| becomes less than epsilon.
>
>
>> Andi Scharfstein <····@synchron.org> writes:
>>
>>
>>> (if (< (- (funcall f b) (funcall f a)) 1L-7)
>>> m
>>> (if (< (funcall f m) 0)
>>> (bisect m b f)
>>> (bisect a m f)))
>>>
>>>That is, if my boundaries a and b differ by a slight enough margin, the
>>>function should terminate. This works nicely for values larger than
>>>2.3841858E-7, but fails for smaller ones. I am guessing that this is some
>>>kind of precision limitation for a type I am using, but I am not able to
>>>see where this might be the case. I'd be grateful if somebody could
>>>enlighten me in this respect... thanks in advance :)
--
Everyman has three hearts;
one to show the world, one to show friends, and one only he knows.
Hi,
> This comment reminded me of when I had to do such problems. There are
> other cases, i.e. where f is not flat between a and b, but f(a)
> happens to equal f(b). Several more checks are needed in this
> approach.
True, I didn't think of that (since I didn't write the procedure to deal
with the general case, but only the sine function). Maybe a better approach
would be to check whether -1E-7 < f(m) < 1E-7, with m := (a+b)/2? This
might cause problems with functions that approximate 0 very closely, but
don't actually reach it, but it's the best I can think of. How did your
solution to this problem look like?
--
Bye: Andi S.
On 18 Jan 2005 14:35:36 GMT, <····@synchron.org> wrote:
> Hi,
>
>> This comment reminded me of when I had to do such problems. There are
>> other cases, i.e. where f is not flat between a and b, but f(a)
>> happens to equal f(b). Several more checks are needed in this
>> approach.
>
> True, I didn't think of that (since I didn't write the procedure to deal
> with the general case, but only the sine function). Maybe a better approach
> would be to check whether -1E-7 < f(m) < 1E-7, with m := (a+b)/2? This
> might cause problems with functions that approximate 0 very closely, but
> don't actually reach it, but it's the best I can think of. How did your
> solution to this problem look like?
One of the features of the sine function is that it is periodic, and
for every b= a + n*2pi, f(a) = f(b).
The solution to the general problem is to use a different technique,
due to the pitfalls we have mentioned.
--
Everyman has three hearts;
one to show the world, one to show friends, and one only he knows.
Hi,
> One thing that bothers me about this approach is that
> the you are not measuring the difference between
> a and b. Might it be the case that b becomes
> very close to a, but in a region where the function f
> is very steep? If f is sufficiently steep you might
> not be able to get (- (funcall f b) (funcall f a))
> less than epsilon.
>
>|b - a| might become much less than epsilon before
>|f(b) - f(a)| becomes less than epsilon.
This is true, but I don't think it applies in this particular case -
|b - a| is definitely still larger than epsilon when |f(b) - f(a)| stops
changing (and by this I mean changing at all, no change even over large
intervals).
Anyway, can you think of any better way to measure when I have to stop? I
think I have an error approximation somewhere for the bisection method
which I might be able to use, but I'd still be happier to know that my
current problem can be resolved as well.
--
Bye: Andi S.
start quoting Andi Scharfstein :
> That is, if my boundaries a and b differ by a slight enough margin, the
> function should terminate. This works nicely for values larger than
> 2.3841858E-7, but fails for smaller ones. I am guessing that this is some
> kind of precision limitation for a type I am using, but I am not able to
> see where this might be the case. I'd be grateful if somebody could
> enlighten me in this respect... thanks in advance :)
>
Common Lisp uses single-precision floats by default, so you're probably in
the situation that you're comparing against a value that's smaller than
epsilon.
short-float-epsilon is 5.960465e-8 in SBCL, which is in the right ballpark;
floating-point precision decreases as your values increase, and Pi is quite
far away from epsilon/0.
I'm not entirely sure how to fix this, although others probably will be.
You could check whether a and b are equal, which might work though I'm not
sure it's defined. You could also get higher precision by using
double-floats.
Hi,
On 2005-01-17 01:28:34 +0100, Svein Ove Aas <·········@aas.no> said:
> You could also get higher precision by using double-floats.
Yes, I think this is my real question: How do I get Lisp to use higher
precision numbers in the relevant computation? I've tried
'ration'alizing them, but this didn't yield the desired outcome. What
leaves me a little stumped is that the precision *does* get higher with
more iterations, i.e. I get an arbitrarily close approximation to pi -
it's just the comparison that does not seem to want to work, even
though the necessary precision could be provided. Anybody know why this
might be the case?
--
Bye: Andi S.
start quoting Andi Scharfstein :
> Hi,
> On 2005-01-17 01:28:34 +0100, Svein Ove Aas <·········@aas.no> said:
>
>> You could also get higher precision by using double-floats.
>
> Yes, I think this is my real question: How do I get Lisp to use higher
> precision numbers in the relevant computation? I've tried
> 'ration'alizing them, but this didn't yield the desired outcome. What
> leaves me a little stumped is that the precision *does* get higher with
> more iterations, i.e. I get an arbitrarily close approximation to pi -
> it's just the comparison that does not seem to want to work, even
> though the necessary precision could be provided. Anybody know why this
> might be the case?
>
Is there any chance you might be using rationals for the computation and
floating-point for the comparison?
That aside, for many/most functions, if one of the arguments is
double-precision, the result will also be. To get double precision, you can
type the inputs as "Nd0" to get the number N in double-precision floating
point.
You can also programmatically convert a real number to a given
floating-point precision by passing it to the float function with an
appropriate template, like so: (float N 1.0d0). In this case, the value of
the template doesn't matter; only its type does.
I'm not sure how to convert complex numbers like this; it's probably in the
hyperspec somewhere.
Lastly, you can change *read-default-float-format* to double-float before
reading/compiling your source files, in which case you get single-floats by
using "1.0f0" or similar. I'm not sure whether that's a good idea, though.
Hi,
> Is there any chance you might be using rationals for the computation
> and floating-point for the comparison?
Indeed, that's what was causing the error, as I conveniently found out just
before reading your post (by playing around with the interpreter for about
half an hour)... nevertheless, a big thank you to all who replied!
Your input was very helpful, and I'm sure the parts which were not directly
applicable to my current problem will come in handy at some later time :)
--
Bye: Andi S.
Andi Scharfstein wrote:
> Hi,
> I hope this is the right group to post in; if it isn't, redirections are
> very welcome.
>
> My problem is this: I've tried computing the value of pi using the
> bisection method (for math class, as you might have guessed), and tried to
> stop after computing a certain number of significant digits. The code to
> check this looks something like this (where 'f' is #'sin):
>
> (if (< (- (funcall f b) (funcall f a)) 1L-7)
> m
> (if (< (funcall f m) 0)
> (bisect m b f)
> (bisect a m f)))
>
> That is, if my boundaries a and b differ by a slight enough margin, the
> function should terminate. This works nicely for values larger than
> 2.3841858E-7, but fails for smaller ones. I am guessing that this is some
> kind of precision limitation for a type I am using, but I am not able to
> see where this might be the case. I'd be grateful if somebody could
> enlighten me in this respect... thanks in advance :)
>
Can you be certain that f(b) >= f(a) always? If not, then
f(b)-f(a)<0<1L-7...
Furthermore, as you've discovered, you want a relative epsilon not an
absolute epsilon. Try this:
(let* ((fa (funcall f a))
(fb (funcall f b))
(epsilon (* 1L-7 (if (zerop fa) 1 (abs fa)))) )
(if (< (abs (- fa fb)) epsilon) ...
That should work for f(a) <= 0.
David Sletten
Hi,
sorry if this post comes out formatted all wrong, I'm forced to using
Google Groups for posting...
> Can you be certain that f(b) >= f(a) always? If not, then
> f(b)-f(a)<0<1L-7...
Yes, this is an invariant that is maintained throughout the course of
the function call (possibly excepting the very first one, that's when
it gets sorted out by replacing f by (lambda (x) (- (funcall f x))) if
necessary).
> Furthermore, as you've discovered, you want a relative epsilon not an
> absolute epsilon. Try this:
> (let* ((fa (funcall f a))
> (fb (funcall f b))
> (epsilon (* 1L-7 (if (zerop fa) 1 (abs fa)))) )
> (if (< (abs (- fa fb)) epsilon) ...
>
> That should work for f(a) <= 0.
I've tried that quickly and have not yet got it working, I'll try it
again when I get some free time later this day. A quick note to some
posters (I'll try replying individally later on) who seem to think
that my program terminates to quickly, the opposite is the case -
infinite (well, stack overflow) recursion.
OK, gotta go now, I'll be back later. Thanks so far for all your
replies!
--
Bye: Andi S.