From: Raymond Toy <toy.raymond@gm...>  20080411 01:59:56

I reported this bug with clisp, but ecl has the same issue: (log #c(1.000000001d0 1d5)) > #C(1.0500000857748893d9 9.999999989666666d6) But with clisp's longfloat, the true answer is closer to: (setf (longfloatdigits) 256) (log (coerce #c(1.000000001d0 1d5) '(complex longfloat))) > #C(1.0500000821378709167497674762372261576369237933558916554150445317294087 7702943L9 9.9999999896666666683134966416887953589434016059498897741294881863894719391 568L6) CMUCL returns #C(1.0500000821378709d9 9.999999989666666d6) I'm guessing here, but I suspect ecl is computing the real part as log(sqrt(x^2+y^2)). This loses accuracy if x is near 1 and y is small (or vice versa). Ray 
From: Juan Jose GarciaRipoll <jjgarcia@us...>  20080413 09:01:05

The following does not imply anything. I will look into this, but I presume one needs special functions to handle the high precision case, as in CMUCL. I also wonder how many digits of accuracy this systematically brings in other cases. Juanjo # /afs/ipp/common/bin/matlab v 2007b nojvm < M A T L A B > Copyright 19842007 The MathWorks, Inc. Version 7.5.0.338 (R2007b) August 9, 2007 To get started, type one of these: helpwin, helpdesk, or demo. For product information, visit http://www.mathworks.com. >> format long >> log(1.000000001e0+1e5*1i) ans = 1.050000086326139e09 + 9.999999989666666e06i On Fri, Apr 11, 2008 at 3:59 AM, Raymond Toy <toy.raymond@...> wrote: > I reported this bug with clisp, but ecl has the same issue: > > (log #c(1.000000001d0 1d5)) > > > #C(1.0500000857748893d9 9.999999989666666d6) > > But with clisp's longfloat, the true answer is closer to: > > (setf (longfloatdigits) 256) > > (log (coerce #c(1.000000001d0 1d5) '(complex longfloat))) > > > #C(1.0500000821378709167497674762372261576369237933558916554150445317294087 > 7702943L9 > > 9.9999999896666666683134966416887953589434016059498897741294881863894719391 > 568L6) > > CMUCL returns #C(1.0500000821378709d9 9.999999989666666d6) > > I'm guessing here, but I suspect ecl is computing the real part as > log(sqrt(x^2+y^2)). This loses accuracy if x is near 1 and y is small > (or vice versa). > > Ray  Facultad de Fisicas, Universidad Complutense, Ciudad Universitaria s/n Madrid 28040 (Spain) http://juanjose.garciaripoll.googlepages.com 
From: Raymond Toy <toy.raymond@gm...>  20080413 14:03:14

Juan Jose GarciaRipoll wrote: > The following does not imply anything. I will look into this, but I > presume one needs special functions to handle the high precision case, > as in CMUCL. I also wonder how many digits of accuracy this > systematically brings in other cases. > > The problem is in the calculation of log(sqrt(x^2+y^2)). If x is close to 1 and y is small or vice versa, sqrt(x^2+y^2) is inaccurate. The solution is to change the computation to log(sqrt(x^2+y^2)) = 1/2*log(x^2+y^2) = 1/2*log(1+((x1)*(x+1)+y^2) But don't use the log function to evaluate this. Use log1p(z) = log(1+z). Most C libraries have a log1p function already, but if not, you can just use the series formula since z is small. Ray 
From: Juan Jose GarciaRipoll <jjgarcia@us...>  20080413 18:42:51

On Sun, Apr 13, 2008 at 4:03 PM, Raymond Toy <toy.raymond@...> wrote: > But don't use the log function to evaluate this. Use log1p(z) = > log(1+z). Most C libraries have a log1p function already, but if not, > you can just use the series formula since z is small. Thanks Ray. I did not know about the log1p function. I have changed ECL to support and export it as a lisp function, adding a simplistic implementation for the platforms that do not support it (mingw32 I believe). Now $ ecl eval '(print (log #c(1.000000001d0 1d5)))' eval '(quit)' #C(1.0500000821378709d9 9.999999989666666d6) Thanks again for reporting this issue. Juanjo  Facultad de Fisicas, Universidad Complutense, Ciudad Universitaria s/n Madrid 28040 (Spain) http://juanjose.garciaripoll.googlepages.com 
From: Juan Jose GarciaRipoll <jjgarcia@us...>  20080413 18:45:31

I forgot: the date of these changes is CVS 20080413 20:34 Juanjo  Facultad de Fisicas, Universidad Complutense, Ciudad Universitaria s/n Madrid 28040 (Spain) http://juanjose.garciaripoll.googlepages.com 
From: Raymond Toy <toy.raymond@gm...>  20080414 00:37:37

Juan Jose GarciaRipoll wrote: > On Sun, Apr 13, 2008 at 4:03 PM, Raymond Toy <toy.raymond@...> wrote: > >> But don't use the log function to evaluate this. Use log1p(z) = >> log(1+z). Most C libraries have a log1p function already, but if not, >> you can just use the series formula since z is small. >> > > Thanks Ray. I did not know about the log1p function. I have changed > ECL to support and export it as a lisp function, adding a simplistic > implementation for the platforms that do not support it (mingw32 I > believe). Now > > > My pleasure. Here's an outline for a one possible implementation of log1p. Very likely not as optimized as the C version, but should work with reasonable accuracy. log1p(x) = 2*log1p(x/(1+sqrt(1+x))) This makes x progressively smaller. When x is small enough, you can use the series log1p(x) = 2*sum((x/(2+x))^(2*k+1)/(2*k+1)) for k from 0 to inf. This series converges quite fast, even for x as large as 1. Hope this helps a bit, Ray 
From: Juan Jose GarciaRipoll <jjgarcia@us...>  20080414 08:06:21

On Mon, Apr 14, 2008 at 2:37 AM, Raymond Toy <toy.raymond@...> wrote: > Here's an outline for a one possible implementation of log1p. Very likely > not as optimized as the C version, but should work with reasonable accuracy. Hi, looks nice. The version I found around, and which only activates when the C library does not support log1p (Windows, probably) is one I found googling around, and which compensates the errors in the ordinary logarithm. double log1p(double x) { double u = 1.0 + x; if (u == 1) { return 0.0; } else { return (log(u) * x)/(u  1.0); } } But your implementation seems more stable. The only thing that looks expensive is the step for making "x" small. Wouldn't it suffice to use log1p(x) = series if x <= 1, and log1p(x) = log1p(1/x) + log(x) elsewhere? Juanjo  Facultad de Fisicas, Universidad Complutense, Ciudad Universitaria s/n Madrid 28040 (Spain) http://juanjose.garciaripoll.googlepages.com 
From: Raymond Toy <toy.raymond@gm...>  20080415 00:11:44

Juan Jose GarciaRipoll wrote: > On Mon, Apr 14, 2008 at 2:37 AM, Raymond Toy <toy.raymond@...> wrote: > >> Here's an outline for a one possible implementation of log1p. Very likely >> not as optimized as the C version, but should work with reasonable accuracy. >> > > Hi, looks nice. The version I found around, and which only activates > when the C library does not support log1p (Windows, probably) is one I > found googling around, and which compensates the errors in the > ordinary logarithm. > > double log1p(double x) > { > double u = 1.0 + x; > if (u == 1) { > return 0.0; > } else { > return (log(u) * x)/(u  1.0); > } > } > > But your implementation seems more stable. The only thing that looks > expensive is the step for making "x" small. Wouldn't it suffice to use > log1p(x) = series if x <= 1, and log1p(x) = log1p(1/x) + log(x) > elsewhere? > > Hey, I like your algorithm better. Clean and simple. I don't understand how it works though, but that's ok. The algorithm I gave was from Brent. Making x small takes a little bit of time, but each step the number is halved, so it doesn't take too long. But yeah, you can just use regular old log for x large enough. I don't know where that breakpoint would be. Ray 