Prev: Still learning to think in Lisp...
Next: NYC LOCAL: Tuesday 11 May 2010 Lisp NYC Meet and NYLUG Hack Meet
From: Raymond Toy on 10 May 2010 19:30

What should the value of (expt 2 #c(2d0 0d0)) be? Ccl and sbcl return

#C(4.000000015237235D0 0.0D0)

Cmucl and clisp return #c(4d0 0d0).

This answer is more accurate, but I couldn't find anything in the CLHS
that would say that #c(4d0 0d0) is the right answer. The rule of
float and rational contagion kind of hints at what might happen, but
only applies to floats, not complexes.

What is the expected answer?

The same issue applies to (log 10 10d0). Ccl and sbcl return
1.0000000138867557D0. Cmucl and clisp return 1d0.

Ray
From: Barry Margolin on 10 May 2010 20:55
In article <m1mxw7fi6p.fsf(a)earthlink.net>,
Raymond Toy <rtoy(a)earthlink.net> wrote:

> What should the value of (expt 2 #c(2d0 0d0)) be? Ccl and sbcl return
>
> #C(4.000000015237235D0 0.0D0)
>
> Cmucl and clisp return #c(4d0 0d0).
>
> This answer is more accurate, but I couldn't find anything in the CLHS
> that would say that #c(4d0 0d0) is the right answer. The rule of
> float and rational contagion kind of hints at what might happen, but
> only applies to floats, not complexes.

From the Dictionary entry for EXPT:

] If the base-number is a rational and power-number is an integer, the
] calculation is exact and the result will be of type rational;
] otherwise a floating-point approximation might result.

So whenever the power number is not an integer, the implementation is
allowed to switch to floating point.

>
> What is the expected answer?
>
> The same issue applies to (log 10 10d0). Ccl and sbcl return
> 1.0000000138867557D0. Cmucl and clisp return 1d0.

Floating point is an approximation, and errors can creep in. It looks
like the latter implementations have some special checks for when the
number is an exact power of the base, while the former use the same
algorithm for everything.

--
Barry Margolin, barmar(a)alum.mit.edu
Arlington, MA
*** PLEASE post questions in newsgroups, not directly to me ***
*** PLEASE don't copy me on replies, I'll read them in the group ***
From: Raymond Toy on 10 May 2010 21:32
>>>>> "Barry" == Barry Margolin <barmar(a)alum.mit.edu> writes:

Barry> In article <m1mxw7fi6p.fsf(a)earthlink.net>,
Barry> Raymond Toy <rtoy(a)earthlink.net> wrote:

>>
>> What is the expected answer?
>>
>> The same issue applies to (log 10 10d0). Ccl and sbcl return
>> 1.0000000138867557D0. Cmucl and clisp return 1d0.

Barry> Floating point is an approximation, and errors can creep
Barry> in. It looks like the latter implementations have some
Barry> special checks for when the number is an exact power of the
Barry> base, while the former use the same algorithm for
Barry> everything.

Well, the difference is way bigger than I would have expected. In
fact, I know that sbcl computes (expt 2 #c(2d0 0d0)) as (exp (* #c(2d0
0d0) (log 2))). That (log 2) returns a single-float as expected and
required. This is what causes the loss of accuracy.

Can't say for clisp, but cmucl doesn't do any special checks about the
exact power of the base. Instead, it just basically coerces both args
to expt to the higher precision before computing the result.

Ray
From: Captain Obvious on 11 May 2010 02:56
RT> Well, the difference is way bigger than I would have expected. In
RT> fact, I know that sbcl computes (expt 2 #c(2d0 0d0)) as (exp (* #c(2d0
RT> 0d0) (log 2))). That (log 2) returns a single-float as expected and
RT> required. This is what causes the loss of accuracy.

Well, if you want to work at double accuracy, why don't you explicitly
specify 2 as double-float?

CL-USER> (expt 2 #c(2.0d0 0.0d0))
#C(4.000000015237235d0 0.0d0)
CL-USER> (expt 2.0d0 #c(2.0d0 0.0d0))
#C(4.0d0 0.0d0)

So it is not a practical issue. If you want consistent accuracy, make it
consistent on your side.

As for the theoretic side, yes, I guess it would be nice if SBCL would
analyze both parameters, convert them in a way which produces most accurate
answer and work with them.
But standard doesn't require that (as far as I understand) and SBCL opts not
to do that. That's fine.
This is essentially performance vs. accuracy trade-off, but with unusual
twist:
there is accuracy impact only if programmers do it wrong, but performance
overhead exists for everybody (I guess).

I don't think compiler should penalize people who do it right while helping
people who do it wrong. At least, it is not in SBCL's spirit.


From: Raymond Toy on 11 May 2010 07:36
>>>>> "Captain" == Captain Obvious <udodenko(a)users.sourceforge.net> writes:

RT> Well, the difference is way bigger than I would have expected.
RT> In fact, I know that sbcl computes (expt 2 #c(2d0 0d0)) as
RT> (exp (* #c(2d0 0d0) (log 2))). That (log 2) returns a
RT> single-float as expected and required. This is what causes
RT> the loss of accuracy.

Captain> Well, if you want to work at double accuracy, why don't
Captain> you explicitly specify 2 as double-float?

Yes, I am well aware of how to get the answer I want. My question is
what is the correct answer?

Captain> This is essentially performance vs. accuracy trade-off,
Captain> but with unusual twist:
Captain> there is accuracy impact only if programmers do it
Captain> wrong, but
Captain> performance overhead exists for everybody (I guess).

How do you know there is a performance tradeoff? Did you check to see
if sbcl does something different? AFAICT, sbcl does nothing special
in this case and punts to the generic expt routine which checks the
types of both arguments to figure out what needs to be done. The cost
of returning the accurate answer is very small.

Captain> I don't think compiler should penalize people who do it
Captain> right while helping people who do it wrong. At least, it
Captain> is not in SBCL's spirit.

And who said it was only sbcl that does it this way?

The rule of float and rational contagion says the rational is
converted to a float of the same format. This hints that 2 should be
converted to 2d0. The very last sentence says

In the case of complex numbers, the real and imaginary parts
are effectively handled individually.

Ray
 |  Next  |  Last
Pages: 1 2 3 4 5 6
Prev: Still learning to think in Lisp...
Next: NYC LOCAL: Tuesday 11 May 2010 Lisp NYC Meet and NYLUG Hack Meet