in reply to Decimal Floating Point (DFP) and does Perl needs DFP?

Dear Monks,

Just an update on this thread. I ran the following testcase on a Debian system w/o hardware DFP. I then did a Perl one-liner, that I believe does the same thing. I needed some printf statements to help figure out what I was doing, so I left them in to help any one else. 

/* ------------------------------------------------------------------ 
    To compile:  gcc -o dfpal2 dfpal2.c -L . -lDFPAL
   ------------------------------------------------------------------ 
+*/

#include <stdio.h>
#include <stdlib.h>
#include "dfpal.h"

int main(int32_t argc, char **argv)
{
    int     i, count;
    char    mystring[100];
    char    sn1[100];
    char    sn2[100];
    int32_t init_err, init_os_err;
    char    *err_str = NULL;
    dfpalflag_t st;
    decimal128 n1, n2;                    /* Declare the decimal floin
+ting point numbers */

    printf("Start:\n");
    if (dfpalInit((void *)malloc(dfpalMemSize()))     /* start the lib
+rary */
                       != DFPAL_ERR_NO_ERROR) {
        fprintf(stderr,"DFPAL init error\n");
        printf("Error:\n");
        dfpalEnd(free);
        return (1);
    }
    n1 = dec128FromString(argv[1]);
    n2 = dec128FromString(argv[2]);            /* assign a value by co
+nverting the input string to a demical */
    count = 60000000;
    printf("Initialized: %s \t %s\n",dec128ToString(n1,sn1),dec128ToSt
+ring(n2,sn2) );
    printf("for loop:\n");
    for (i = 0; i < count; i++)
    {
//          printf("Step %d: %s \t %s\n",i,dec128ToString(n1,sn1),dec1
+28ToString(n2,sn2) );
        n1 = dec128Multiply(n1, n2);            /* do the maths */
    //    printf("Step %d: %d \t %d\n",i,i,i );
     }

    printf("Final amount=%s\n",dec128ToString(n1,mystring));    /* pri
+nt the results */

    dfpalEnd(free);                    /* end the library cleanly */
}

[download]
The results were better than I expected. YMMV. 
~/dfpal# time ./dfpal2 10. 1.000001
Start:
Initialized:   10        1.000001
for loop:
Final amount=1141973130130727445029596475.971760

real    0m10.722s
user    0m10.725s
sys     0m0.000s

~/dfpal# time perl -e '$n1=10;$n2=1.000001;for ($i=0;$i<60_000_000;$i+
++) { $n1=$n1*$n2;} print "Final amount=$n1\t$n2\n";'

Final amount=1.14197312449358e+27       1.000001

real    0m8.698s
user    0m8.701s
sys     0m0.000s

~/dfpal#

[download]
I had expected that the software only solution would be several times slower than the binary floating point hardware. Big Surprise!

Comments welcome!

Regards...Ed

"Well done is better than well said." - Benjamin Franklin

Replies are listed 'Best First'.
Re^2: Decimal Floating Point (DFP) and does Perl needs DFP?
by syphilis (Archbishop) on Jan 19, 2015 at 01:52 UTC 
    real    0m10.722s
user    0m10.725s
sys     0m0.000s

    [download]


    With the Math::Decimal128 one-liner, on my Windows 7 box, perl-5.20.0: 
    C:\>perl -MMath::Decimal128=":all" -le "$t = time();$n1=Math::Decimal1
+28->new('10', 0);$n2=Math::Decimal128->new('1000001', -6);for(1..6000
+0000){$n1 *= $n2}print $n1;$t2 = time() - $t;print $t2;"
114197313013072744502959647597176e-5
25

    [download]
    25 seconds is not great, but OTOH it's not too bad for 60 million overloaded XSub calls.

    Cheers,
    Rob
[reply]
[d/l]
[select]

      Hello syphilis ,

      I got 35-37 seconds on a Debian Linux with Perl 5.14.2.

      Regards...Ed

      "Well done is better than well said." - Benjamin Franklin

[reply]
Re^2: Decimal Floating Point (DFP) and does Perl needs DFP?
by Anonymous Monk on Jan 18, 2015 at 16:31 UTC

    I tried two extended double versions of given loop: __float80 (hardware) and __float128 (software emulation). The results were, respectively: 

    1141973130132393144424071168.000000
1141973130132392916052606976.000000
    Difference in speed -- about 60-fold.

    ps. would you mind stating the mathematically correct result (say, to six digits after decimal point)?
[reply]
      would you mind stating the mathematically correct result (say, to six digits after decimal point)?

      According to Math::Decimal128: 
      C:\>perl -MMath::Decimal128=":all" -le "$n1=Math::Decimal128->new('10'
+, 0);$n2=Math::Decimal128->new('1000001', -6);for(1..60000000){$n1 *=
+ $n2}print $n1;"
114197313013072744502959647597176e-5

      [download]
      Or, to 6 decimal places: 1141973130130727445029596475.971760 which is the same as flexvault got.

      With 80-bit long doubles, the calculation came out as: 
      C:\>perl -MMath::LongDouble=":all" -le "$n1=Math::LongDouble->new('10'
