## Re: [Ecls-list] log(complex) inaccurate for some values

 Re: [Ecls-list] log(complex) inaccurate for some values From: Raymond Toy - 2008-04-13 14:03:14 ```Juan Jose Garcia-Ripoll 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+((x-1)*(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 ```

 [Ecls-list] log(complex) inaccurate for some values From: Raymond Toy - 2008-04-11 01:59:56 ```I reported this bug with clisp, but ecl has the same issue: (log #c(1.000000001d0 1d-5)) -> #C(1.0500000857748893d-9 9.999999989666666d-6) But with clisp's long-float, the true answer is closer to: (setf (long-float-digits) 256) (log (coerce #c(1.000000001d0 1d-5) '(complex long-float))) -> #C(1.0500000821378709167497674762372261576369237933558916554150445317294087 7702943L-9 9.9999999896666666683134966416887953589434016059498897741294881863894719391 568L-6) CMUCL returns #C(1.0500000821378709d-9 9.999999989666666d-6) 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 ```
 Re: [Ecls-list] log(complex) inaccurate for some values From: Juan Jose Garcia-Ripoll - 2008-04-13 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 1984-2007 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+1e-5*1i) ans = 1.050000086326139e-09 + 9.999999989666666e-06i On Fri, Apr 11, 2008 at 3:59 AM, Raymond Toy wrote: > I reported this bug with clisp, but ecl has the same issue: > > (log #c(1.000000001d0 1d-5)) -> > > #C(1.0500000857748893d-9 9.999999989666666d-6) > > But with clisp's long-float, the true answer is closer to: > > (setf (long-float-digits) 256) > > (log (coerce #c(1.000000001d0 1d-5) '(complex long-float))) -> > > #C(1.0500000821378709167497674762372261576369237933558916554150445317294087 > 7702943L-9 > > 9.9999999896666666683134966416887953589434016059498897741294881863894719391 > 568L-6) > > CMUCL returns #C(1.0500000821378709d-9 9.999999989666666d-6) > > 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 ```
 Re: [Ecls-list] log(complex) inaccurate for some values From: Raymond Toy - 2008-04-13 14:03:14 ```Juan Jose Garcia-Ripoll 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+((x-1)*(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 ```
 Re: [Ecls-list] log(complex) inaccurate for some values From: Juan Jose Garcia-Ripoll - 2008-04-13 18:42:51 ```On Sun, Apr 13, 2008 at 4:03 PM, Raymond Toy 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 1d-5)))' -eval '(quit)' #C(1.0500000821378709d-9 9.999999989666666d-6) Thanks again for reporting this issue. Juanjo -- Facultad de Fisicas, Universidad Complutense, Ciudad Universitaria s/n Madrid 28040 (Spain) http://juanjose.garciaripoll.googlepages.com ```
 Re: [Ecls-list] log(complex) inaccurate for some values From: Juan Jose Garcia-Ripoll - 2008-04-13 18:45:31 ```I forgot: the date of these changes is CVS 2008-04-13 20:34 Juanjo -- Facultad de Fisicas, Universidad Complutense, Ciudad Universitaria s/n Madrid 28040 (Spain) http://juanjose.garciaripoll.googlepages.com ```
 Re: [Ecls-list] log(complex) inaccurate for some values From: Raymond Toy - 2008-04-14 00:37:37 ```Juan Jose Garcia-Ripoll wrote: > On Sun, Apr 13, 2008 at 4:03 PM, Raymond Toy 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 ```
 Re: [Ecls-list] log(complex) inaccurate for some values From: Juan Jose Garcia-Ripoll - 2008-04-14 08:06:21 ```On Mon, Apr 14, 2008 at 2:37 AM, Raymond Toy 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 ```
 Re: [Ecls-list] log(complex) inaccurate for some values From: Raymond Toy - 2008-04-15 00:11:44 ```Juan Jose Garcia-Ripoll wrote: > On Mon, Apr 14, 2008 at 2:37 AM, Raymond Toy 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 ```