>>>>> "MM" == Martin Maechler <maech...@stat.math.ethz.ch>
>>>>> on Mon, 3 Aug 2009 19:30:24 +0200 writes:
>>>>> "HWB" == Hans W Borchers <hwborch...@googlemail.com>
>>>>> on Mon, 3 Aug 2009 13:15:11 +0000 (UTC) writes:
>>>
HWB> Thanks for pointing out the weak point in this
HWB> computation. I tried out your suggestions and they both
HWB> deliver the correct and accurate result.
HWB> But as a general solution this approach is not
HWB> feasible. We want to provide "complex-step derivatives"
HWB> as a new method for computing exact gradients, for
HWB> example in 'numDeriv::grad' as
HWB> grad(fun, x, method="complex-step")
HWB> and we can not reasonably assume that a user prepares
HWB> his function specially for calling this method.
HWB> I tried out other numerical math software besides
HWB> Matlab, that is Octave, Scilab, Euler and others, and
HWB> they all return the same correct value up to 15
HWB> digits. Should we not expect that R is capable of doing
HWB> the same?
MM> R's "a ^ b" on non-Windows typically uses the C library's
MM> 'cpow()' when that is provided (HAVE_C99_COMPLEX).
MM> Indeed, this seems to use the "general" complex power function
MM> which loses a few bits -- unavoidably.
MM> Our Windows-version of complex a ^ b does so as well.
MM> Here's an R version of what (I believe) once was the C library
MM> cpow(); at least I confirm that it has the same slight
MM> inaccuracy [if you are as this very border line case with '1e-15';
MM> use 1e-12 and you have no problems !! ]
MM> Cpow <- function(a,b)
MM> {
MM> ## Purpose: a ^ b in complex
MM> ## Find bug in complex_pow()
MM> ##
-------------------------------------------------------------------------
MM> ## Author: Martin Maechler, Date: 15 Jan 2000, 21:33
MM> a <- as.complex(a)
MM> b <- as.complex(b)
MM> hypot <- function(x,y)Mod(x + 1i*y)
MM> logr <- log(hypot(Re(a), Im(a)) );
MM> logi <- atan2(Im(a), Re(a));
MM> x <- exp( logr * Re(b) - logi * Im(b) );
MM> y <- logr * Im(b) + logi * Re(b);
MM> x * complex(re= cos(y), im= sin(y))
MM> }
MM> ----------------
MM> So, yes we could check for the special case of "^2" and use
MM> multiplication and then for " ^ n " and ... ...
MM> and only otherwise call cpow(x,y) {or our own C-version of
MM> that}.
MM> I'm looking into implenting that now.
MM> Expect to hear more about it within 24 hours.
I have now committed a change to R-devel (not sure if to be
back-ported to R-patched)
which uses ~ log2(n) multiplications for z^n when
n is integer (and |n| <= 2^16 ; that cut-off being somewhat arbitrary).
Along the same line, I've looked what we do for x^y when x,y
are double. Till now, we had only special cased the cases
y == 0 (, 1), 2; and after some simple tests, I've decided to
use the log(n) #{multiplications} algorithm, whenever |n| <= 256.
Thanks also to Tom Short for investigating what other free
packages use.
Martin Maechler, ETH Zurich
HWB> Martin Becker <martin.becker <at> mx.uni-saarland.de>
HWB> writes:
>>>
>>> Dear Ravi,
>>>
>>> the inaccuracy seems to creep in when powers are
>>> calculated. Apparently, some quite general function is
>>> called to calculate the squares, and one can avoid the
>>> error by reformulating the example as follows:
>>>
>>> rosen <- function(x) { n <- length(x) x1 <- x[2:n] x2 <-
>>> x[1:(n-1)] sum(100*(x1-x2*x2)*(x1-x2*x2) + (1-x2)*(1-x2))
>>> }
>>>
>>> x0 <- c(0.0094, 0.7146, 0.2179, 0.6883, 0.5757, 0.9549,
>>> 0.7136, 0.0849,
HWB> 0.4147, 0.4540)
>>> h <- c(1.e-15*1i, 0, 0, 0, 0, 0, 0, 0, 0, 0) xh <- x0 + h
>>>
>>> rx <- rosen(xh) Re(rx) Im (rx)
>>>
>>> I don't know which arithmetics are involved in the
>>> application you mentioned, but writing some auxiliary
>>> function for the calculation of x^n when x is complex and
>>> n is (a not too large) integer may solve some of the
>>> numerical issues. A simple version is:
>>>
>>> powN <- function(x,n) sapply(x,function(x)
>>> prod(rep(x,n)))
>>>
>>> The corresponding summation in 'rosen' would then read:
>>>
>>> sum(100*powN(x1-powN(x2,2),2) + powN(1-x2,2))
>>>
>>> HTH,
>>>
>>> Martin
>>>
MM> ______________________________________________
HWB> R-devel@r-project.org mailing list
HWB> https://stat.ethz.ch/mailman/listinfo/r-devel
MM> ______________________________________________
MM> R-devel@r-project.org mailing list
MM> https://stat.ethz.ch/mailman/listinfo/r-devel
______________________________________________
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
