mpf_cmp

[Toine Bommelijn]

Hi Paul and Torbjörn,

I disagree that the numbers are slightly different in value, the calculation
is: truncate((1/12 * 1) * 100) / 100
First part is: (1/12 * 1) * 100) = 8.3333333333333333333333333 .... this is
a number that cannot be exactly represented with a binary representation.
I then use the following function to truncate the 8.3333... number, this is
the source of my truncate function:

<start>
            //Truncate is a tricky operation -> GNUmp is much more accurate
in
            //its calculations but floats still have an error, if the result
            //of a calculation is "124.999999999999999999999999999999999993"
            //is the result of truncate "124" and not 125 as expected, to
            //prevent this error -> do the following
            //1. Find the nearest integer,  = NI
            //2. Calculate absolute difference, Diff = abs(NI - float)
            //3. Is this difference smaller than the number-of-DIGITS
precision
            //   set by the user then its a rounding error, truncate(
float ) = NI

            //Round to the nearest upper integer
            if( mpf_sgn(Top->_GFloat) > 0 ) {
                mpf_ceil(_TmpGFloat, Top->_GFloat);
                }
            else {
                mpf_floor(_TmpGFloat, Top->_GFloat);
                }

            //Calculate the absolute difference
            mpf_sub(_TmpGFloat, _TmpGFloat, Top->_GFloat);

            if( mpf_sgn(_TmpGFloat) < 0 ) {
                mpf_neg(_TmpGFloat, _TmpGFloat);
                }

            //Calculate the smallest possible with the current number of
digits
            mpf_t Error;
            mpf_init( Error );
            CalcSmallestValue( _NumDig, Error );

            if( mpf_cmp(_TmpGFloat, Error) < 0 ) {
                //The difference between the two numbers is smaller than the
                //current precision in digits -> so its a rounding error
                //
                //Again calculate the nearest integer to store it as a value
                if( mpf_sgn(Top->_GFloat) > 0 ) {
                    mpf_ceil(_TmpGFloat, Top->_GFloat);
                    }
                else {
                    mpf_floor(_TmpGFloat, Top->_GFloat);
                    }
                mpz_set_f(_TmpGInt, _TmpGFloat);
                }
            else {
                //Truncate the float
                mpf_trunc(_TmpGFloat, Top->_GFloat);
                //Convert the value to an integer
                mpz_set_f(_TmpGInt, _TmpGFloat);
                }

            mpf_clear( Error );
            Top->PutGIntValue(_TmpGInt);
<end>

Then I divide the INTEGER value 8 by 100 -> the result is 0.8, I think this
is a floating-point value that can be exactly represented by a
floating-point representation like GNUmp uses, its not a rational number
like 0.83333333333...

Just before the comparison I added a call to function "gmp_asprintf", it
says that both arguments are 0.08 although their internal representation is
different, then I expect the "mpf_cmp" function to work in a similar way
(unless it was documented as behaving differently).
I traced through the source of "mpf_cmp" and the call to "mpn_cmp" inside
"mpf_cmp" returns 0 -> meaning identical values, I still find it strange
that "mpf_cmp" returns the value 1

<start - part of "mpf_cmp">
  if (usize > vsize)
    {
      cmp = mpn_cmp (up + usize - vsize, vp, vsize);
      if (cmp == 0)
 return usign;
    }
<end - part of "mpf_cmp">

In the coming Trinc-Prolog there is the possibility to switch run-time
between normal 64 bit float calculations and using the GNUmp DLL. I
performed the same calculations with normal 64 bit arithmetic and it
succeeded where GNUmp fails. That is another argument in favor of incorrect
"mpf_cmp" behavior.

Regards,
Toine Bommelijn

ps: thanks for the quick response in this discussion