So the details of floating-point math is black magic as far as I'm
concerned but I'm hoping maybe someone can explain to me the reason
between decode-float and integer-decode-float (other than the obvious
one that the first value is an float in one case and an integer in the
other). When should you use one or the other?
Related to that, I was wondering why in Allegro:
(= (multiple-value-bind (significand expt sign)
(integer-decode-float least-positive-double-float)
(* (scale-float (float significand 0.0d0) expt) sign))
least-positive-double-float) ==> T
but:
(= (multiple-value-bind (significand expt sign)
(decode-float least-positive-double-float)
(* (scale-float significand expt) sign))
least-positive-double-float) ==> NIL
I was expecting them both to return T.
Then, for grins, I tried the same thing in SBCL and got:
(= (multiple-value-bind (significand expt sign)
(integer-decode-float least-positive-double-float)
(* (scale-float (float significand 0.0d0) expt) sign))
least-positive-double-float) ==> ERROR
and:
(= (multiple-value-bind (significand expt sign)
(decode-float least-positive-double-float)
(* (scale-float significand expt) sign))
least-positive-double-float) ==> T
-Peter
--
Peter Seibel ·····@javamonkey.com
Lisp is the red pill. -- John Fraser, comp.lang.lisp
Peter Seibel <·····@javamonkey.com> writes:
> Then, for grins, I tried the same thing in SBCL and got:
>
> (= (multiple-value-bind (significand expt sign)
> (integer-decode-float least-positive-double-float)
> (* (scale-float (float significand 0.0d0) expt) sign))
> least-positive-double-float) ==> ERROR
Whoops, that's clearly wrong. I have a fix locally that removes the
spurious typecheck in scale-float, making this return T. Apologies.
Christophe
--
http://www-jcsu.jesus.cam.ac.uk/~csr21/ +44 1223 510 299/+44 7729 383 757
(set-pprint-dispatch 'number (lambda (s o) (declare (special b)) (format s b)))
(defvar b "~&Just another Lisp hacker~%") (pprint #36rJesusCollegeCambridge)
Peter Seibel wrote:
> So the details of floating-point math is black magic as far as I'm
> concerned but I'm hoping maybe someone can explain to me the reason
> between decode-float and integer-decode-float (other than the obvious
> one that the first value is an float in one case and an integer in the
> other). When should you use one or the other?
>
> Related to that, I was wondering why in Allegro:
>
> (= (multiple-value-bind (significand expt sign)
> (integer-decode-float least-positive-double-float)
> (* (scale-float (float significand 0.0d0) expt) sign))
> least-positive-double-float) ==> T
>
> but:
>
> (= (multiple-value-bind (significand expt sign)
> (decode-float least-positive-double-float)
> (* (scale-float significand expt) sign))
> least-positive-double-float) ==> NIL
>
> I was expecting them both to return T.
Both *do* return <T> in LispWorks.
CL-USER 5 >
(= (multiple-value-bind (significand expt sign)
(decode-float least-positive-double-float)
(* (scale-float significand expt) sign))
least-positive-double-float)
T
CL-USER 6 >
Carl Taylor
>>>>> "Peter" == Peter Seibel <·····@javamonkey.com> writes:
Peter> So the details of floating-point math is black magic as far as I'm
Peter> concerned but I'm hoping maybe someone can explain to me the reason
Peter> between decode-float and integer-decode-float (other than the obvious
Peter> one that the first value is an float in one case and an integer in the
Peter> other). When should you use one or the other?
When you want a float result or an integer result? Sorry, I don't
have a better answer than that.
But integer-decode-float can be used to output IEEE floats in binary
format, with some care. You can't do that easily with just
decode-float. Also allows you to some exact arithmetic if you want
to.
Peter> Related to that, I was wondering why in Allegro:
Peter> (= (multiple-value-bind (significand expt sign)
Peter> (integer-decode-float least-positive-double-float)
Peter> (* (scale-float (float significand 0.0d0) expt) sign))
Peter> least-positive-double-float) ==> T
Peter> but:
Peter> (= (multiple-value-bind (significand expt sign)
Peter> (decode-float least-positive-double-float)
Peter> (* (scale-float significand expt) sign))
Peter> least-positive-double-float) ==> NIL
Peter> I was expecting them both to return T.
That's my expectation as well. CMUCL returns T for both. What does
(* (scale-float s e) sign) return in the latter case? Just curious.
Ray
Raymond Toy <···········@ericsson.com> writes:
> That's my expectation as well. CMUCL returns T for both. What does
> (* (scale-float s e) sign) return in the latter case? Just curious.
(multiple-value-bind (m e s)
(decode-float least-positive-double-float)
(* (scale-float m e) s)) ==> 0.0d0
FWIW, a binary search reveals that the scale-float/decode-float
identity breaks down between 2.2250738585072014d-308 and
2.225073858507201d-308. That is:
(multiple-value-bind (m e s)
(decode-float 2.2250738585072014d-308)
(* (scale-float m e) s)) ==> 2.2250738585072014d-308
(multiple-value-bind (m e s)
(decode-float 2.225073858507201d-308)
(* (scale-float m e) s)) ==> 0.0d0
As far as I can tell, the identity fails for all double-floats between
2.225073858507201d-308 and least-positive-double-float.
-Peter
--
Peter Seibel ·····@javamonkey.com
Lisp is the red pill. -- John Fraser, comp.lang.lisp
From: ··········@YahooGroups.Com
Subject: Re: decode-float vs integer-decode-float
Date:
Message-ID: <REM-2004jul23-008@Yahoo.Com>
> From: Peter Seibel <·····@javamonkey.com>
> the scale-float/decode-float identity breaks down between
> 2.2250738585072014d-308 and 2.225073858507201d-308.
What you show there are merely DECIMAL printed representations of the
floating-point values, which cannot possibly be exactly equal to those
values. To make your observation meaningful, you need to find out what
the actual binary floating-point value is inside the machine. In
theory, integer-decode-float should provide the needed info. The
following is in CMU Common Lisp 18b.
(setq x1 2.2250738585072014d-308) ;Note: conversion from decimal to binary
(integer-decode-float x1)
4503599627370496
-1074
1
(setq x1i *)
(format nil "~2r" *)
"10000000000000000000000000000000000000000000000000000"
Aha, that's an exact power of two, which means that a slightly smaller
floating-point value will have an exponent of -1075. If -1074 is the
most-negative exponent for which normalized double-floats are possible,
then it makes sense that things would start to break down for any
slightly smaller value.
To verify that we can re-build floating-point values from the result of
integer-decode-float, I tried this:
(scale-float (float x1i 1.0d0) -1074)
(= * x1) ;==> T
Yup, it worked. Now to do the same with your slightly-smaller value:
(setq x2 2.225073858507201d-308) ;Decimal to binary conversion again
(integer-decode-float x2)
9007192812290047
-1075
1
Indeed the exponent is now one less than before.
(setq x2i *)
(format nil "~2r" *)
"11111111111111111111001111111111111111111111111111111"
And the bit pattern is now almost the next power of two, but not really
very close at all. I guess your binary search missed the transition
point by a half a mile! It would be fun to build a floating-point value
exactly equal to that power minus one, i.e. bit pattern of all ones.
But first, let's verify that we can re-build the floating-point value
you actually started with.
(scale-float (float x2i 1.0d0) -1075)
(= * x2) ;==> NIL
Hmm, this is serious. I would think that integer-decode-float would
just pull out the bit pattern of the mantissa and build an integer out
of it, and float would just convert that bit pattern from the integer
back to a floating-point mantissa, so then we'd get the same exact
value re-built? What's wrong? So what floating-point value got re-built
anyway??
(setq x2bad **)
(integer-decode-float x2bad)
9007191201677311
-1075
1
(format nil "~2r" *)
"11111111111111111111000011111111111111111111111111111"
Indeed, if we can believe the result from integer-decode-float, the
mantissa has two 1-bits changed to 1-bits. Bug in CMUCL??
So let's see if we can re-build a floating-point value exactly twice
x2, by using an exponent of -1074 instead of -1075:
(scale-float (float x2i 1.0d0) -1074)
(integer-decode-float *)
9007192812290047
-1074
1
(= * x2i) ;==> T
(format nil "~2r" **)
"11111111111111111111001111111111111111111111111111111"
Yeah, no problem with that bit pattern, providing the exponent is in
the good range, so FLOAT is working fine, just SCALE-FLOAT sometimes
makes a mistake? One more experiment:
(setq x2t2 (scale-float (float x2i 1.0d0) -1074))
(format nil "~2r" (integer-decode-float *))
"11111111111111111111001111111111111111111111111111111"
(scale-float x2t2 -1)
(format nil "~2r" (integer-decode-float *))
"11111111111111111111000011111111111111111111111111111"
If we can trust integer-decode-float to tell us what's really in the
floating-point mantissa, then clearly there's something going wrong
inside scale-float. But can we really trust integer-decode-float??
Does CMUCL provide a way to directly inspect the internal
representation of floating-point values, to compare with what
integer-decode-float is telling us?
Anyway, let's see whether we can construct the exact power of two,
minus one, as mantissa, with exponent -1075. Knowing the problems we're
having above, I'll use exponent -1074 at first, then try to bump it by
-1.
(< (expt 2 52) x2i (expt 2 53)) ;==> T
(setq x3t2 (scale-float (float (+ -1 (expt 2 53)) 1.0d0) -1074))
(format nil "~2r" (integer-decode-float x3t2))
"11111111111111111111111111111111111111111111111111111"
Looks good, we built the desired mantissa.
(setq x3 (scale-float x3t2 -1))
(format nil "~2r" (integer-decode-float x3))
"11111111111111111111001111111111111111111111111111111"
And SCALE-FLOAT seems to have trashed the mantissa again.
(scale-float x3 1)
(format nil "~2r" (integer-decode-float *))
"11111111111111111111001111111111111111111111111111111"
Well, clearly there's a bug in scale-float!!
Regardless of whether integer-decode-float is telling the correct bit
pattern or not, scale-float by -1 followed by +1 should return to the
original value, and it clearly doesn't. Let's try the other way, +1
followed by -1:
(scale-float (scale-float x3 1) -1)
(format nil "~2r" (integer-decode-float *))
"11111111111111111111000011111111111111111111111111111"
Yes, I'd definitely call that a bug in scale-float. If it's impossible
to exactly represent a particular unnormalized floating-point value
then it would make sense for scale-float to give some other value
instead, but here the unnormalized value with mantissa having all ones
except two zero-bits does exist, we started with it this time, and
still we got a different value back, so it wasn't avoiding a value that
can't be represented, no excuse for giving a different value.
Now in your version of CL, instead of that weird pattern with two
additional zero bits, you got all zeroes, a zero mantissa. What version
of CL were you using? Hmm, in your original article you said you're
using Allegro, right? Maybe it also has a big in scale-float, but a
different bug, or maybe it was simply designed to avoid the bug by
declaring that all floating-point values too small to be represented as
normalized values are forced to exactly zero?
I see you also have SBCL available there. Maybe you'll try some of my
experiments above in both Allegro and SBCL to see if you can diagnose
where the discrepancy is creeping in, perhaps a bug or design decision
in scale-float in all three implementations, different in each.
··········@YahooGroups.Com writes:
>> From: Peter Seibel <·····@javamonkey.com>
>> the scale-float/decode-float identity breaks down between
>> 2.2250738585072014d-308 and 2.225073858507201d-308.
>
> What you show there are merely DECIMAL printed representations of the
> floating-point values, which cannot possibly be exactly equal to those
> values. To make your observation meaningful, you need to find out what
> the actual binary floating-point value is inside the machine. In
> theory, integer-decode-float should provide the needed info. The
> following is in CMU Common Lisp 18b.
> [...]
> "11111111111111111111001111111111111111111111111111111"
>
> And the bit pattern is now almost the next power of two, but not really
> very close at all. I guess your binary search missed the transition
> point by a half a mile!
An alternative hypothesis for which there is extremely good support
(i.e. verging on certainty) is that there are bugs in the cmucl 18b
printer and reader for floating point. As you should be aware, it is
possible to ensure print/read consistency between printed decimal
representations of floating point numbers and the numbers themselves,
even if the printed representation in floating point format is not
exactly equal to the floating point value; however, errors in this
will manifest themselves as the kind of discrepancies you are seeing.
Christophe
--
http://www-jcsu.jesus.cam.ac.uk/~csr21/ +44 1223 510 299/+44 7729 383 757
(set-pprint-dispatch 'number (lambda (s o) (declare (special b)) (format s b)))
(defvar b "~&Just another Lisp hacker~%") (pprint #36rJesusCollegeCambridge)
From: ··········@YahooGroups.Com
Subject: Re: decode-float vs integer-decode-float
Date:
Message-ID: <REM-2004jul28-003@Yahoo.Com>
> From: Christophe Rhodes <·····@cam.ac.uk>
> > "11111111111111111111001111111111111111111111111111111"
> > And the bit pattern is now almost the next power of two, but not really
> > very close at all. I guess your binary search missed the transition
> > point by a half a mile!
> An alternative hypothesis for which there is extremely good support
> (i.e. verging on certainty) is that there are bugs in the cmucl 18b
> printer and reader for floating point.
Ah, indeed that would explain why, looking only at the posted decimal
floating-point representations, and what happens when those are read
back into CMUCL, but not having access to whatever numbers he
originally got from the binary search, the bit pattern is so very far
from 2**n - 1. So maybe his binary search was just fine, and internally
he really did have exactly (2**n - 1) * 2**(-thousandsomething), which
even printed best-possible in decimal notation (I suppose I should
manually check that sometime), but reading it back in trashed the
value.
So it's a good thing I factored that problem out of the subsequent
experiments I did, relying on integer-decode-float to view a
floating-point value, and (scale-float (float mantissa 1.0d0) exponent)
to re-build a floating-point value, instead of print and read
respectively. In fact, given the bug I found in scale-float the moment
a normalized value is scaled smaller to force it to be unnormalized,
it's possible that any bug in READ of floating-point values is actually
caused by that bug in scale-float.
Hmm, I guess I'll manually check the posted decimal float value:
2.225073858507201d-308 = (2225073858507201 * 10**(-15)) * 10**(-308)
When read in, internally it's 9007192812290047 * 2**(-1075)
So let's scale everything up by 2**1075 and see if we get that 900...:
2225073858507201 * 2**1075 * 10**(-15) * 10**(-308)
(* 2225073858507201 (expt 2 1075))
90071992547409904492334735...0503456236639879168
I can just look at that and see it doesn't match!
(floor * (expt 10 (+ 15 308)))
9007199254740990
(format nil "~2b" *)
"11111111111111111111111111111111111111111111111111110"
(ceiling *** (expt 10 (+ 15 308)))
9007199254740991
(format nil "~2b" *)
"11111111111111111111111111111111111111111111111111111"
So indeed the binary search, and the decimal print of it which was
posted, were correct, appx. (2**n-1) * 2**(-1075)
But reading that value back into CMUCL hits the bug, apparently caused
by the bug I found in scale-float. So a riddle: Is there any
way in CMUCL to create a double-float whose mantissa is all one-bits
and whose (base2) exponent is exactly -1075? Well, apparently that
original binary search somehow produced that value. But how??
>>>>> "RobertMaas" == RobertMaas <··········@YahooGroups.Com> writes:
RobertMaas> Ah, indeed that would explain why, looking only at the posted decimal
RobertMaas> floating-point representations, and what happens when those are read
RobertMaas> back into CMUCL, but not having access to whatever numbers he
I think this floating-point printer/reader bug is fixed in the latest
CMUCL coming soon to an ftp site near you.
RobertMaas> So indeed the binary search, and the decimal print of it which was
RobertMaas> posted, were correct, appx. (2**n-1) * 2**(-1075)
RobertMaas> But reading that value back into CMUCL hits the bug, apparently caused
RobertMaas> by the bug I found in scale-float. So a riddle: Is there any
RobertMaas> way in CMUCL to create a double-float whose mantissa is all one-bits
RobertMaas> and whose (base2) exponent is exactly -1075? Well, apparently that
RobertMaas> original binary search somehow produced that value. But how??
What's wrong with (scale-float (float #b1<52 or 53>1 1d0) -1075)?
There is (kernel:make-double-float hi lo) where hi is the top 32-bits
of an IEEE double-float and lo is the low 32-bits. I'll leave it as
an exercise for the reader to convert the desired FP number into this
format.
Ray
··········@YahooGroups.Com writes:
> What you show there are merely DECIMAL printed representations of
> the floating-point values, which cannot possibly be exactly equal to
> those values.
Out of context, but 2 divided 10 last time I checked :-)
1.100000000000000088817841970012523233890533447265625-ly y'rs,
mwh
--
I'm okay with intellegent buildings, I'm okay with non-sentient
buildings. I have serious reservations about stupid buildings.
-- Dan Sheppard, ucam.chat (from Owen Dunn's summary of the year)
From: ··········@YahooGroups.Com
Subject: Re: decode-float vs integer-decode-float
Date:
Message-ID: <REM-2004jul30-003@Yahoo.Com>
> > What you show there are merely DECIMAL printed representations of
> > the floating-point values, which cannot possibly be exactly equal to
> > those values.
> From: Michael Hudson <···@python.net>
> Out of context, but 2 divided 10 last time I checked :-)
Which means that any exact binary fraction can be represented exactly
in base ten, if you're willing to use as many digits in the
representation as bits in the origianl binary fraction. But the decimal
representation I was referring to, namely 2.225073858507201d-308,
doesn't have anywhere near enough digits to do that, so it can't
possibly be an exact binary fraction such as an internal floating-point
value.
> 1.100000000000000088817841970012523233890533447265625
There you have those very many extra digits needed to get an exact
representation, not at all like the values that were posted.
(Also you have a value slightly larger than 1.1, rather than a value
slightly smaller than an exact (negative-)power of two, but that
difference is moot. The same principle applies to both.)
So maybe my original wording was slightly confusing, sorry.
Those particular binary values can't be represented exactly
using the number of decimal digits originally shown,
would need *lots* more digits instead.
··········@YahooGroups.Com writes:
> > > What you show there are merely DECIMAL printed representations of
> > > the floating-point values, which cannot possibly be exactly equal to
> > > those values.
> > From: Michael Hudson <···@python.net>
> > Out of context, but 2 divided 10 last time I checked :-)
>
> Which means that any exact binary fraction can be represented exactly
> in base ten, if you're willing to use as many digits in the
> representation as bits in the origianl binary fraction. But the decimal
> representation I was referring to, namely 2.225073858507201d-308,
> doesn't have anywhere near enough digits to do that, so it can't
> possibly be an exact binary fraction such as an internal floating-point
> value.
>
> > 1.100000000000000088817841970012523233890533447265625
>
> There you have those very many extra digits needed to get an exact
> representation, not at all like the values that were posted.
> (Also you have a value slightly larger than 1.1, rather than a value
> slightly smaller than an exact (negative-)power of two, but that
> difference is moot. The same principle applies to both.)
>
> So maybe my original wording was slightly confusing, sorry.
> Those particular binary values can't be represented exactly
> using the number of decimal digits originally shown,
> would need *lots* more digits instead.
Yes, I'd been fighting floating point code when I wrote my article and
was feeling particularly pedantic that day. I'll try to keep that in
check :-)
Cheers,
mwh
PS: I tried to reply to this privately... that's not a very useful
email address you have in your From: line.
--
<MFen> want to write my requirements for me?
<radix> Sure!
<radix> "show a dancing monkey in the about box"
-- from Twisted.Quotes
From: ··········@YahooGroups.Com
Subject: Re: decode-float vs integer-decode-float
Date:
Message-ID: <REM-2004jul31-002@Yahoo.Com>
> From: Michael Hudson <···@python.net>
> PS: I tried to reply to this privately... that's not a very useful
> email address you have in your From: line.
It's my only defense against spam. If you take the time to join the
Yahoo group, you can e-mail me at that address. If you mechanically
collect addresses seen in Usenet articles, to build a mailing list of
30 million people, and blindly mail to them all, and don't spend the
time reading the 20 million non-delivery notices that result, you won't
bother me.
You can try sending e-mail to my regular Yahoo! Mail address, where
your e-mail will be accepted, but it'll be mixed in with over 2500 spam
per month and I'll never see your e-mail because I just don't have time
to browse so many messages looking for a needle in a haystack.
If you know of any free WebMail or mail-forwarding service that blocks
spam, rather than accepting it to a Bulk folder, so if your e-mail is
mistakenly classified as spam you'll know it, let me know. Spam has
killed e-mail as a useful service on the InterNet.
On Tue, 13 Jul 2004 22:23:23 GMT, Peter Seibel <·····@javamonkey.com>
wrote:
>So the details of floating-point math is black magic as far as I'm
>concerned but I'm hoping maybe someone can explain to me the reason
>between decode-float and integer-decode-float (other than the obvious
>one that the first value is an float in one case and an integer in the
>other). When should you use one or the other?
>
I'll take a shot.
Suppose you want to store stuff to a file using WRITE-BYTE.
To store a float you could use INTEGER-DECODE-FLOAT and then
write however many bytes it takes to store the three results
that come back.
But if you used DECODE-FLOAT you get back a float as one of
the values, so now you're into an infinite recursion...
Of course this assumes that the decoded value for every float can be
undecoded back into the float without any loss of information which
I assume is the intent of these instructions -- I've never actually
used them.
JP Massar <······@alum.mit.edu> writes:
> On Tue, 13 Jul 2004 22:23:23 GMT, Peter Seibel <·····@javamonkey.com>
> wrote:
>
>>So the details of floating-point math is black magic as far as I'm
>>concerned but I'm hoping maybe someone can explain to me the reason
>>between decode-float and integer-decode-float (other than the obvious
>>one that the first value is an float in one case and an integer in the
>>other). When should you use one or the other?
>>
>
> I'll take a shot.
>
> Suppose you want to store stuff to a file using WRITE-BYTE. To store
> a float you could use INTEGER-DECODE-FLOAT and then write however
> many bytes it takes to store the three results that come back.
>
> But if you used DECODE-FLOAT you get back a float as one of
> the values, so now you're into an infinite recursion...
So I'm with you this far. Now my question is why would I ever want to
use DECODE-FLOAT. (I'm sure there's some obvious situation--if I did a
bunch of numeric programming it would probably became obvious. I'm
just trying to short circuit that by finding someone who has already
hit the obvious use cases for each function.)
-Peter
--
Peter Seibel ·····@javamonkey.com
Lisp is the red pill. -- John Fraser, comp.lang.lisp
>>>>> "Peter" == Peter Seibel <·····@javamonkey.com> writes:
Peter> So I'm with you this far. Now my question is why would I ever want to
Peter> use DECODE-FLOAT. (I'm sure there's some obvious situation--if I did a
Peter> bunch of numeric programming it would probably became obvious. I'm
Peter> just trying to short circuit that by finding someone who has already
Peter> hit the obvious use cases for each function.)
FWIW, cmucl uses decode-float in its implementation of logb which is
used to compute the sqrt of a complex number accurately. However,
decode-float is used to get the exponent, not the fraction, so
integer-decode-float would have worked as well, except the exponent
would be different.
Ray
Peter Seibel <·····@javamonkey.com> writes:
> So I'm with you this far. Now my question is why would I ever want to
> use DECODE-FLOAT. (I'm sure there's some obvious situation--if I did a
> bunch of numeric programming it would probably became obvious. I'm
> just trying to short circuit that by finding someone who has already
> hit the obvious use cases for each function.)
Adding a lot of numbers :-)
If you have to add up a lot of floating point numbers, perhaps even of
varying magnitudes, you risk loosing precission. When you add a small
number to a big one, it might loose most of its significant digits in
the process, like in
* (+ 1.0d13 pi)
1.0000000000003142d+13
So the thing to do is to sort the numbers, and start adding together
those of smaller magnitude, and then those of larger magnitude. If
there are a _lot_ of numbers, you might divide the set into sets of a
few hundreds, add them up, and then add the results together.
All this #'SORTing costs time, and in particular if there are a *lot*
of numbers. one possible solution to this is to use #'DECODE-FLOAT,
and add numbers according to their magnitude. Like
(defun sum-it-up (huge-list)
(let ((ht (make-hash-table :test #'eql)))
(dolist (x huge-list)
(multiple-value-bind (m e s)
(decode-float x)
(let ((y (get-hash e ht)))
(cond (y (setf (get-hash e ht) (+ y x)))
(T (setf (get-hash e ht) x))))))
;...etc...
))
(up to errors). One can refine this idea a lot more, of course.
Mario S. Mommer <········@yahoo.com> writes:
> Peter Seibel <·····@javamonkey.com> writes:
>> So I'm with you this far. Now my question is why would I ever want to
>> use DECODE-FLOAT. (I'm sure there's some obvious situation--if I did a
>> bunch of numeric programming it would probably became obvious. I'm
>> just trying to short circuit that by finding someone who has already
>> hit the obvious use cases for each function.)
>
> Adding a lot of numbers :-)
>
> If you have to add up a lot of floating point numbers, perhaps even of
> varying magnitudes, you risk loosing precission. When you add a small
> number to a big one, it might loose most of its significant digits in
> the process, like in
>
> * (+ 1.0d13 pi)
>
> 1.0000000000003142d+13
>
> So the thing to do is to sort the numbers, and start adding together
> those of smaller magnitude, and then those of larger magnitude. If
> there are a _lot_ of numbers, you might divide the set into sets of a
> few hundreds, add them up, and then add the results together.
>
> All this #'SORTing costs time, and in particular if there are a *lot*
> of numbers. one possible solution to this is to use #'DECODE-FLOAT,
> and add numbers according to their magnitude. Like
>
> (defun sum-it-up (huge-list)
> (let ((ht (make-hash-table :test #'eql)))
> (dolist (x huge-list)
> (multiple-value-bind (m e s)
> (decode-float x)
> (let ((y (get-hash e ht)))
> (cond (y (setf (get-hash e ht) (+ y x)))
> (T (setf (get-hash e ht) x))))))
> ;...etc...
> ))
>
> (up to errors). One can refine this idea a lot more, of course.
So just to make sure I understood you I filled in the "etc." and got
this (note you can use the default value argument to GETHASH to avoid
the COND):
(defun sum-it-up (huge-list)
(let ((by-exponent (make-hash-table :test #'eql)))
(dolist (x huge-list)
(multiple-value-bind (m e s)
(decode-float x)
(incf (gethash e by-exponent 0) (* m s))))
(loop for e being the hash-key using (hash-value m) of by-exponent
summing (scale-float m e))))
But wouldn't this work pretty much as well:
(defun isum-it-up (huge-list)
(let ((by-exponent (make-hash-table :test #'eql)))
(dolist (x huge-list)
(multiple-value-bind (m e s)
(integer-decode-float x)
(incf (gethash e by-exponent 0) (* m s))))
(loop for e being the hash-key using (hash-value m) of by-exponent
summing (scale-float (float m 0.0d0) e))))
So I'm still trying to understand when to use DECODE-FLOAT rather than
INTEGER-DECODE-FLOAT. I guess one difference between these two
functions is that I have to hard-wire the INTEGER-DECODE-FLOAT one
with a particular type of float (in the argument to FLOAT) instead of
letting normal float type contagion do its thing. So that's a
difference I can maybe imagine caring about.
-Peter
--
Peter Seibel ·····@javamonkey.com
Lisp is the red pill. -- John Fraser, comp.lang.lisp
--=-=-=
Content-Type: text/plain
On 8806 day of my life Peter Seibel wrote:
> Now my question is why would I ever want to use DECODE-FLOAT.
Here is sample of "fast" approximation of SQRT:
(defun approx-sqrt (x)
(declare (type float x))
(multiple-value-bind (m e s) (decode-float x)
(assert (plusp s)) ; May be, ASSERT is not good choice, but...
(scale-float
(* (+ m 1.0)
(if (oddp e)
#.(* 0.5 (sqrt 2.0))
0.5))
(ash e -1))))
CL-USER> (fast-sqrt 63.0)
7.9375
CL-USER> (sqrt 63.0)
7.937254
CL-USER> (fast-sqrt 65.0)
8.529475
CL-USER> (sqrt 65.0)
8.062258
CL-USER> (approx-sqrt 0.1)
0.31819803
CL-USER> (sqrt 0.1)
0.31622776
CL-USER> (approx-sqrt 64.0)
8.485281
Maximal error is at 2^k numbers. For (2^k)-1 it is very precise.
Understanding of the approach is left as exercise :)
Of course, modern CPUs with FPUs has fast built-in SQRT. But when you
don't have FPU, you can use fast approximation (if you don't need
precise answer).
I have wrote the code just now (specially for you, Peter :) and tested
it slightly, so it may be buggy.
--
Ivan Boldyrev
Is 'morning' a gerund?
--=-=-=
Content-Type: application/pgp-signature
-----BEGIN PGP SIGNATURE-----
Version: GnuPG v1.3.6 (GNU/Linux)
iD8DBQBA+hoN4rmsj66VbhcRAvqQAJ42TBXZ3N6wHSFcvF1/b6U97DR1ywCePouU
dvLgvtxfL4d4ER+w50EnjSo=
=CQ+P
-----END PGP SIGNATURE-----
--=-=-=--
--=-=-=
Content-Type: text/plain
On 8809 day of my life Ivan Boldyrev wrote:
> Here is sample of "fast" approximation of SQRT:
>
> (defun approx-sqrt (x)
Warning: I wasn't aware of FLOAT-RADIX. The code assumes
that FLOAT-RADIX always returns 2.
Though it is possible to implement the function in protable manner.
Code in ieee.lisp that I posted before is based on this assumption
too.
--
Ivan Boldyrev
Sorry my terrible English, my native language is Lisp!
--=-=-=
Content-Type: application/pgp-signature
-----BEGIN PGP SIGNATURE-----
Version: GnuPG v1.3.6 (GNU/Linux)
iD8DBQBA+mga4rmsj66VbhcRAtY7AJ0eNvzQ/plncUmGTmDk7WDW+Z6QWwCbBwL0
GxT8xYq/ncrLtRdzO1UFLy0=
=wy73
-----END PGP SIGNATURE-----
--=-=-=--
Peter Seibel <·····@javamonkey.com> writes:
> Now my question is why would I ever want to use DECODE-FLOAT.
It seems that DECODE-FLOAT has better parity with SCALE-FLOAT, just as
frexp() and ldexp() do in C/C++.� There's no need to convert the
integral significand back into a float, and there's no need to clean
up the returned significand to indicate a float value of zero.
As an example of using DECODE-FLOAT, I can share the following code
that I adapted from Erik Naggum's suggestions�:
(defun epsilon (float)
(typecase float
(short-float short-float-epsilon)
(single-float single-float-epsilon)
(double-float double-float-epsilon)
(long-float long-float-epsilon)))
(defun negative-epsilon (float)
(typecase float
(short-float short-float-negative-epsilon)
(single-float single-float-negative-epsilon)
(double-float double-float-negative-epsilon)
(long-float long-float-negative-epsilon)))
(defun scaled-epsilon (float &optional (operation '+))
"Return the smallest number that would return a value different from
FLOAT if OPERATION were applied to FLOAT and this number. OPERATION
should be either + or -, and defauls to +."
(multiple-value-bind (significand exponent)
(decode-float float)
(multiple-value-bind (1.0-significand 1.0-exponent)
(decode-float (float 1.0 float))
(scale-float
(if (and (eq operation '-)
(= significand 1.0-significand))
(negative-epsilon float)
(epsilon float))
(- exponent 1.0-exponent)))))
Footnotes:
� http://www.dinkumware.com/manuals/reader.aspx?b=p/&h=math.html
� http://groups.google.com/groups?selm=3283568052702575KL2065E%40naggum.no
--
Steven E. Harris
From: ··········@YahooGroups.Com
Subject: Re: decode-float vs integer-decode-float
Date:
Message-ID: <REM-2004jul23-007@Yahoo.Com>
> From: Peter Seibel <·····@javamonkey.com>
> maybe someone can explain to me the reason between decode-float and
> integer-decode-float (other than the obvious one that the first value
> is an float in one case and an integer in the other). When should you
> use one or the other?
Presumably you're going to use this information in some further
processing, right? So which form of information would be easiest for
that further processing to use? For example, I wanted to convert a
single floating-point value to the corresponding lower and upper bounds
as rational numbers. decode-float would have been a pain to work with,
but integer-decode-float made it easy: I get the mantissa, which is an
integer, and just subtract or add 1/2 to get an endpoint. (If the
original value is exactly a power of 2, a slightly different
calculation is needed for the lower bound, namely subtracting 1/4, but
I didn't bother for my first draft of the code since even in that case
my wrong lower bound gave a true lower bound.)
On the other hand, if your algorithm wants to normalize a number to be
close to 1.0, for example as a starting point for Newton's method, or
to print the number in "scientific notation", then regular decode-float
would probably be what you want.