You should not concentrate on the comparison to Python (there are 1
situation were Python performed well: so what?); rather test it again C. I
think this was the motive of the original poster.

My intentiona has been to point out that there are situation were Python is
not necessarily slow (isn't there a study where Python is compared to Common
Lisp and where that magic "80% of C speed" occurs again?); and at the same
time there are1000 situations were Python does not fulfill peoples
expectation (see the Python newsgroup). But there are situations were Python
is fast enough (I can read satellite binary data in Python in a second and
plot the earth data measurements with DISLIN very quickly). But I go not
around and say Python is as fast as C or Fortran; but I do also not say
Common Lisp is slow. I know the arguments from the Forth people: Only 10% of
the mankind is capable of good programming and most of the people are just
stupid and do not understand how to program well in Forth and they always
try to program C in Forth and so and so and so on...(consult your deja news
deposit).

I do not see why people are always that upset when it comes to Common Lisp
and C: If I want C then I must use C (it is hard to circumvent this) and if
I want Lisp then I must use Lisp.

I hope you do not see it as an act of impoliteness but I will close the
discussion here. I am not interested to go into a heated debate; but I would
have no problem to quarrel face to face. Don't get it wrong: for me usnet is
more or like entertainment.

  What nonsense!  If you are so hung up on performance, a simple solution
  is actually to read the file with your favorite speed-freak "language"
  and have it produce something that you can use either directly or read as
  binary floats.  On your own hardware, it is fantastically unlikely that
  two different programming languages use different floating-point formats
  of the same width, so it is safe to use binary floats between programs.

| I just only wrote my experience and nothing more.

  But then people need to be aware of how much you do not want to become
  good at what you do.  You have made this point yourself several times,
  even thinking that being a professional programmer is a bad thing (for
  you, obviously, but phrased so it would appear to be a general insult to
  all professional programmers).  Your experience is that of an unwilling,
  uninterested permanent newbie.  You will therefore never venture beyond
  what the language offers and has optimized for short programs.

  Your grievances are based in ignorance and an unwillingness to understand
  that what Common Lisp offers without much effort is the ability to write
  short, correct programs, _not_ short, fast programs.  From what you
  write, it is impossible _not_ to conclude that you have failed to grasp
  how to write software to begin with, you are only able to make certain
  things "work" to your satisfaction and other things not, where the
  criteria appear to be quite randomly adopted as opposed to acquired
  though a thoughtful process with a particular goal in mind.

| Don't get it wrong: for me usnet is more or like entertainment.

  Of course it is.  People need to keep this in mind, however, when they
  decide to waste their time trying to help you with anything at all, or
  when they decide to waste their time to correct your incorrect and
  ignorant sputtering of unfounded opinions.  Of course, there is much more
  "entertainment value" in annoying people than in trying to think on your
  own. so people should know what you expect to get out of your postings.

///
--
  In a fight against something, the fight has value, victory has none.
  In a fight for something, the fight is a loss, victory merely relief.

I agreed with your post in general, but

Christopher C. Stacy <cst...@dtpq.com> writes:

I would rewrite this sound bite to say:

 If efficiency is the most important thing to you, and if you're
 not willing to do any more than your language requires you to do
 in order to get that efficiency, then Lisp is not for you.

    Erik>   This is why most reasonable implementations of floating-point hardware
    Erik>   offer much more internal precision, so that intermediate results do not
    Erik>   introduce lots of precision-related errors.  If, as a compiler writer,
    Erik>   you are very good at this stuff, you use _only_ the longer internal forms
    Erik>   for as long as you possibly can, storing only single- or double-float
    Erik>   numbers when you have to.  Common Lisp supports a long-float format that
    Erik>   should be mapped to this longer form, which is usually an 80-bit format
    Erik>   instead of the 64-bit double-float format used by default.  (_Computing_
    Erik>   with single-float these days is just plain nuts, even though it may make
    Erik>   sense to store values of that type.)  It appears reasonable to expect
    Erik>   that double-floats will read and print correctly when the underlying
    Erik>   representation used for the computation is such a long-float, but then it
    Erik>   would be wasteful to believe the same for the (internal) long-floats.

But doing this has its own problems.  If all operations were done
with, say, 80-bit floats, rounding of the final result when stored in
a 64-bit float can produce different results than if all the
computations were done with 64-bit floats.

This shows up on x86 implementations of Lisp where sometimes
double-float-epsilon isn't.

Fortunately, most of my numeric work doesn't really care about these
rounding issues. :-)

