[R] FAQ 7.31

Peter Dalgaard p.dalgaard at biostat.ku.dk
Fri May 6 09:53:07 CEST 2005


Robin Hankin <r.hankin at noc.soton.ac.uk> writes:

> 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)
> 
> 
> Does the List have any comments to make?

Shouldn't matter in this sort of case, since the added term is bound
to underflow relative to the partial sum. 

The things to watch out for is that mathematical equality doesn't
imply numerical equality (3 * 0.1 != 0.3), and sometimes granularity
effects in correction terms. We've been bitten by code of the type

   x(n+1) = x(n) + f(x(n))

a couple of times in the numerical code, typically in Newton-type
iterations where the correction term alternates by +/- a few units in
the last place, and the program goes into an infinite loop. Optimizing
compilers can introduce such effects in otherwise functioning code by
applying mathematical transformations (see above...) to the code,
rearranging terms, moving terms outside of loops, folding constants,
etc.


-- 
   O__  ---- Peter Dalgaard             Blegdamsvej 3  
  c/ /'_ --- Dept. of Biostatistics     2200 Cph. N   
 (*) \(*) -- University of Copenhagen   Denmark      Ph: (+45) 35327918
~~~~~~~~~~ - (p.dalgaard at biostat.ku.dk)             FAX: (+45) 35327907




More information about the R-help mailing list