Re: [R] FAQ 7.31

From: Prof Brian Ripley <ripley_at_stats.ox.ac.uk>
Date: Fri 06 May 2005 - 18:00:30 EST

On Fri, 6 May 2005, Robin Hankin wrote:

> Jacobi's theta functions crop up here and there in various branches of
> mathematics,
> such as unimodular and elliptic functions.
>
> They have Taylor expansions, but the powers increase as the square of n, as
> in
>
> 1 + z + z^4 + z^9 + z^16 + z^25 + . . .
>
> so they converge very quickly, provided |z|<1
>
> The following toy function shows how I'm implementing these objects. I just
> add terms until
> they make no difference:
>
>
> f <- function(z, maxiter=10){
> out.old <- 1
> for(n in 1:maxiter){
> out.new <- out.old + z^(n*n)
> if(identical(out.new, out.old)){
> return(out.new)
> }
> out.old <- out.new
> }
> stop("not converged")
> }
>
> [NB this is not a real theta function! Real theta functions are more
> complicated. f() just shows the issues]
>
> I worry about the use of identical() here, because I am comparing two floats
> for equality
> as discussed in FAQ 7.31. Perhaps all.equal() would be better:
>
> all.equal(1e99,1e99+1e83, tol=.Machine$double.eps)

Well, probably 2*.Machine$double.eps since each of two numbers could be one bit from the truth.

The registers on ix86 machines have more precision than the representation stored in RAM. So it is possible that out.new has more precision because it has been kept in registers and out.old has been retrieved from RAM and so has less. The result may well be continued looping.

I think the scenario in the previous para is quite unlikely (it would need a very clever C compiler), but it is theoretically possible and may become likely with semi-compiled versions of R. The C-level equivalent is quite common.

-- 
Brian D. Ripley,                  ripley@stats.ox.ac.uk
Professor of Applied Statistics,  http://www.stats.ox.ac.uk/~ripley/
University of Oxford,             Tel:  +44 1865 272861 (self)
1 South Parks Road,                     +44 1865 272866 (PA)
Oxford OX1 3TG, UK                Fax:  +44 1865 272595

______________________________________________
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Received on Fri May 06 18:08:01 2005

This archive was generated by hypermail 2.1.8 : Fri 03 Mar 2006 - 03:31:37 EST