[Rd] print(...,digits=2) behavior

[Petr Savicky]

On Mon, Feb 07, 2011 at 09:56:18AM +0100, Martin Maechler wrote:
[...]
> Namely, one important purpose of the test is to ensure that e.g.
> 
>   print(3.597, digits = 3)
> 
> is printed as  3.6   and not 3.60
> 
> Now I have had -- since 1997 at least -- an R version of
> scientific() for more easy experimentation.
> Here's an updated version of that:

[...]
> 
> Now, I'm waiting for comments/suggestions ..

In a previous email, i discussed the case print(9.991, digits=3).
This case is already handled in your R level experimental function
scientific(). I noticed this only after sending the previous email.

I would like to suggest a modification of the algorithm. The purpose
is to determine the minimum number of digits, such that the printed number
has the same value as the number rounded to exactly getOption("digits")
digits, but shorter representation.

Algorithm. First, the mantissa of the number rounded to "digits" digits
is computed as an integer by the formula

  a <- floor(alpha*10^(digits - 1) + 0.5)

The obtained integer number a satisfies

  10^(digits-1) <= a <= 10^digits

The case a = 10^digits corresponds to the situation that we have to
increase the exponent. This may be tested either before the remaining
part of the code or as the last step as below.

For example, if alpha = 3.597 and digits = 3, we get a = 360. The
question, whether (digits - 1) digits are sufficient for printing the
number, is equivalent to the condition that "a" is divisible by 10.
Similarly, (digits - 2) digits are sufficient if and only if "a" is
divisible by 100. This may be tested using the following code

    nsig <- digits
    for (j in 1:nsig) {
        a <- a / 10
        if (a == floor(a)) {
            nsig <- nsig - 1
        } else {
            break
        }
    }
    ## nsig == 0 if and only if we started with a == 10^digits
    if (nsig == 0) {
        nsig <- 1
        kpower <- kpower + 1
    }

This code uses double precision for a, since values up to 10^digits
may occur and digits may be 15 or more. The algorithm is not exact
for digits = 15, but works reasonably. I suggest to use a different,
slower and more accurate algorithm, if getOption("digits") >= 16.

Please, find in an attachment two versions of the above algorithm. One
is a modification of your R level function from scientific.R and the
other is a patch against src/main/format.c, which i tested.

I suggest to consider the following test cases. The presented output
is obtained using the patch.

  example1 <- c(
  7.94999999999999,
  7.95000000000001,
  8.04999999999999,
  8.05000000000001)
  for (x in example1) print(x, digits=2)
  
  [1] 7.9
  [1] 8
  [1] 8
  [1] 8.1
  
  example2 <- c(
  3.59949999999999,
  3.59950000000001,
  3.60049999999999,
  3.60050000000001)
  for (x in example2) print(x, digits=4)
  
  [1] 3.599
  [1] 3.6
  [1] 3.6
  [1] 3.601
  
  example3 <- c(
  1.00000049999999,
  1.00000050000001)
  for (x in example3) print(x, digits=7)
  
  [1] 1
  [1] 1.000001
  
  example4 <- c(
  9.99999949999999,
  9.99999950000001)
  for (x in example4) print(x, digits=7)
  
  [1] 9.999999
  [1] 10

I appreciate comments.

Petr Savicky.

-------------- next part --------------
###--- R function that does the same as 'scientific' in

###--- /usr/local/R-0.50/src/main/format.c
###--- ~~~~~~~~~~~||||||~~~~~~~~~~~~~~~~~~

scientific1 <- function(x, digits = getOption('digits')) ## eps = 10^-(digits)
{
  ##-   /* for 1 number  x , return
  ##-    *      sgn    = 1_{x < 0}  {0/1}
  ##-    *      kpower = Exponent of 10;
  ##-    *      nsig   = min(digits, #{significant digits of alpha})
  ##-    *
  ##-    * where  |x| = alpha * 10^kpower   and  1 <= alpha < 10
  ##-    */

  eps <- 10 ^(-digits)

  if (x == 0) {
    kpower <- 0
    nsig <- 1
    sgn <- 0
  } else {
    if(x < 0) {
      sgn <- 1; r <- -x
    } else {
      sgn <- 0; r <- x
    }
    kpower <- floor(log10(r));##-->      r = |x| ;  10^k <= r

    if (kpower <= -308) ## close to denormalization -- added for R 2.x.y
        alpha <- (r * 1e+30) / 10^(kpower+30)
    else
        alpha <- r / 10^kpower

    ## "a" integer, 10^(digits-1) <= a <= 10^digits

    a <- floor(alpha*10^(digits - 1) + 0.5)
    nsig <- digits
    for (j in 1:nsig) {
        a <- a / 10
        if (a == floor(a)) {
            nsig <- nsig - 1
        } else {
            break
        }
    }
    if (nsig == 0) {
        nsig <- 1
        kpower <- kpower + 1
    }
  }
  left <- kpower + 1
  c(sgn = sgn, kpower = kpower, nsig = nsig,
    left = left, right = nsig - left,
    sleft = sgn + max(1,left))
}

-------------- next part --------------
diff --minimal -U 5 -r R-devel/src/main/format.c R-devel-print/src/main/format.c
--- R-devel/src/main/format.c	2010-04-27 17:52:24.000000000 +0200
+++ R-devel-print/src/main/format.c	2011-02-12 11:07:37.000000000 +0100
@@ -167,28 +167,26 @@
 	    alpha = (r * 1e+30)/pow(10.0, (double)(kp+30));
 	}
 	else
 	    alpha = r / pow(10.0, (double)kp);
 
-	/* make sure that alpha is in [1,10) AFTER rounding */
-
-	if (10.0 - alpha < eps*alpha) {
-	    alpha /= 10.0;
-	    kp += 1;
-	}
-	*kpower = kp;
-
-	/* compute number of digits */
-
-	*nsig = R_print.digits;
-	for (j = 1; j <= *nsig; j++) {
-	    if (fabs(alpha - floor(alpha+0.5)) < eps * alpha) {
-		*nsig = j;
-		break;
-	    }
-	    alpha *= 10.0;
-	}
+        /* alpha integer, 10^(digits-1) <= alpha <= 10^digits */
+        alpha = floor(alpha*pow(10.0, R_print.digits - 1.0) + 0.5);
+        *nsig = R_print.digits;
+        for (j = 1; j <= R_print.digits; j++) {
+            alpha /= 10.0;
+            if (alpha == floor(alpha)) {
+                (*nsig)--;
+            } else {
+                break;
+            }
+        }
+        if (*nsig == 0) {
+            *nsig = 1;
+            kp += 1;
+        }
+        *kpower = kp;
     }
 }
 
 /*
    The return values are