[ clisp-Bugs-1934968 ] log(complex) sometimes inaccurate

[SourceForge.net <no...@so...>]

Bugs item #1934968, was opened at 2008-04-04 22:51
Message generated for change (Comment added) made by rtoy
You can respond by visiting: 
https://sourceforge.net/tracker/?func=detail&atid=101355&aid=1934968&group_id=1355

Please note that this message will contain a full copy of the comment thread,
including the initial issue submission, for this request,
not just the latest update.
Category: None
Group: None
Status: Open
Resolution: None
Priority: 5
Private: No
Submitted By: Raymond Toy (rtoy)
Assigned to: Nobody/Anonymous (nobody)
Summary: log(complex) sometimes inaccurate

Initial Comment:
(log #c(1.000000001d0 1d-5)) -> 
#C(1.0500000863261394d-9 9.999999989666666d-6)

But consider:

(setf (long-float-digits) 256)

(log (coerce #c(1.000000001d0 1d-5) '(complex long-float))) ->

#C(1.05000008213787091674976747623722615763692379335589165541504453172940877702943L-9
   9.9999999896666666683134966416887953589434016059498897741294881863894719391568L-6)

so the real part is a bit off.

----------------------------------------------------------------------

>Comment By: Raymond Toy (rtoy)
Date: 2008-04-07 14:56

Message:
Logged In: YES 
user_id=28849
Originator: YES

That's because sbcl uses cmucl's implementation (I think), which, in turn,
comes from an article by Kahan.

Nothing says you should use naive implementations of the special
functions. :-)

Clisp's sinh(x) isn't consistent either:  (= (sinh x) (* 1/2 (- (exp x)
(exp (- x))))) returns NIL for small values of x. :-)

----------------------------------------------------------------------

Comment By: Sam Steingold (sds)
Date: 2008-04-07 13:46

Message:
Logged In: YES 
user_id=5735
Originator: NO

at least clisp is consistent:
(= (realpart (log #c(1.000000001d0 1d-5))) (log (abs #c(1.000000001d0
1d-5)))) ==> T
as opposed to sbcl which returns NIL.
:-)

----------------------------------------------------------------------

Comment By: Raymond Toy (rtoy)
Date: 2008-04-07 13:36

Message:
Logged In: YES 
user_id=28849
Originator: YES

I think this is only a problem if the realpart is near 1 and the imagpart
is small (or vice versa).

The computation of log(abs(z)) loses precision.  It might be better
computed as 

1/2*log(x^2+y^2) = 1/2*log(1+((x-1)*(x+1)+y^2), using a special log(1+x)
function for small x.

----------------------------------------------------------------------

You can respond by visiting: 
https://sourceforge.net/tracker/?func=detail&atid=101355&aid=1934968&group_id=1355