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.
###--- 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))
}
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
______________________________________________
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
