Opened 9 years ago

Closed 9 years ago

#679 closed defect (fixed)

(expt 2 #c(-2d0 -1d0)) inaccurate

Reported by: rtoy Owned by: rme
Priority: normal Milestone:
Component: project hosting Version: 1.4
Keywords: Cc:

Description

CCL 1.4 and 1.5 (32 bit) on Mac OS X 10.5 gives

#C(0.19230972534099303d0 -0.1597403190784087d0)

The correct answer should be closer to

#C(0.19230972430417584 -0.15974031883619208)

as produced by cmucl and clisp.

Perhaps a contagion problem?

Change History (9)

comment:1 Changed 9 years ago by gb

  • Owner set to gb

I'm not really sure whether you mean 32-bit PPC or 32-bit X86 CCL, but in every combination that I tried I got the result that you say is more correct.

My best guess is that you reversed the results above, and meant to say that the more correct answer is:

#C(0.19230972534099303d0 -0.1597403190784087d0)

and are saying that

#C(0.19230972430417584d0 -0.15974031883619208d0)

as returned by CCL is less accurate than CL requires it to be.

I'd agree that the answer returned by

(expt 2.0d0 #c(-2d0 -1d0))

is more accurate, but I'm not sure that I see a sequence of CL contagion rules that lead to the more accurate answer. Do you ?

comment:2 Changed 9 years ago by rtoy

Oops. This is 32-bit x86 ccl.

And I got it backwards. Sorry about that!

The correct answer should be cos(log(2))/4 - i*sin(log(2))/4. CCL says

? (/ (cos (log 2d0)) 4)
0.19230972534099303D0

I don't know of any contagion rule that directly leads to the more accurate answer. CLHS 12.1.4.4 says the result is the max precision of the argument, but doesn't say anything about the precision. It kind of hints that, but doesn't say it.

Note also that (log 2 10d0)}} returns the same answer as (/ (log 2) (log 10d0)) but isn't the same as (log 2d0 10d0).

CCL also produces quite accurate values for (expt 2 (float 1/3 1d0)). So I was expecting the same for the (expt 2 #c(-2d0 -1d0)).

comment:3 Changed 9 years ago by rme

The problem is that in the general case, (expt b e) works by computing (exp (* e (log b))); if b is rational, LOG will return a single-float result.

LOG has a similar issue: it does a change of base (/ (log-e x) (log-e b)) and log-e will return a single-float when passed a rational argument. I'm wondering if it would be better to apply the CL contagion rules in LOG before doing the change of base. (see also ticket:340).

comment:4 follow-up: Changed 9 years ago by gb

I don't see how application of CL contagion rules makes (LOG 2) into (LOG 2.0d0) here. I sort of wish that I did, because I wouldn't want to have to say that if you want a more accurate answer you should ask for one:

(expt 2.0d0 #c(-2d0 -1d0)) 

rather than depending on that happening as an artifact of the algorithm used by EXPT.

comment:5 Changed 9 years ago by rme

The contagion rules (http://www.lispworks.com/documentation/HyperSpec/Body/12_aab.htm) are not that detailed when it comes to complex values, so maybe I am making unwarranted assumptions.

Given: (expt 2 #c(-2d0 -1d0))

by 12.1.5.2, this becomes: (expt #c(2 0) #c(-2d0 -1d0))

Although I do not see it stated in the spec, I would expect that the rational complex would be converted to a (complex double-float) in this example, by analogy to 12.1.4.1.

If this conversion is done in EXPT (a "primitive" CL function), then the implementation, which involves calling LOG, will get the right type of argument.

That's how I was looking at it, anyway.

comment:6 Changed 9 years ago by rme

At the risk of making myself look more foolish than usual, what about something like this patch?

Index: l0-float.lisp
===================================================================
--- l0-float.lisp	(revision 13688)
+++ l0-float.lisp	(working copy)
@@ -839,6 +839,16 @@
            #+64-bit-target
            (%single-float-expt (%short-float b) (%short-float e))
            ))
+	((typep (realpart e) 'double-float)
+	 (let ((upgraded-base (if (complexp b)
+				(coerce b '(complex double-float))
+				(float b 0d0))))
+	   (exp (* e (log upgraded-base)))))
+	((typep (realpart b) 'double-float)
+	 (let ((upgraded-exponent (if (complexp e)
+				    (coerce e '(complex double-float))
+				    (float e 0d0))))
+	   (exp (* upgraded-exponent (log b)))))
         (t (exp (* e (log b))))))

comment:7 in reply to: ↑ 4 Changed 9 years ago by rtoy

Replying to gb:

I don't see how application of CL contagion rules makes (LOG 2) into (LOG 2.0d0) here. I sort of wish that I did, because I wouldn't want to have to say that if you want a more accurate answer you should ask for one:

(expt 2.0d0 #c(-2d0 -1d0)) 

rather than depending on that happening as an artifact of the algorithm used by EXPT.

You're definitely not allowed to convert (log 2) to (log 2d0) in general. But for the two arg log, it's not so clear. You're not required to implement (log x base) exactly as (/ (log x) (log base)).

comment:8 Changed 9 years ago by rme

  • Owner changed from gb to rme

comment:9 Changed 9 years ago by rme

  • Resolution set to fixed
  • Status changed from new to closed

(In [13769]) In EXPT, try to avoid loss of precision by promoting the type of the base to a double-float when the exponent is a DOUBLE-FLOAT or (COMPLEX DOUBLE-FLOAT).

Fixes ticket:679.

Note: See TracTickets for help on using tickets.