Re^5: subtraction issue

[syphilis]

it takes n fixed-point decimal digits to exactly represent an n-bit fixed-point fractional binary number

Thank you for correcting me ... but shouldn't that be "n+1 fixed-point decimal digits" ?

use strict;
use warnings;
use Math::MPFR qw(:mpfr);

Rmpfr_set_default_prec(1000);

my $d = 1 / 11;
my $x = Math::MPFR->new($d);

print "$x\n", length($x), "\n";

# (print() removes trailing zeros.)
# Outputs:
9.09090909090909116141432377844466827809810638427734375e-2
58

Take out the "." and the "e-2" and we're left with 54 decimal digits.

Cheers,
Rob

[pryrt]

That's confusing to me, because the math says that for a fraction that's exactly representable with d digits after the binary point, f = m * 2**-d = m * 5**d / (5**d * 2**d) = m * 5**d / 10**d, which is exactly representable with d digits after the decimal point, because m * 5**d is an integer, and an integer over 10**d only needs d digits after the decimal point: 17 / 10**5 = 0.00017.

I will try to find some time this weekend to see what's going on. 1/11 isn't exactly representable in binary, but assuming that mpfr is using the exact value of the double precision, I would think it would fit in the 52 digits after the point. Oh, wait! The 52 is actually the digits of the fraction of the mantissa after removing enough powers of 2 to get the mantissa between 1 and 2. So internally, it would be 1.4545... * 2**-4. So we've got 52 binary digits after the double floating point, plus four more digits due the power of two, which makes d = 58: the 54 SIG figs you show, plus the 2 zeroes from e-2, plus apparently 2 zeroes after the final 5 shown. (This would have been so much easier if I were at a computer with perl, not just my tablet and my emulator in my skull...)

[syphilis]

... assuming that mpfr is using the exact value of the double precision ...

It is - unless there's a bug at work.

which makes d = 58

Here's a couple of values that, according to mpfr, require 60 decimal digits:

1 / 526
1.90114068441064637711435114653113487293012440204620361328125e-3

1 / 535
1.86915887850467297461032334382480257772840559482574462890625e-3

Might be time I stopped fiddling and started thinking. (Have you got a spare one of those in-skull emulators that you could send over ?)

Update: Fiddling is much easier, especially on a Sat'dy night. I've found numerous doubles that apparently require 122 decimal digits for exact representation, but haven't found any that require more than 122. For example:

25 / (8360484 ** 15)
3.66833261843476555514935702258199343058161821865425601769198458144384
+41059633175331427349000780680559003454762069283595612e-103

Cheers,
Rob

[pryrt]

I can actually find one a lot longer, using the structure. Start with ou = 1 + ULP => 1 + 2**-52 -- that's the farthest separation between two active bits possible in a double. Then, shift it as far to the right of the decimal point as you can: x = ou * 2**-1022 = 2**-1022 + 2**-1074. The MPFR shows 303 digits. But I don't think that's all the digits. I actually expect 2**-1022 to take 1022 digits (including leading zeros, not including the "0.") after a fixed decimal point. And 2**-1074 should take 1074 digits (including leading zeros) after a fixed decimal point, so the sum of those two should really go to the 10**-1074th place. Since it starts at the 308th digit after the fixed decmial point, I actually expect 765 digits.

use warnings;
use strict;
use Data::IEEE754::Tools 0.016 qw/:all/;
use Math::MPFR qw/:mpfr/;
use POSIX qw/floor ceil log10 log2/;

Rmpfr_set_default_prec(1000);

# easy to come up with one that takes many digits to express:
my $x = nextUp( POS_NORM_SMALLEST() );
print "POS_NORM_SMALLEST = 2**-1022\n";
print "POS_NORM_SMALLEST+1ULP = (1 + 2**-52)*(2**-1022)\n";
print "\n";
print "(      1  )           \n";
print "(1 + -----)*(2**-1022)\n";
print "(    2**52)           \n";
print "\n";
printf "IEEE754 Hex String: '%s'\n", hexstr754_from_double($x);
printf "%-30a %-30s %-30s\n", $x, to_hex_floatingpoint($x), to_dec_flo
+atingpoint($x);
my $s=Math::MPFR->new($x);
printf "[%d digits long]: %s\n", length($s)-6, $s;

my $d = ceil(-log10($x));   # start digit
printf "For _FIXED POINT_ notation,\n";
printf "\tstarts %d digits after the decimal point\n", $d;
my $ulp = ulp($x);          # value of last bit (whether or not it is 
+set)
printf "long version of ulp: %s\n", Math::MPFR->new($ulp);
my $l2 = log2($ulp);
printf "log2(ulp) = %d, so I expect it to fit within %d digits of the 
+fixed decimal point\n", $l2, -$l2;

__OUTPUT__
POS_NORM_SMALLEST = 2**-1022
POS_NORM_SMALLEST+1ULP = (1 + 2**-52)*(2**-1022)

(      1  )
(1 + -----)*(2**-1022)
(    2**52)

IEEE754 Hex String: '0010000000000001'
0x1.0000000000001p-1022        +0x1.0000000000001p-1022       +0d1.000
+0000000000002p-1022
[303 digits long]: 2.2250738585072018771558785585789482407880088486837
+041956131300312119688603996006965297904292212628858639037013670281908
+017171296072711910355127227413175152199055740043138804567803233377539
+881639177387328959246074229270113078053813397081653361296447449529789
+5212189790907838525833659018517896187998851504e-308
For _FIXED POINT_ notation,
        starts 308 digits after the decimal point

Since I don't know MPFR, I tried to look at it using Math::BigFloat, but I couldn't get it to stop prematurely rounding at about 30 digits after the fixed point, even with accuracy(1000) and precision(-1000). Apparently, I don't know BigFloat very well, either. I'm pretty busy the rest of the weekend, but I know what my coffee breaks at work next week are going to include: writing my own simplistic string-based fixed-point bigs, that can simply add two bigs, or divide a big by two, without all the overhead of modules I don't understand that seem to be rounding prematurely.

[BrowserUk]

ssuming that mpfr is using the exact value of the double precision ... It is

I thought that the 'mp' in mpfr stood for multi-precision?