There is another benefit to using the extra precision:  the regions
of instability in numeric algorithms will be smaller.  Naive users
may be far less likely to get bogus answers.

    Joe> Yes, but it won't be less accurate.  If all the operations done on
    Joe> the 64-bit floats are exact, then the 80-bit floats will be exact, too.
    Joe> Some inexact operations on 64-bit floats will be exact on 80-bit
    Joe> floats, and this can only improve the precision.  Operations that
    Joe> are inexact on 80-bit floats will be inexact on 64-bit floats, but
    Joe> the resulting error will be smaller.

    Joe> There is another benefit to using the extra precision:  the regions
    Joe> of instability in numeric algorithms will be smaller.  Naive users
    Joe> may be far less likely to get bogus answers.

I'm not a numerical analyst and don't pretend to be, but I disagree
with this statement.  If the extra precision was always available, I'd
agree.  But usually you'll have to convert that 80-bit float to a
64-bit float and that's where you lose.

Somewhere on the web there's an article by Kahan about fused
multiply/add operations that some machines provide.  Not exactly the
same, but I think it hints at the subtleties.  Instead of rounding
twice, it rounds once.  He gives some examples where the use of the
fused mac in finding the roots of a quadratic.  In some cases, the
fused mac gives the square root of a negative number instead of the
positive number that would have happened with IEEE rounding.

    Joe> Yes, but it won't be less accurate.  If all the operations
    Joe> done on the 64-bit floats are exact, then the 80-bit floats
    Joe> will be exact, too.  Some inexact operations on 64-bit floats
    Joe> will be exact on 80-bit floats, and this can only improve the
    Joe> precision.  Operations that are inexact on 80-bit floats will
    Joe> be inexact on 64-bit floats, but the resulting error will be
    Joe> smaller.

    Joe> There is another benefit to using the extra precision: the
    Joe> regions of instability in numeric algorithms will be smaller.
    Joe> Naive users may be far less likely to get bogus answers.

    Ray> I'm not a numerical analyst and don't pretend to be, but I
    Ray> disagree with this statement.  If the extra precision was
    Ray> always available, I'd agree.  But usually you'll have to
    Ray> convert that 80-bit float to a 64-bit float and that's where
    Ray> you lose.

    Ray> Somewhere on the web there's an article by Kahan about fused
    Ray> multiply/add operations that some machines provide.

I think you mean
  http://www.cs.berkeley.edu/~wkahan/MxMulEps.pdf
IIRC, Kahan also argues against Java requiring *all* operations to be
performed in 64-bit IEEE arithmetic.

    Ray> Not exactly the same, but I think it hints at the subtleties.
    Ray> Instead of rounding twice, it rounds once.  He gives some
    Ray> examples where the use of the fused mac in finding the roots
    Ray> of a quadratic.  In some cases, the fused mac gives the
    Ray> square root of a negative number instead of the positive
    Ray> number that would have happened with IEEE rounding.

Kahan argues that since Java does not conform to IEEE 754 (it omits
exception flags and directed rounding) it is dangerously inadequate
for floating point.

On page 30 of the paper
http://www.cs.berkeley.edu/~wkahan/JAVAhurt.pdf
Kahan explicitly recommends (in boldface type, no less)
``Use Double Precision.  When the naive program A(x) is run in
arithmetic twice as precise as the data X and the desired result,
it cannot be harmed by roundoff.''

On page 43, Kahan recommends that except for certain unusual situations
such as gargantuan arrays, all floating point arithmetic should be
performed at the widest finite precision that is not too slow or too
narrow.

Matlab is a numerical computation package that is quite popular.
Matlab versions prior to 5.0 used the full precision available
in the hardware (64-bit mantissa) and rounded to double precision
when storing the results in memory.  Upon version 5 on Windows,
Matlab changed to using double-precision (53-bit mantissa) in
the hardware.

Kahan argues that this is a mistake and the extra precision should
be turned back on.  To quote Kahan, ``Why not get different results
that are better results whenever the hardware affords the opportunity?''