+);$n2=Math::LongDouble->new('1.000001');for(1..60000000){$n1 *= $n2}p
+rint $n1;"
1.14197313013239314e+027

      [download]
      Or, printing to 28 decimal digits (which implies a level of precision that we don't have): 
      C:\>perl -MMath::LongDouble=":all" -le "$n1=Math::LongDouble->new('10'
+);$n2=Math::LongDouble->new('1.000001');for(1..60000000){$n1 *= $n2}p
+rint LDtoSTRP($n1, 28);"
1.141973130132393144424071168e+027

      [download]
      This figure is equivalent to the first one you quoted.

      With __float128: 
      C:\>perl -MMath::Float128=":all" -le "$n1=Math::Float128->new('10');$n
+2=Math::Float128->new('1.000001');for(1..60000000){$n1 *= $n2}print $
+n1;"
1.14197313013072744502959647306684e+27

      [download]
      This differs from your second figure - and I don't think there's any valid reason that they should differ.
      I'm inclined to think that mine is correct (surprise :-) because it matches up pretty well with the _Decimal128 calculation (where the respective precisions are reasonably close to each other).
      Might you have assigned a long double or double value that got cast to a __float128, instead of assigning an actual __float128 value ?

      As LanX points out, rounding is involved. The above one-liner programs all employ rounding to nearest with ties to even - the different results of course arising from the different precisions involved.

      Cheers,
      Rob
[reply]
[d/l]
[select]
        > rounding is involved.

        Calculating 10 * 1.000001 ** 6e7 should minimize the propagation of rounding errors, if exponentiation is sufficiently well implemented.

        Cheers Rolf

        PS: Je suis Charlie!

[reply]

        Right you are about the __float128 case. After rewriting the constant as (__float128)1000001 / 1e6, my result now matches your calculation.

        I'm inclined to think the bc-calculated value is accurate. But verifying that is rather laborious, for 1000001**6e7 is a number with ca 360M decimal digits.

[reply]
[d/l]
      > would you mind stating the mathematically correct result (say, to six digits after decimal point)?

      This should be correct for at least 50 digits after decimal point:

       
1141973130130727445029596475.97165373568561208619839490929389898529956
+04918087368059784750017004880305063884032742217053556637130

      [download]

      calculated using Math::BigFloat (which operates decimal with configurable precision) and taking advantage of exponentiation laws described here.

      (NB: BigFloat's own ->bpow() method is far slower)

      use strict;
use warnings;
use Math::BigFloat;

my $exp=6e7;
Math::BigFloat->precision(-100);
my $x=Math::BigFloat->new('1.000001');
my $val=Math::BigFloat->new('1');
$\="\n";


for my $bit (reverse split //,sprintf '%b',$exp) {
    # print "$bit ",$x;
    $val->bmul($x) if $bit;
    $x->bpow(2);
}

$val->bmul(10);

print "Result : $val";

my $syphilis=Math::BigFloat->new('1141973130130727445029596475.971760'
+);

print "diff: ", $val - $syphilis;

      [download]

      took about a second on a netbook.

      Compilation started at Tue Jan 20 02:19:30

/usr/bin/perl -w /home/lanx/pm/big_expo.pl 
Result : 1141973130130727445029596475.97165373568561208619839490929389
+89852995604918087368059784750017004880305063884032742217053556637130
diff: -0.0001062643143879138016050907061010147004395081912631940215249
+982995119694936115967257782946443362870

Compilation finished at Tue Jan 20 02:19:31

      [download]

      The results from the C-libs are as you can see only accurate till the third digit after decimal point.

      I'm not a big expert on BigFloat, if you see a problem, all feedback welcome. :)

      Cheers Rolf

      PS: Je suis Charlie!

      update

      you can easily increase the accuracy to 1000 digits after decimal point by setting precision to -1050, only takes a second.

[reply]
[d/l]
[select]
        This should be correct for at least 50 digits after decimal point:
        1141973130130727445029596475.9716537356856120861983949092938989852995604918087368059784750017004880305063884032742217053556637130

        Yes, it's correct for 69 digits after the decimal point. The following, provided by Math::MPFR, appears to be correct for 129 digits - as increasing the mpfr precision above 450 bits doesn't alter the 129 digit value:
        11419731301307274450295964759.716537356856120861983949092938989852995604918087368059784750017004844856171227209226348596180640102

        And that agrees (to at least 129 digits) with LanX's script when M::BF precision is set to -1050.

        Cheers,
        Rob
[reply]
      Someone please correct me, but taken that each operation causes up to one bit of error and that we perform 60e6 operations ...

      We'll need calculate with 60.000.000 bits to be 100% sure of avoiding rounding errors.

      So a brute force calculation should be ruled out.

      Furthermore what is mathematical correct?

      Though it's taught in school, any rounding that maps 0.5 to 1 is with 50% probability wrong.

      Cheers Rolf

      PS: Je suis Charlie!

[reply]
      Update.

      In above, s/52606976/something else/ — a cast to __float80 was used for the lack of proper printf format. But the error is magnitudes higher than that, so... Regarding the expected result: 

      $ echo 'scale=80; 10 * (1.000001^78125)^768' | bc
1141973130130727445029596475.971653735685612
      Looks like the 60M figure was carefully picked :-)

[reply]

In Section Seekers of Perl Wisdom