Re: [R] exactly representable numbers
Robin Hankin said the following on 9/11/2006 3:52 AM: > Hi > > Given a real number x, I want to know how accurately R can represent > numbers near x. > > In particular, I want to know the infinum of exactly representable > numbers greater than x, and the supremum of exactly representable > numbers > less than x. And then the interesting thing is the difference > between these two. > > > I have a little function that does some of this: > > > f <- function(x,FAC=1.1){ >delta <- x > while(x+delta > x){ >delta <- delta/FAC > } > return(delta*FAC) > } > > But this can't be optimal. > > Is there a better way? > > > I believe this is what .Machine$double.eps is. From ?.Machine double.eps: the smallest positive floating-point number 'x' such that '1 + x != 1'. It equals 'base^ulp.digits' if either 'base' is 2 or 'rounding' is 0; otherwise, it is '(base^ulp.digits) / 2'. See also .Machine$double.neg.eps. Is this what you need? --sundar __ 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 and provide commented, minimal, self-contained, reproducible code.
Re: [R] exactly representable numbers
Hi Sundar thanks for this. But I didn't make it clear that I'm interested in extreme numbers such as 1e300 and 1e-300. Then > f(1e300) [1] 7.911257e+283 is different from 1e300*.Machine$double.eps [I'm interested in the gap between successive different exactly representable numbers right across the IEEE range] rksh On 11 Sep 2006, at 10:50, Sundar Dorai-Raj wrote: > > > Robin Hankin said the following on 9/11/2006 3:52 AM: >> Hi >> Given a real number x, I want to know how accurately R can >> represent numbers near x. >> In particular, I want to know the infimum of exactly representable >> numbers greater than x, and the supremum of exactly representable >> numbers >> less than x. And then the interesting thing is the difference >> between these two. >> I have a little function that does some of this: >> f <- function(x,FAC=1.1){ >>delta <- x >> while(x+delta > x){ >>delta <- delta/FAC >> } >> return(delta*FAC) >> } >> But this can't be optimal. >> Is there a better way? > > I believe this is what .Machine$double.eps is. From ?.Machine > > double.eps: the smallest positive floating-point number 'x' such that > '1 + x != 1'. It equals 'base^ulp.digits' if either 'base' > is 2 or 'rounding' is 0; otherwise, it is > '(base^ulp.digits) > / 2'. > > See also .Machine$double.neg.eps. Is this what you need? > > --sundar -- Robin Hankin Uncertainty Analyst National Oceanography Centre, Southampton European Way, Southampton SO14 3ZH, UK tel 023-8059-7743 __ 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 and provide commented, minimal, self-contained, reproducible code.
Re: [R] exactly representable numbers
On 9/11/2006 6:15 AM, Robin Hankin wrote: > Hi Sundar > > > thanks for this. But I didn't make it clear that I'm interested in > extreme numbers > such as 1e300 and 1e-300. > > Then > > > f(1e300) > [1] 7.911257e+283 > > is different from > > 1e300*.Machine$double.eps > > > [I'm interested in the gap between successive different exactly > representable > numbers right across the IEEE range] I'm not sure the result you're looking for is well defined, because on at least the Windows platform, R makes use of 80 bit temporaries as well as 64 bit double precision reals. I don't know any, but would guess there exist examples of apparently equivalent formulations of your question that give different answers because one uses the temporaries and the other doesn't. But in answer to your question: a bisection search is what I'd use: you start with x+delta > x, and you know x+0 == x, so use 0 and delta as bracketing points. You should be able to find the value in about 50-60 bisections if you start with delta == x, many fewer if you make use of the double.eps value. Here's my version: not tested too much. f <- function(x) { u <- x l <- 0 mid <- u/2 while (l < mid && mid < u) { if (x < x + mid) u <- mid else l <- mid mid <- (l + u)/2 } u } > f(1e300) [1] 7.438715e+283 > 1e300 + 7.438715e+283 > 1e300 [1] TRUE > 1e300 + 7.438714e+283 > 1e300 [1] FALSE Duncan Murdoch > > > > rksh > > > > > > On 11 Sep 2006, at 10:50, Sundar Dorai-Raj wrote: > >> >> Robin Hankin said the following on 9/11/2006 3:52 AM: >>> Hi >>> Given a real number x, I want to know how accurately R can >>> represent numbers near x. >>> In particular, I want to know the infimum of exactly representable >>> numbers greater than x, and the supremum of exactly representable >>> numbers >>> less than x. And then the interesting thing is the difference >>> between these two. >>> I have a little function that does some of this: >>> f <- function(x,FAC=1.1){ >>>delta <- x >>> while(x+delta > x){ >>>delta <- delta/FAC >>> } >>> return(delta*FAC) >>> } >>> But this can't be optimal. >>> Is there a better way? >> I believe this is what .Machine$double.eps is. From ?.Machine >> >> double.eps: the smallest positive floating-point number 'x' such that >> '1 + x != 1'. It equals 'base^ulp.digits' if either 'base' >> is 2 or 'rounding' is 0; otherwise, it is >> '(base^ulp.digits) >> / 2'. >> >> See also .Machine$double.neg.eps. Is this what you need? >> >> --sundar > > -- > Robin Hankin > Uncertainty Analyst > National Oceanography Centre, Southampton > European Way, Southampton SO14 3ZH, UK > tel 023-8059-7743 > > __ > 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 > and provide commented, minimal, self-contained, reproducible code. __ 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 and provide commented, minimal, self-contained, reproducible code.
Re: [R] exactly representable numbers
On Mon, 11 Sep 2006, Duncan Murdoch wrote: > On 9/11/2006 6:15 AM, Robin Hankin wrote: > > Hi Sundar > > > > > > thanks for this. But I didn't make it clear that I'm interested in > > extreme numbers > > such as 1e300 and 1e-300. That's not relevant, unless you are interested in extremely small numbers. > > Then > > > > > f(1e300) > > [1] 7.911257e+283 (That is inaccurate) > > is different from > > > > 1e300*.Machine$double.eps Yes, since 1e300 is not a power of two. However, Sundar is right in the sense that this is an upper bound for normalized numbers. f(1) is .Machine$double.neg.eps, but so it is for all 1 <= x < 2. This gives you the answer: .Machine$double.neg.eps * 2^floor(log2(x)) Similarly for going below (but carefully as you get an extra halving on the powers of two). These results hold for all but denormalized numbers (those below 1e-308). > > [I'm interested in the gap between successive different exactly > > representable > > numbers right across the IEEE range] > > I'm not sure the result you're looking for is well defined, because on > at least the Windows platform, R makes use of 80 bit temporaries as well > as 64 bit double precision reals. I don't know any, but would guess > there exist examples of apparently equivalent formulations of your > question that give different answers because one uses the temporaries > and the other doesn't. Not at R level. For something to get stored in a real vector, it will be a standard 64-bit double. > But in answer to your question: a bisection search is what I'd use: > you start with x+delta > x, and you know x+0 == x, so use 0 and delta as > bracketing points. You should be able to find the value in about 50-60 > bisections if you start with delta == x, many fewer if you make use of > the double.eps value. Here's my version: not tested too much. > > f <- function(x) { >u <- x >l <- 0 >mid <- u/2 >while (l < mid && mid < u) { > if (x < x + mid) u <- mid > else l <- mid > mid <- (l + u)/2 >} >u > } > > > f(1e300) > [1] 7.438715e+283 > > 1e300 + 7.438715e+283 > 1e300 > [1] TRUE > > 1e300 + 7.438714e+283 > 1e300 > [1] FALSE > > > Duncan Murdoch > > > > > > > > > rksh > > > > > > > > > > > > On 11 Sep 2006, at 10:50, Sundar Dorai-Raj wrote: > > > >> > >> Robin Hankin said the following on 9/11/2006 3:52 AM: > >>> Hi > >>> Given a real number x, I want to know how accurately R can > >>> represent numbers near x. > >>> In particular, I want to know the infimum of exactly representable > >>> numbers greater than x, and the supremum of exactly representable > >>> numbers > >>> less than x. And then the interesting thing is the difference > >>> between these two. > >>> I have a little function that does some of this: > >>> f <- function(x,FAC=1.1){ > >>>delta <- x > >>> while(x+delta > x){ > >>>delta <- delta/FAC > >>> } > >>> return(delta*FAC) > >>> } > >>> But this can't be optimal. > >>> Is there a better way? > >> I believe this is what .Machine$double.eps is. From ?.Machine > >> > >> double.eps: the smallest positive floating-point number 'x' such that > >> '1 + x != 1'. It equals 'base^ulp.digits' if either 'base' > >> is 2 or 'rounding' is 0; otherwise, it is > >> '(base^ulp.digits) > >> / 2'. > >> > >> See also .Machine$double.neg.eps. Is this what you need? > >> > >> --sundar > > > > -- > > Robin Hankin > > Uncertainty Analyst > > National Oceanography Centre, Southampton > > European Way, Southampton SO14 3ZH, UK > > tel 023-8059-7743 > > > > __ > > 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 > > and provide commented, minimal, self-contained, reproducible code. > > __ > 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 > and provide commented, minimal, self-contained, reproducible code. > -- Brian D. Ripley, [EMAIL PROTECTED] 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, UKFax: +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 and provide commented, minimal, self-contained, reproducible code.
Re: [R] exactly representable numbers
On 9/11/2006 8:20 AM, Prof Brian Ripley wrote: > On Mon, 11 Sep 2006, Duncan Murdoch wrote: > >> On 9/11/2006 6:15 AM, Robin Hankin wrote: >> > Hi Sundar >> > >> > >> > thanks for this. But I didn't make it clear that I'm interested in >> > extreme numbers >> > such as 1e300 and 1e-300. > > That's not relevant, unless you are interested in extremely small numbers. > >> > Then >> > >> > > f(1e300) >> > [1] 7.911257e+283 > > (That is inaccurate) > >> > is different from >> > >> > 1e300*.Machine$double.eps > > Yes, since 1e300 is not a power of two. However, Sundar is right in the > sense that this is an upper bound for normalized numbers. > > f(1) is .Machine$double.neg.eps, but so it is for all 1 <= x < 2. > This gives you the answer: .Machine$double.neg.eps * 2^floor(log2(x)) I'm not sure what is going wrong, but that is too small (on my machine, at least): > f1 <- function(x) .Machine$double.neg.eps * 2^floor(log2(x)) > f1(1e300) [1] 7.435085e+283 > 1e300 + f1(1e300) == 1e300 [1] TRUE Notice the difference in the 3rd decimal place from the empirical answer from my bisection search below. > > Similarly for going below (but carefully as you get an extra halving on > the powers of two). > > These results hold for all but denormalized numbers (those below 1e-308). > > >> > [I'm interested in the gap between successive different exactly >> > representable >> > numbers right across the IEEE range] >> >> I'm not sure the result you're looking for is well defined, because on >> at least the Windows platform, R makes use of 80 bit temporaries as well >> as 64 bit double precision reals. I don't know any, but would guess >> there exist examples of apparently equivalent formulations of your >> question that give different answers because one uses the temporaries >> and the other doesn't. > > Not at R level. For something to get stored in a real vector, it will be > a standard 64-bit double. I don't think that's a proof, since R level code can call C functions, and there are an awful lot of callable functions in R, but I don't have a counter-example. Duncan Murdoch > >> But in answer to your question: a bisection search is what I'd use: >> you start with x+delta > x, and you know x+0 == x, so use 0 and delta as >> bracketing points. You should be able to find the value in about 50-60 >> bisections if you start with delta == x, many fewer if you make use of >> the double.eps value. Here's my version: not tested too much. >> >> f <- function(x) { >>u <- x >>l <- 0 >>mid <- u/2 >>while (l < mid && mid < u) { >> if (x < x + mid) u <- mid >> else l <- mid >> mid <- (l + u)/2 >>} >>u >> } >> >> > f(1e300) >> [1] 7.438715e+283 >> > 1e300 + 7.438715e+283 > 1e300 >> [1] TRUE >> > 1e300 + 7.438714e+283 > 1e300 >> [1] FALSE >> >> >> Duncan Murdoch >> >> > >> > >> > >> > rksh >> > >> > >> > >> > >> > >> > On 11 Sep 2006, at 10:50, Sundar Dorai-Raj wrote: >> > >> >> >> >> Robin Hankin said the following on 9/11/2006 3:52 AM: >> >>> Hi >> >>> Given a real number x, I want to know how accurately R can >> >>> represent numbers near x. >> >>> In particular, I want to know the infimum of exactly representable >> >>> numbers greater than x, and the supremum of exactly representable >> >>> numbers >> >>> less than x. And then the interesting thing is the difference >> >>> between these two. >> >>> I have a little function that does some of this: >> >>> f <- function(x,FAC=1.1){ >> >>>delta <- x >> >>> while(x+delta > x){ >> >>>delta <- delta/FAC >> >>> } >> >>> return(delta*FAC) >> >>> } >> >>> But this can't be optimal. >> >>> Is there a better way? >> >> I believe this is what .Machine$double.eps is. From ?.Machine >> >> >> >> double.eps: the smallest positive floating-point number 'x' such that >> >> '1 + x != 1'. It equals 'base^ulp.digits' if either 'base' >> >> is 2 or 'rounding' is 0; otherwise, it is >> >> '(base^ulp.digits) >> >> / 2'. >> >> >> >> See also .Machine$double.neg.eps. Is this what you need? >> >> >> >> --sundar >> > >> > -- >> > Robin Hankin >> > Uncertainty Analyst >> > National Oceanography Centre, Southampton >> > European Way, Southampton SO14 3ZH, UK >> > tel 023-8059-7743 >> > >> > __ >> > 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 >> > and provide commented, minimal, self-contained, reproducible code. >> >> __ >> 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 >> and provide commented, minimal, self-contained, reproducible code. >> > __ R-help@stat.math.ethz
Re: [R] exactly representable numbers
On Mon, 11 Sep 2006, Duncan Murdoch wrote: > On 9/11/2006 8:20 AM, Prof Brian Ripley wrote: > > On Mon, 11 Sep 2006, Duncan Murdoch wrote: > > > > > On 9/11/2006 6:15 AM, Robin Hankin wrote: > > > > Hi Sundar > > > > > > > > > > > > thanks for this. But I didn't make it clear that I'm interested in > > > > extreme numbers > > > > such as 1e300 and 1e-300. > > > > That's not relevant, unless you are interested in extremely small numbers. > > > > > > Then > > > > > > > > > f(1e300) > > > > [1] 7.911257e+283 > > > > (That is inaccurate) > > > > > > is different from > > > > > > > > 1e300*.Machine$double.eps > > > > Yes, since 1e300 is not a power of two. However, Sundar is right in the > > sense that this is an upper bound for normalized numbers. > > > > f(1) is .Machine$double.neg.eps, but so it is for all 1 <= x < 2. > > This gives you the answer: .Machine$double.neg.eps * 2^floor(log2(x)) > > I'm not sure what is going wrong, but that is too small (on my machine, at > least): > > > f1 <- function(x) .Machine$double.neg.eps * 2^floor(log2(x)) > > f1(1e300) > [1] 7.435085e+283 > > 1e300 + f1(1e300) == 1e300 > [1] TRUE > > Notice the difference in the 3rd decimal place from the empirical answer from > my bisection search below. I wasn't going into that much detail: just what accuracy are we looking for here? .Machine$double.neg.eps isn't that accurate (as its help page says). > > Similarly for going below (but carefully as you get an extra halving on the > > powers of two). > > > > These results hold for all but denormalized numbers (those below 1e-308). > > > > > > > > [I'm interested in the gap between successive different exactly > > > > representable > > > > numbers right across the IEEE range] > > > > > > I'm not sure the result you're looking for is well defined, because on at > > > least the Windows platform, R makes use of 80 bit temporaries as well as > > > 64 bit double precision reals. I don't know any, but would guess there > > > exist examples of apparently equivalent formulations of your question that > > > give different answers because one uses the temporaries and the other > > > doesn't. > > > > Not at R level. For something to get stored in a real vector, it will be a > > standard 64-bit double. > > I don't think that's a proof, since R level code can call C functions, and > there are an awful lot of callable functions in R, but I don't have a > counter-example. Oh, it is proof: the storage is declared as C double, and extended precision double only exists on the chip. Since R has to look up 'x' whenever it sees the symbol, it looks at the stored value, not on a register. A compiler could do better, but the current code cannot. -- Brian D. Ripley, [EMAIL PROTECTED] 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, UKFax: +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 and provide commented, minimal, self-contained, reproducible code.
Re: [R] exactly representable numbers
On 9/11/2006 9:01 AM, Prof Brian Ripley wrote: > On Mon, 11 Sep 2006, Duncan Murdoch wrote: > >> On 9/11/2006 8:20 AM, Prof Brian Ripley wrote: >> > On Mon, 11 Sep 2006, Duncan Murdoch wrote: >> > >> > > On 9/11/2006 6:15 AM, Robin Hankin wrote: >> > > > Hi Sundar >> > > > >> > > > >> > > > thanks for this. But I didn't make it clear that I'm interested in >> > > > extreme numbers >> > > > such as 1e300 and 1e-300. >> > >> > That's not relevant, unless you are interested in extremely small numbers. >> > >> > > > Then >> > > > >> > > > > f(1e300) >> > > > [1] 7.911257e+283 >> > >> > (That is inaccurate) >> > >> > > > is different from >> > > > >> > > > 1e300*.Machine$double.eps >> > >> > Yes, since 1e300 is not a power of two. However, Sundar is right in the >> > sense that this is an upper bound for normalized numbers. >> > >> > f(1) is .Machine$double.neg.eps, but so it is for all 1 <= x < 2. >> > This gives you the answer: .Machine$double.neg.eps * 2^floor(log2(x)) >> >> I'm not sure what is going wrong, but that is too small (on my machine, at >> least): >> >> > f1 <- function(x) .Machine$double.neg.eps * 2^floor(log2(x)) >> > f1(1e300) >> [1] 7.435085e+283 >> > 1e300 + f1(1e300) == 1e300 >> [1] TRUE >> >> Notice the difference in the 3rd decimal place from the empirical answer from >> my bisection search below. > > I wasn't going into that much detail: just what accuracy are we looking > for here? .Machine$double.neg.eps isn't that accurate (as its help page > says). > >> > Similarly for going below (but carefully as you get an extra halving on the >> > powers of two). >> > >> > These results hold for all but denormalized numbers (those below 1e-308). >> > >> > >> > > > [I'm interested in the gap between successive different exactly >> > > > representable >> > > > numbers right across the IEEE range] >> > > >> > > I'm not sure the result you're looking for is well defined, because on at >> > > least the Windows platform, R makes use of 80 bit temporaries as well as >> > > 64 bit double precision reals. I don't know any, but would guess there >> > > exist examples of apparently equivalent formulations of your question >> > > that >> > > give different answers because one uses the temporaries and the other >> > > doesn't. >> > >> > Not at R level. For something to get stored in a real vector, it will be a >> > standard 64-bit double. >> >> I don't think that's a proof, since R level code can call C functions, and >> there are an awful lot of callable functions in R, but I don't have a >> counter-example. > > Oh, it is proof: the storage is declared as C double, and extended > precision double only exists on the chip. Since R has to look up 'x' > whenever it sees the symbol, it looks at the stored value, not on a > register. A compiler could do better, but the current code cannot. It's a proof of something, but "apparently equivalent formulations" is a vague requirement. I would guess that if I did come up with a counter-example I would have to push the bounds a bit, and it would be questionable whether the formulations were equivalent. For example, it would clearly be outside the bounds for me to write a function which did the entire computation in C (which would likely give a different answer than R gives). But what if some package already contains that function, and I just call it from R? Duncan Murdoch __ 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 and provide commented, minimal, self-contained, reproducible code.
Re: [R] exactly representable numbers
Hi Duncan [snip] On 11 Sep 2006, at 12:12, Duncan Murdoch wrote: > Here's my version: not tested too much. > > f <- function(x) { >u <- x >l <- 0 >mid <- u/2 >while (l < mid && mid < u) { > if (x < x + mid) u <- mid > else l <- mid > mid <- (l + u)/2 >} >u > } > thanks for this. Wouldn't it be a good idea to have some function that returns "the smallest exactly representable number strictly greater than x"? Or does this exist already? Robin Hankin Uncertainty Analyst National Oceanography Centre, Southampton European Way, Southampton SO14 3ZH, UK tel 023-8059-7743 __ 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 and provide commented, minimal, self-contained, reproducible code.
Re: [R] exactly representable numbers
On 9/11/2006 11:01 AM, Robin Hankin wrote: > Hi Duncan > > > [snip] > > On 11 Sep 2006, at 12:12, Duncan Murdoch wrote: > >> Here's my version: not tested too much. >> >> f <- function(x) { >>u <- x >>l <- 0 >>mid <- u/2 >>while (l < mid && mid < u) { >> if (x < x + mid) u <- mid >> else l <- mid >> mid <- (l + u)/2 >>} >>u >> } >> > > > > thanks for this. Wouldn't it be a good idea to have some function > that returns "the smallest exactly representable number strictly > greater than x"? > > Or does this exist already? I don't know if it exists. I wouldn't have much use for it, but if you would, maybe it's worth writing. If you try to do it based on my code above, here are some bugs I've noticed since sending that: - it doesn't work for negative x - it doesn't work for denormal x (e.g. 1.e-310) I'm not sure what goes wrong in the latter case, but it reports numbers which are unnecessarily large: > f(1.e-310) [1] 4.940656e-324 > 1.e-310 == 1.e-310 + 4e-324 [1] FALSE > 1.e-310 == 1.e-310 + 3e-324 [1] TRUE So the right answer is somewhere between 3e-324 and 4e-324, but my function says something bigger. I suspect the best way to do this accurately is to look at the bit patterns of the stored numbers, and add a 1 to the least significant bit, but that's a bit too much work for me. Duncan Murdoch > > > > Robin Hankin > Uncertainty Analyst > National Oceanography Centre, Southampton > European Way, Southampton SO14 3ZH, UK > tel 023-8059-7743 __ 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 and provide commented, minimal, self-contained, reproducible code.
Re: [R] exactly representable numbers
> thanks for this. Wouldn't it be a good idea to have some function > that returns "the smallest exactly representable number strictly > greater than x"? > > Or does this exist already? I didn't read the whole thread, but yes, this exists. It's in the standard c library, header "math.h" and called "nextafter" Shouldn't be too hard to get it into R Greetings Marc __ 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 and provide commented, minimal, self-contained, reproducible code.