>>>>> "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?
R's "a ^ b" on non-Windows typically uses the C library's
'cpow()' when that is provided (HAVE_C99_COMPLEX).
Indeed, this seems to use the "general" complex power function
which loses a few bits -- unavoidably.
Our Windows-version of complex a ^ b does so as well.
Here's an R version of what (I believe) once was the C library
cpow(); at least I confirm that it has the same slight
inaccuracy [if you are as this very border line case with '1e-15';
use 1e-12 and you have no problems !! ]
Cpow <- function(a,b)
{
## Purpose: a ^ b in complex
## Find bug in complex_pow()
## -------------------------------------------------------------------------
## Author: Martin Maechler, Date: 15 Jan 2000, 21:33
a <- as.complex(a)
b <- as.complex(b)
hypot <- function(x,y)Mod(x + 1i*y)
logr <- log(hypot(Re(a), Im(a)) );
logi <- atan2(Im(a), Re(a));
x <- exp( logr * Re(b) - logi * Im(b) );
y <- logr * Im(b) + logi * Re(b);
x * complex(re= cos(y), im= sin(y))
}
----------------
So, yes we could check for the special case of "^2" and use
multiplication and then for " ^ n " and ... ...
and only otherwise call cpow(x,y) {or our own C-version of
that}.
I'm looking into implenting that now.
Expect to hear more about it within 24 hours.
Martin Maechler, ETH Zurich
HWB> Hans W. Borchers
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
>>
______________________________________________
HWB> R-devel@r-project.org mailing list
HWB> https://stat.ethz.ch/mailman/listinfo/r-devel
______________________________________________
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