This paper also talks about the `Fused Multiply Accumulate' instruction
(FMAC) which computes
   (fp-round (+ s (* r q)))
rather than
   (fp-round (+ s (fp-round (* r q))))

(assume fp-round is a function that maps the exact result into the nearest
result that can be represented in a floating point number).

Kahan notes that FMAC in general produces more accurate answers, but that
there are subtle problems with using it indiscriminately.  In particular,
when evaluating an expression like this:  (- (* u v) (* r q)), do you
want to round the (* u v) or the (* r q) prior to using the FMAC
instruction?  It almost never matters, but if u = r and v = q, FMAC can
produce a non-zero answer whereas rounding both produces a zero.

    Joe> The article to which you refer,
    Joe> ``Matlab's Loss is Nobody's Gain''
    Joe> is available at
    Joe> http://www.cs.berkeley.edu/~wkahan/MxMulEps.pdf

This is new to me.  I haven't looked at Kahan papers in several years.

    Joe> Kahan argues that this is a mistake and the extra precision should
    Joe> be turned back on.  To quote Kahan, ``Why not get different results
    Joe> that are better results whenever the hardware affords the opportunity?''

I have no problem with using as much precision as possible and I
always use doubles when possible.

I do think it is a small problem in lisp where

   (= (1+ double-float-epsilon) 1)

may or may not be T depending on whether there was a store between the
comparison or not.  And if you change double-float-epsilon to be some
other number, there will be cases were a number smaller than
double-float-epsilon that, when added to 1, isn't equal to 1.

  Hm.  Neither Allegro CL, CMUCL, or CLISP manage to distinguish 1d0 from
  (1+ double-float-epsilon), but with an identical epsilon, LispWorks does.

  (scale-float double-float-epsilon (float-precision double-float-epsilon))
  should in theory be identical to (1+ double-float-epsilon), but couriously,
  (1- (scale-float ...)) yields a value twice as big as the epsilon.  This
  is really weird."

///
--
  In a fight against something, the fight has value, victory has none.
  In a fight for something, the fight is a loss, victory merely relief.

(let ((epsilon double-float-epsilon))
  (not (= (float 1 epsilon) (+ (float 1 epsilon) epsilon))))

on a Sparc, HP, SGI, and Powerpc, and they all returned true, as required.

I think that this is a case where too much precision gets you incorrect
accuracy.  I'm fully prepared to take this as a "bug" in our software, where
we don't "dumb down" enough for the architecture, but I do place the blame
on the Intel architecture itself for not following IEEE-754 fp specs.
The workaround for such a bug is to define a new, larger double-float-epsilon
for x86 archiectures only.

Historically, the epsilons have been a problem for us.  We were called
on their values being too large (i.e. not being "the smallest value for
which ...") and in the same 5.0.1 patch which we gave more exact float
printing and reading, we also changed the float epsilons:

5.0.1:

USER(1): (inspect double-float-epsilon)
double-float = 2.220446049250313d-16 [#x3cb00000 00000000] @ #x400cfca
[1i] USER(2): (inspect single-float-epsilon)
single-float = 1.1920929e-7 [#x34000000] @ #x400e192
[2i] USER(3):

6.x:

CL-USER(1): (inspect double-float-epsilon)
double-float = 1.1102230246251568d-16 [#x3ca00000 00000001] @ #x4025fc2
[1i] CL-USER(2): (inspect single-float-epsilon)
single-float = 5.960465e-8 [#x33800001] @ #x4024862
[2i] CL-USER(3):

Unfortunately, though the single-float-epsilon is good for x86,
the double-float-epsilon is not right for x86 because of the 80-bit
internal representation, which messes up the rounding.

Allegro CL 6.1 on Windows with the latest patches seems to have no problems.
(= 1.0d0 (+ 1.0d0 double-float-epsilon))  gives T

The next smaller float, (double-float (/ 4503599627370497 (expt 2 105))),
gives NIL

  That _is_ the problem.  It should return false.

| Did you use double-float-negative-epsilon for the subtractive case?

  Look, I subtract 1 from the equivalent of (1+ double-float-epsilon).

///
--
  In a fight against something, the fight has value, victory has none.
  In a fight for something, the fight is a loss, victory merely relief.

(defun check-epsilon (putative-epsilon)
  (not (= (float 1 putative-epsilon) (+ (float 1 putative-epsilon)
putative-epsilon))))

(check-epsilon putative-epsilon)  => NIL

You can compute double-float-epsilon as follows:

(defun compute-double-float-epsilon ()
  (- 1.0d0 (* (- (/ 4.0d0 3.0d0) 1.0d0) 3.0d0)))

(check-epsilon (compute-double-float-epsilon)) => T

(integer-decode-float double-float-epsilon)
 4503599627370497
-105
1

(integer-decode-float (find-epsilon))
4503599627370496
-104
1

The formula for computing double-float-epsilon is from Kahan.
It works as follows:
  4/3 = 1.0101010101...
  The rounded quotient differs from the exact answer by 1/3 ULP
  (unit in the last place stored).
  Subtracting 1 from the rounded quotient is exact.
  Multiplying this result by 3 is exact, so the result of the
  multiplication differs from 1 by  1 ULP.

CL-USER(1): 2.1510571102112408d-16
2.1510571102112408d-16
CL-USER(2): (integer-decode-float *)
8725724278030336
-105
1
CL-USER(3): (not (= (float 1 **) (+ (float 1 **) **)))
T
CL-USER(4):

I suspect that the correct epsilon for x86 is only a couple of bits
different than the current one, to counteract any too-eager rounding by
the 80-bit hardware.

I wrote a program that searches for the epsilon on my machine.
I get the answer

1.1107651257113995d-16

(integer-decode-float *)
4505798650626049
-105
1

(setq foo 1.1107651257113995d-16)

(not (= (float 1 foo) (+ (float 1 foo) foo)))
T

The next float down is:
(* 4505798650626048 (expt 2 -105))

(setq bar (double-float (* 4505798650626048 (expt 2 -105))))
1.1107651257113993d-16

(not (= (float 1 bar) (+ (float 1 bar) bar)))
NIL

So unless I've blown it again (today's typo-braino factor is rather high),
double-float-epsilon on my machine is
1.1107651257113995d-16

whereas the binding of double-float-epsilon is
1.1102230246251568d-16

Unfortunately, it isn't!  (This is interesting.)

So let me check my work again, in case I'm being dumb.
(setq computed-epsilon   1.1107651257113995d-16)
(setq computed-epsilon-x 1.1107651257113993d-16)

(not (= (float 1 computed-epsilon) (+ (float 1 computed-epsilon)
computed-epsilon)))
T

(not (= (float 1 computed-epsilon-x) (+ (float 1 computed-epsilon-x)
computed-epsilon-x))
NIL

Ok, so far so good.  Let's look at the binary value of computed-epsilon.
(format t "~b" (integer-decode-float computed-epsilon))
10000000000100000000000000000000000000000000000000001

This is just bizarre!  Mathematically, when we add 1 to
computed epsilon, we are performing this operation:

1 + 4505798650626049/40564819207303340847894502572032

Now in order to add these, they need a common denominator,
so converting we get

40564819207303340847894502572032 + 4505798650626049
---------------------------------------------------
         40564819207303340847894502572032

This would be the exact answer if we could represent it.
However, we can't (the numerator has too many digits).
In binary, the numerator is:
10000000000000000000000000000000000000000000000000000
    10000000000100000000000000000000000000000000000000001

but we can only keep 52 digits beyond the most-significant bit.
The 53rd bit is a 1, so we can't tell if we need to round up
or down (1 is half-way) unless we look at more bits.  Since we
round to even, we ought to round down if *all* the remaining
bits are zero, but round up if *any* of them are 1.

Ignoring for just a moment that 1 bit just a few digits down
from the 53rd bit, we actually have what we'd expect.  The
very last bit is 1, so we round up giving us a rounded
numerator of
  10000000000000000000000000000000000000000000000000001

If that very last bit were a zero, we'd round down giving
us a numerator of
  10000000000000000000000000000000000000000000000000000

which would give us a result of 1.0

We'd expect that any precision float epsilon would have
a mantissa of  10000....0001  (with the appropriate number
of zeroes) because the next smaller mantissa ought to be
 rounded down.

But what's with that bit we were ignoring?!  Apparently
the rounding algorithm is ignoring it too!

  The canonical way to get back to an floating point number from the
  integer-decode-float values is a lot simpler and less computationally
  intensive: (scale-float (float <mantissa> 1d0) <exponent>).  In all
  likelihood, it is implemented with a simple hardware-supported conversion
  from large integer mantissa to large floating point value and adjusting
  the exponent rather than building a rational number out of two bignums on
  my way to the double-float.

  Interesting analysis.  I doubt that the above has any influence, so I
  have not worked out whether it had.

///
--
  In a fight against something, the fight has value, victory has none.
  In a fight for something, the fight is a loss, victory merely relief.

    Lieven> I don't know how feasible this is, but one way to keep numerical
    Lieven> analysts happy would be to have an option to save the full 80-bit
    Lieven> internal representations whenever fp registers have to spilled into
    Lieven> memory.

CMUCL used to have an 80-bit long-float.

Another option, implemented in CMUCL, is to allow the user to set the
precision bits in the FP control register.  I think this tells the
hardware to perform rounding at the desired point.  So if you said
round at 53 bits (double-float), all basic ops are rounded to
double-float precision, even if the numbers are 80-bit floats.

I've seen a number of messages that suggest that extra precision
carried by my floating point hardware may be to blame.  For
example:

  "Raymond Toy" <t...@rtp.ericsson.se> wrote in message
  news:4nlmbfkgu2.fsf@rtp.ericsson.se...
  >
  > If all operations were done with, say, 80-bit floats, rounding
  > of the final result when stored in a 64-bit float can produce
  > different results than if all the computations were done with
  > 64-bit floats.

  * Raymond Toy
  | I do think it is a small problem in lisp where
  |
  |    (= (1+ double-float-epsilon) 1)
  |
  | may or may not be T depending on whether there was a store between the
  | comparison or not.

  "Duane Rettig" <du...@franz.com> wrote
  >
  > I think that this is a case where too much precision gets you incorrect
  > accuracy.  I'm fully prepared to take this as a "bug" in our software, where
  > we don't "dumb down" enough for the architecture, but I do place the blame
  > on the Intel architecture itself for not following IEEE-754 fp specs.

  "Duane Rettig" <du...@franz.com> wrote in message news:4ofg9ztkx.fsf@beta.franz.com...
  >
  > It is likely that the reason for such inaccuracies is due to the
  > mixing of in-register (i.e. 80-bit) operations with memory
  > (64-bit) operations.

What's puzzling to me is that Kahan recommends almost always using the widest
floating point precision available (provided it's fast enough), and he has
nothing bad to say about the Intel floating point architecture.  In fact,
in ``The Baleful Effect of Computer Benchmarks upon Applied Mathematics,
Physics and Chemistry''  available at:
http://www.cs.berkeley.edu/~wkahan/ieee754status/baleful.ps

he includes the Intel 386, 486, Pentium and P6 among the conforming
implementations of IEEE 754.

My understanding is that the Intel chips perform all calculations in
IEEE Double Extended format, and that the result may be rounded to
IEEE Double format if stored in memory.  Kahan defends this practice:

    `The idea of an Extended format has been amply vindicated by its
     use in Hewlett-Packard's financial calculators, which all perform
     all arithmetic and financial functions to there more sig. decimals
     than they store or display.'

