more useful options PerlMonks

### Re^7: subtraction issue

by syphilis (Chancellor)
 on Jan 28, 2017 at 06:59 UTC ( #1180519=note: print w/replies, xml ) Need Help??

in reply to Re^6: subtraction issue

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

Replies are listed 'Best First'.
Re^8: subtraction issue
by pryrt (Priest) on Jan 28, 2017 at 18:30 UTC

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.

# easy to come up with one that takes many digits to express

I eventually realized that if I increased the precision to 4000 bits (which is overkill) I would uncover values that require more digits:
```3 / (3566377 ** 47)
0x1.7f1efda69bf79p-1022
3.32997129135359853606999222704654320128433273140834060375717919336765
+633480476074271117603537227797195965402328724231998706388357189066271
+703262139769146889312212556637321188791634569121935603887648273923627
+861919579821592345988060019740748868666186270141425088422184655360375
+823488548150443554012541921413793579099043424106801075148744029407167
+466158885097383679527382985052678222745861830039508807658877054796192
+051679570983765540503801105423327949689320479645855971334458658126010
+779989548774133737798898372623216628041645234343358029716664630265490
+827644527215647502548686221202499548660963335999734105043799288928491
+124484505682527752151459948594916997840488088218758704043473479179413
+158767647452761573740845092915600087271621454476644430542364716529846
+19140625e-308
This(excluding "." and "e-308") is 767 digits.

Cheers,
Rob

Looks lie it addresses what we've been fiddling with.

Cheers,
Rob
Re^8: subtraction issue
by BrowserUk (Pope) on Jan 28, 2017 at 13:45 UTC
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?

With the rise and rise of 'Social' network sites: 'Computers are making people easier to use everyday'
Examine what is said, not who speaks -- Silence betokens consent -- Love the truth but pardon error.
"Science is about questioning the status quo. Questioning authority". The enemy of (IT) success is complexity.
In the absence of evidence, opinion is indistinguishable from prejudice.
I thought that the 'mp' in mpfr stood for multi-precision?

Aaah ... but when a double is assigned to an X-bit (where X > 53) precision mpfr object, the first 53 bits are set the same as the double, and the following bits are all 0.
(To get an X-bit representation of (eg) 1.3, you would need to assign the string '1.3', not the double 1.3.)
A small demo:
```use strict;
use warnings;
use Math::MPFR qw(:mpfr);

my \$prec_53   = Rmpfr_init2(53);
my \$prec_4000 = Rmpfr_init2(4000);

my \$d = 1.3;

Rmpfr_set_d(\$prec_53, \$d, MPFR_RNDN);
# The 53 bits of \$prec_53 are now set to:
# 10100110011001100110011001100110011001100110011001101

Rmpfr_set_d(\$prec_4000, \$d, MPFR_RNDN);
# The first 53 bits are as for \$prec_53
# The remaining 3947 bits are 0.

# But we'll run a check:

die "Things aint right"
if \$prec_53 != \$prec_4000;

print "\$prec_53\n";
# prints 1.3

print "\$prec_4000\n";
# prints 1.3000000000000000444089209850062616169452667236328125
So if I want to assign (from a string) the value held by that double to a 53-bit precision object I can just assign the string "1.3".
But if I want to assign (from a string) the same value to a 4000-bit object then I need to assign the latter string.

I assume that the latter value is an exact representation of the value of the double because after the 53rd decimal digit, all (1153) following decimal digits are 0.

Cheers,
Rob
I assume that the latter value is an exact representation of the value of the double because after the 53rd decimal digit, all (1153) following decimal digits are 0.

Two thoughts:

1. The following extract from the mpfr docs suggests that when mpfr does calculations, it does them to full (infinite) precision and only truncates back to the variable's precision (or the current default) once the calculation is complete:
The semantics of a calculation in MPFR is specified as follows: Compute the requested operation exactly (with “infinite accuracy”), and round the result to the precision of the destination variable, with the given rounding mode. The MPFR floating-point functions are intended to be a smooth extension of the IEEE 754 arithmetic. The results obtained on a given computer are identical to those obtained on a computer with a different word size, or with a different compiler or operating system.

MPFR does not keep track of the accuracy of a computation. This is left to the user or to a higher layer (for example the MPFI library for interval arithmetic). As a consequence, if two variables are used to store only a few significant bits, and their product is stored in a variable with large precision, then MPFR will still compute the result with full precision.

Which I believe means that calculations that end up in 53-bit precision variables may (will often) yield different results from calculations done using IEEE-754 64-bit floats; because of greater accuracy at intermediate points.

2. MPFR seems to produce far more digits of "precision" than it ought.
```In your example above 1 / 526, you posted a result of: 1.9011406844106463 7711435114653113487293012440204620361328125e-3
But the actual result (to 100 digits of precision is : 1.9011406844106463 87832699619771863117870722433460076045627376425855513307984790874524714828897338403e-3
As you can see, only the first 17 digits of the MPFR output are accurate.

As calculated with perl's doubles                    : 1.901140684410646 40000e-003

Ditto your 1 / 535 example: 1.869158878504672 97461032334382480257772840559482574462890625e-3
Actual (to 100 digits)    : 1.869158878504672 897196261682242990654205607476635514018691588785046728971962616822429906542056074766e-3
This time only 16 digits of the MPFR output match.
With perl's doubles       : 1.86915887850467 300000e-003
```

I'm not sure what that all means, but I did validate the 1/526 example manually, as well as calculating them with Math::BigFloat.

With the rise and rise of 'Social' network sites: 'Computers are making people easier to use everyday'
Examine what is said, not who speaks -- Silence betokens consent -- Love the truth but pardon error.
"Science is about questioning the status quo. Questioning authority". The enemy of (IT) success is complexity.
In the absence of evidence, opinion is indistinguishable from prejudice.

Create A New User
Node Status?
node history
Node Type: note [id://1180519]
help
Chatterbox?
and all is quiet...

How do I use this? | Other CB clients
Other Users?
Others musing on the Monastery: (6)
As of 2018-05-21 19:57 GMT
Sections?
Information?
Find Nodes?
Leftovers?
Voting Booth?
World peace can best be achieved by:

Results (161 votes). Check out past polls.

Notices?