Re: [R] exactly representable numbers

2006-09-11 Thread Sundar Dorai-Raj


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

2006-09-11 Thread Robin Hankin
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

2006-09-11 Thread Duncan Murdoch
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

2006-09-11 Thread Prof Brian Ripley
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

2006-09-11 Thread Duncan Murdoch
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

2006-09-11 Thread Prof Brian Ripley
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

2006-09-11 Thread Duncan Murdoch
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

2006-09-11 Thread Robin Hankin
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

2006-09-11 Thread Duncan Murdoch
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

2006-09-11 Thread Marc Halbrügge

> 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.