Kahan cautions against the `stopped clock paradox' where `inferior
computers might get exactly correct results when superior ones get
merely excellent approximations'.

  "Duane Rettig" <du...@franz.com> wrote in message news:4ofg9ztkx.fsf@beta.franz.com...
  > But it seems like they have also seen the same problem, because
  > their new SSE and SSE2 units are 32 and 64 bits, instead of the
  > 80-bit internal representations of the fp unit.

Kahan says ``Without an Extended format some ostensibly straightforward
double computations are prone to malfunction unless carried out in devious
ways known only to experts;  and matrix computations upon vast arrays of
double data degrade too rapidly as increasing dimensions engender
worsened roundoff.''  Narrowing the internal precision of the floating
point unit is never going to make things better.  You can emulate the
narrowing on a higher-precision fp unit by storing to memory between
operations (or using a mode bit if available).  You cannot easily
emulate higher precision on a narrow fp unit.

Now back to double-float-epsilon....
Double-float-epsilon as defined by the hyperspec is `the smallest
positive [double] float ..., such that the following expression is
true when evaluated:
  (not (= (float 1 epsilon) (+ (float 1 epsilon) epsilon)))

The hyperspec does not require any floating point to conform to IEEE
specifications,  but given any floating point specification, there
is a corresponding `float-epsilon' for it.  In particular,
double-float-epsilon for IEEE 754 Double format is:

10000000000000000000000000000000000000000000000000001
times 2^-105

IEEE 754 Single epsilon is
100000000000000000000001
times 2^-47

and Intel's 80-bit internal floats (conforming to IEEE 754
double-extended) epsilon is:
1000000000000000000000000000000000000000000000000000000000000001
times 2^-127

Now what happens if we add 1.0d0 and double-float-epsilon?
There are 8 cases corresponding to whether the operands are
in double format (memory) or extended format (registers),
and whether the answer is left in a register or saved to memory.
When the float is brought in to the floating point unit it
is extended by appending 11 zeroes to the mantissa and adjusting
the exponent appropriately.  (I'll be ignoring the exponents)
So double-float-epsilon in a register looks like this:
1000000000000000000000000000000000000000000000000000100000000000

When we add this to 1.0, the correct sum is
100000000000000000000000000000000000000000000000000001000000000000000000000 00000000000000000000000000000010000
0000000

which, when rounded to fit in the registers is:
1000000000000000000000000000000000000000000000000000010000000000

(so we lose that bit out there at position 105)