Hello syphilis,

Now, your answer is what I was looking for!

This isn't about finding fault with the current Perl implementation, but about the possibility of improving the future of Perl5.

Were you able to view the referenced article? If not, I'll try to get a pdf of the article for you.

IBM added 54 DFP instructions and the XLC compiler uses them. I haven't looked at 'C' code in more than 20 years, but the idea of getting 34 digits of decimal accuracy when using Perl is down right 'sexy'. What was interesting in the wiki article was that the hardware was 76 times faster than the software for the demo. The wiki is from power6, and power7 was twice the speed of power6, and power8 is more than 3 times the performance of power7. But what is also interesting is that the software implementation is supported by Linux. So if some OpenPower hardware company developed a DFP plug-in, then Perl could have hardware Decimal Floating Point in the future on less expensive hardware. Developers could use the software versions for design/development/testing, and then production could use the hardware.

You gave me a lot to look into, so thanks for your work on this subject.

Regards...Ed

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

Were you able to view the referenced article?

Yes, but the only IBM (I think) processor that I have is a ppc64 on an old mac box - now running Debian Wheezy with gcc-4.6.3.
I think that pre-dates the hardware improvements that you're talking about. Certainly, the _Decimal64/128 operations on that box are generally slower than on my (newer) Intel (Windows) and AMD (Ubuntu) boxes.
My ppc64 processor does, however, utilise IBM's favoured "Densely Packed Decimal" format instead of Intel's/AMD's preferred "Binary Integer Decimal" format.
DPD format is really nice - I can encode/decode that fairly quickly using a hash lookup, whereas BID format necessitates that the value be calculated.
I assume that each of my processors has some special capacity to handle DFP operations, and that gcc is making use of that capability - but there's no guarantee that either of those assumptions is correct.

the idea of getting 34 digits of decimal accuracy when using Perl is down right 'sexy'

I agree, and Math::Decimal128 will give you that - but only for the basic arithmetic operations at the moment.
Perl will improve markedly when someone gets around to (optionally) providing DFP processing built into perl. (But I think the time for that is not yet quite right.)

You've actually now got me wishing I had a power8 box to play with.
Which OS should it run to best test/experience these IBM capabilities ? ... AIX ?
Is power8 bigendian ?

Cheers,
Rob

my $c = 1/3;
[download]

I suppose this will remove some of the problems associated with floating point calculations, but probably not all of them.

Perl 6 has a built-in Rat (rational number) class

It'll be interesting to see how well that performs and what limitations it has. (Are denominator and numerator limited to ivsize ?)
Perl5 does provide access to rational arithmetic via Math::BigRat (which is slow) and Math::GMPq (which is fast, but requires the gmp library).

Cheers,
Rob

Comparing a 2014.08 Rakudo with Perl 5.20 running the rat_mul_div_cancel micro benchmark showed "Rat (rational number) support is 6x faster than Perl 5".

To quote the user doc for Rat:

"To prevent the numerator and denominator from becoming pathologically large, the denominator is limited to 64 bit storage. On overflow of the denominator a Num (floating-point number) is returned instead."

and

"If you want arbitrary precision arithmetic with rational numbers, use the FatRat type instead."

»ö« ^{. o O ( "the celebrity tell-all of the Perl-6 cult?" )}

Hello Laurent_R,

Yes, but it's still in the software. It is one of the features of Perl 6 that I look forward to...Ed

Regards...Ed

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

The only advantage of decimal calculation is that the resulting problems are easier understood because humans have 10 fingers.

But the problems won't go away, eg 1/3 stays periodic in decimal and binary representation.

It's not C ... it's math!

See also Humans have too many fingers

Update

Some languages offer support for rational numbers, i e 1/3 would be represented by two integers 1 and 3 instead of 0.333333...3, but this still keeps the representation of irrational numbers like the root of 2 or pi open.

LanX,

Did you look at the article I quoted?

It has an example of a comparison of Binary Floating Point ( BFP ) and a DFP using both hardware and software to show:

# ./dfp_hw 10 1.000001 60000000
double  fund= 10.0000000000 interest= 1.000000999999999917733362053700
Decimal fund= 10.0000000000 interest= 1.000001000000000000000000000000
Print final funds
  double  fund = 1141973124493563816969240576.0000000000
  Decimal fund = 1141973130130727445029596475.9717600000
[download]

That's a big difference! The difference starts after '11419731'...

Also, my under graduate degree was in math!

Regards...Ed

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

> Did you look at the article I quoted?

No it didn't load, but I never considered leaving my money for 60,000,000 years on a bank account.

Did you read my post about calculating in cents right away? (The one I linked)

Maybe a lexical option for decimal precision could be a reasonable feature request, such that for precision 2 the literal 10.01 is always internally represented as 1001.

Cheers Rolf

PS: Je suis Charlie!

They example you used cheats by using a number that's rational in decimal but irrational in binary. If they had picked a number that's irrational in both (e.g. 1000001/999999 instead of 1000001/1000000), both the BFP and the DFP calculations would suffer from the same problem.

DFP helps in some cases, but it might not be the best solution.

It's not C ... it's math!

Hi LanX, why do you keep shouting at perlmonks?

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]

~/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]

Comments welcome!

Regards...Ed

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

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

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

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

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

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!

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 :-)