Re^8: subtraction issue

[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?

[syphilis]

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

[BrowserUk]

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.

[syphilis]

In your example above 1 / 526, you posted a result of: 1.9011406844106463 7711435114653113487293012440204620361328125e-3

Yes, in this instance (as also in earlier posts) I'm using mpfr to express the exact base 10 equivalent of the value held by the double.
(And that value held by the double is 1/526, calculated to 53 bits of precision.)

But the actual result (to 100 digits of precision is : 1.9011406844106463 87832699619771863117870722433460076045627376425855513307984790874524714828897338403e-3

Yes, this is a much better approximation of 1/526 than that provided by the double, and a different figure is therefore expected.

We can, of course, use mpfr to give us an approximation of 1/526 to whatever precision we want (so long as the precision does not exceed MAX_PREC):

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

my $x = Rmpfr_init2(328);
Rmpfr_set_d($x, 1.0, MPFR_RNDN);
$x /= 526;

print $x;

# Outputs:
# 1.901140684410646387832699619771863117870722433460076045627376425855
+513307984790874524714828897338403e-3
# which is the same as the figure you've provided

Cheers,
Rob