Bug 14768 - Improve accuracy of R_pow_di
Summary: Improve accuracy of R_pow_di
Status: CLOSED FIXED
Alias: None
Product: R
Classification: Unclassified
Component: Accuracy (show other bugs)
Version: R 2.14.1
Hardware: All All
: P1 minor
Assignee: R-core
URL:
Depends on:
Blocks:
 
Reported: 2011-12-30 22:31 UTC by Derek Jones
Modified: 2011-12-31 12:48 UTC (History)
2 users (show)

See Also:


Attachments

Note You need to log in before you can comment on or make changes to this bug.
Description Derek Jones 2011-12-30 22:31:18 UTC
File: arithmetic.c
Function: R_pow_di

The existing code:

	if (n < 0) { n = -n; x = 1/x; }
	for(;;) {
	    if(n & 01) xn *= x;
	    if(n >>= 1) x *= x; else break;
	}

should be changed to that proposed below to reduce the error in the result.
For n == 10 an extra digit of accuracy is achieved.

Proposed change:

        is_neg=FALSE;
	if (n < 0) { n = -n; is_neg=TRUE; }
	for(;;) {
	    if(n & 01) xn *= x;
	    if(n >>= 1) x *= x; else break;
	}
        if (is_neg) xn = 1 / xn;
Comment 1 Duncan Murdoch 2011-12-31 00:34:02 UTC
For n==10, the code given should produce identical results to the existing code.

I don't believe it would ever produce "an extra digit" of accuracy.  A claim like that would need support.

Please discuss this on R-devel rather than submitting an unsubstantiated bug report.
Comment 2 Derek Jones 2011-12-31 00:57:46 UTC
(In reply to comment #1)
> For n==10, the code given should produce identical results to the existing
> code.
> 
> I don't believe it would ever produce "an extra digit" of accuracy.  A claim
> like that would need support.
> 
> Please discuss this on R-devel rather than submitting an unsubstantiated bug
> report.

My mistake.  I should have said n == -10

Also see the discussion at:
http://shape-of-code.coding-guidelines.com/2011/12/30/initial-impressions-of-rangelab/
Comment 3 Martin Maechler 2011-12-31 09:21:57 UTC
(In reply to comment #2)
> (In reply to comment #1)
> > For n==10, the code given should produce identical results to the existing
> > code.
> > 
> > I don't believe it would ever produce "an extra digit" of accuracy.  A claim
> > like that would need support.
> > 
> > Please discuss this on R-devel rather than submitting an unsubstantiated bug
> > report.
> 
> My mistake.  I should have said n == -10
> 
> Also see the discussion at:
> http://shape-of-code.coding-guidelines.com/2011/12/30/initial-impressions-of-rangelab/

Hmm, this seems quite convincing to me.
Why can't I reopen it - even if logged into bugzilla?  [Simon]?
Consider it "assigned to me" for now

Martin
Comment 4 Martin Maechler 2011-12-31 12:48:29 UTC
(In reply to comment #3)
> (In reply to comment #2)
> > (In reply to comment #1)
> > > For n==10, the code given should produce identical results to the existing
> > > code.
> > > 
> > > I don't believe it would ever produce "an extra digit" of accuracy.  A claim
> > > like that would need support.
> > > 
> > > Please discuss this on R-devel rather than submitting an unsubstantiated bug
> > > report.
> > 
> > My mistake.  I should have said n == -10
> > 
> > Also see the discussion at:
> > http://shape-of-code.coding-guidelines.com/2011/12/30/initial-impressions-of-rangelab/
> 
> Hmm, this seems quite convincing to me.
> Why can't I reopen it - even if logged into bugzilla?  [Simon]?
> Consider it "assigned to me" for now
> 
> Martin

I have found that  R_pow_di() is not used for R's arithmetic at all, currently
(for accuracy reasons: We've found in the past, and I have confirmed now on one platform, that using C's  pow(x, (double)n)  *is* more accurate in some cases, and never less accurate, even for n as small as 'n == 3'.);
We still use it e.g., inside the code for round(x, digits).

Hence the prioritity now minimal at P1.