BigFloat Precision

Willman023 has asked for the wisdom of the Perl Monks concerning the following question:

Hello,

Question about the precision of BigFloat. In the code below I'm getting a rounding error, or perhaps I'm losing precision even though I am using all BigFloat arithmetic functions and casting all doubles(32 bit) into BigFloats for calculations. To illustrate this I have provided sample output and code. In the first output I get double precision, where as in the second trial I use all BigFloat functions in hopes of getting higher precision, but instead either gross rounding or loss of precision is occuring. Should I not be using the BigFloat functions? If this makes sense to anyone, can anyone tell me what I'm doing wrong, or even offer some ideas. Thanks!

~bW

Sample Output with uncommented line:
c = 10:         0.63331008207472
c = 16:         0.574732628084873
c = 27:         0.507598958744487
c = 46:         0.439468804563385
c = 77:         0.376031147677405
c = 129:        0.316746525404527
c = 359:        0.05 #rounding
c = 599:        1.1  #rounding
c = 1000:       0.04 #rounding

With commented line:
c = 10 pi = 1.6
c = 16 pi = 0.57473262808487273781458803647565021197900281124204984995
+6301387578855115765977011850452313041695
c = 27 pi = 0.5
c = 46 pi = 0.44
c = 77 pi = 0.1
c = 129 pi = 0.32
c = 359 pi = 0.05


#invariant_dist.pl
use strict;
use Math::BigFloat;
use Math::Trig;

my @c = (10, 16, 27, 46, 77, 129, 359, 599, 1000);
my $rho = -5;
my @lambda;
for(0..@c-1) {
    $lambda[$_] = abs($c[$_]-($rho*sqrt($c[$_])));
} 
my $mu = 1;
my @pi = MMCQueueInvDisbn(\@lambda,$mu,\@c);
my $i = 0;
my $sum = 0;
my $mean = 0;

# Compute the invariant distribution for an M/M/c queue
sub MMCQueueInvDisbn { 
    my $lambda = shift;
    my $mu = shift;
    my $c = shift;
    my $rho;
    my $G = Math::BigFloat->new(0);
    my @pi;
    
    # Lets be efficent here!
    my @fact;
    foreach(0..$c[@c-1]) {
        $fact[$_] = 0;
    }
    
    for(my $i=0; $i < @c; $i++) {
        printf($i."\n");
        $rho = $lambda[$i] / $mu;
        $G = MMCQueueNormConst($lambda[$i],$mu,$c[$i],\@fact);
        for (0..$c[$i]-1) {
            $pi[$i][$_] = 0;
        }
        #for n > 250, n! = Inf so $pi = NaN so using BigFloat
        $pi[$i][$c[$i]] = $rho**$c[$i] / (Math::BigFloat->new(Factoria
+l($c[$i], \@fact))->bstr()*$G->bstr());
        #Some Fucking Rounding is occuring
        #$pi[$i][$c[$i]] = (Math::BigFloat->new($rho)->bpow(Math::BigF
+loat->new($c[$i])))->bdiv(Math::BigFloat->new(Factorial($c[$i], \@fac
+t))->bmul($G));
        printf("c = $c[$i]:\t $pi[$i][$c[$i]]\n"); 
    }
    
    return @pi;
}

# Compute the normalization constant for an M/M/c/c queue
sub MMCQueueNormConst { 
    my $lambda = shift;
    my $mu = shift;
    my $c = shift;
    my $fact = shift;
    
    my $rho = $lambda / $mu;
    my $G = Math::BigFloat->new(0);
    for (0..$c) {
        my $num = Math::BigFloat->new(Math::BigFloat->new($rho)->bpow(
+Math::BigFloat->new($_)));
        $G->badd($num->bdiv(Factorial($_, \@$fact)));
    }
    return $G; 
}

# Factorial subroutine
sub Factorial {
    my $n = shift;
    my $fact = shift;
    my $N = Math::BigFloat->new($n);
    
    if ($n == 0) { 
        return 1;
    } elsif($n > 300) {
        # Or could use Sterling's Approximation - trying BigFloat libr
+ary first
        $N->bfac();
        return $N;
    } elsif (@$fact[$n] > 1) {
        return @$fact[$n]; 
    } elsif ($n < 0) { 
        die "$!\n";
    } else {
        $N = $N * Factorial($N-1, \@$fact);
        @$fact[$n] = $N;
        return $N;
    }
}
[download]

Janitored by Arunbear - added readmore tags, as per Monastery guidelines

Comment on BigFloat Precision Download Code

Replies are listed 'Best First'.
Re: BigFloat Precision by TimToady (Parson) on Apr 24, 2005 at 03:00 UTC
If I recall correctly, BigFloat has no way of knowing how much precision you want, which is why it lets you specify it as part of the interface. I see no calls to a precision method in your code, so you're probably getting some kind of default precision.	[reply]
Re^2: BigFloat Precision by belg4mit (Prior) on Apr 26, 2005 at 17:22 UTC
Default precision is "whatever the reult of the math operation is" and hence can grow quite large. Because PDL only offers doubles I tried using BigFloat to fill a Math::Matrix which was then exponentiated to a high power. Without a defined precision the program slowed down for each progressive iteration. Woohoo! `-- I'm not belgian but I play one on TV. On dit que je parle comme un belge aussi.`	[reply]
Re^2: BigFloat Precision by Willman023 (Scribe) on Apr 24, 2005 at 06:47 UTC
Ok, maybe thats what I don't know how to do. I tried the following pre-division: `$numerator->precision(999999999999999999999999999999999999999999999999 +99999999999999999999999999999999999 9999999999999999999999999999999999999999999999999999999999999999999999 +99999999999999999999999999999999999999999999999999999); $denominator->precision(9999999999999999999999999999999999999999999999 +99999999999999999999999999999999999 9999999999999999999999999999999999999999999999999999999999999999999999 +9999999999999999999999999999999999999999999999999999999);` [download] After which the values for numerator and denominator = 0. I tried to set the global precision to this as well, and seemingly everything becomes a NaN. I'm using Perl 5.8 btw, I dunno if that helps. Thanks, bW	[reply] [d/l]
Re: BigFloat Precision by eibwen (Friar) on Apr 23, 2005 at 12:41 UTC
I've noticed that you don't consistently use Math::BigFloat, eg: `$rho = $lambda[$i] / $mu;` There may be an additional problem with your code, but this inconsistent usage is suspect.	[reply] [d/l]
Re^2: BigFloat Precision by belg4mit (Prior) on Apr 26, 2005 at 17:23 UTC
In my expeience this is not a problem however it can be rectified and the code simplified by using the bignum pragma instead. `-- I'm not belgian but I play one on TV. On dit que je parle comme un belge aussi.`	[reply]
Re: BigFloat Precision by gam3 (Curate) on Apr 23, 2005 at 12:54 UTC
It seems to me that your error is happening up at the top of your script in the square root on @c. `for(0..@c-1) { # $lambda[$_] = abs($c[$_]-($rhosqrt($c[$_]))); $lambda[$_] = abs($c[$_]-($rho Math::BigFloat->new($c[$_])->bsqr +t())); }` [download] -- gam3 A picture is worth a thousand words, but takes 200K.	[reply] [d/l]
Re^2: BigFloat Precision by Willman023 (Scribe) on Apr 23, 2005 at 21:22 UTC
Ok gam3 and eibwen, I took out all the regular math you found and replaced with Math::BigFloats, that must have slipped past me. After a few hours of debugging I'm convinced I'm losing precision in the division that I'm doing. The values in the following equation are all BigFloats and the numerator and an even larger denominator are pretty big. When they divide they should produce a very small number, but instead it looks rounded. Perhaps I'm losing the remainder? So it may just be that I don't know how to interpret the documentation on the bdiv() function and how to get the full fraction value, instead of just the quotient. Sample Output: numerator: 2187057715053006176499296408615417037129128897214238256113698344404908 +33433694463191774 1033900096370541032981352897479414902891779671690952835880834484674318 +8026.335842416468 8113979341672817820980901270335818142717085694006013286169841932603829 +32770387188042899 0115508942455532934618389216110784616173476066841503276383837003379903 +94156866919077627 9824783435928019360618178076885689415041227681844670764900853378260631 +99991559227379758 3721140943359293289412998382173652701723265352101625829774791587832397 +09014794285875292 6315552361783352460136562979750515887006847684035410231950364818914502 +59071165536628948 2609717574837339905541901751286406418328318293008678171684392671085100 +45199013199185147 5002542698170815681092942517642292616577517786222417104149554986200865 +05926041669295421 7421195236361801616322976530913722541245775661381844217848774307041052 +85817913486193962 1614955079372880536990393780066744832735297424321264402399647634986599 +22373666181492732 7512892111811467257097624814635258985073512630492058901597404864371649 +85151098402011553 402802193447837576587046088366626941267521 denominator: 1500000000000000000000000000000000000000000000000000000000000000000000 +00000000000000000 0000000000000000000000000000000000000000000000000000000000000000000000 +00000 c = 77 pi = 0.1 #<-- Rounded or something Revised Code: my $numerator = Math::BigFloat->new($rho)->bpow(Math::BigFloat->new($c +[$i])); my $denominator = Math::BigFloat->new(Factorial($c[$i], \@fact))->bmul +($G); printf("numerator: %s\n", $numerator ); printf("denominator: %s\n", $denominator ); $pi[$i][$c[$i]] = $numerator->copy()->bdiv($denominator); #$pi[$i][$c[$i]] = (Math::BigFloat->new($rho)->bpow(Math::BigFloat->ne +w($c[$i])))->bdiv(Math::BigFloat->new(Factorial($c[$i], \@fact))->bmu +l($G) ); printf("c = $c[$i] pi = $pi[$i][$c[$i]]\n"); [download] Thanx ~bW	[reply] [d/l]
Re: BigFloat Precision by gam3 (Curate) on Apr 24, 2005 at 19:17 UTC
Try this code. Math::BigFloat->new('2187057715053006176499296408615417037129128897214 +238256113698344404908334336944631917741033900096370541032981352897479 +4149028917796716909528358808344846743188026.3358424164688113979341672 +817820980901270335818142717085694006013286169841932603829327703871880 +428990115508942455532934618389216110784616173476066841503276383837003 +379903941568669190776279824783435928019360618178076885689415041227681 +844670764900853378260631999915592273797583721140943359293289412998382 +173652701723265352101625829774791587832397090147942858752926315552361 +783352460136562979750515887006847684035410231950364818914502590711655 +366289482609717574837339905541901751286406418328318293008678171684392 +671085100451990131991851475002542698170815681092942517642292616577517 +786222417104149554986200865059260416692954217421195236361801616322976 +530913722541245775661381844217848774307041052858179134861939621614955 +079372880536990393780066744832735297424321264402399647634986599223736 +661814927327512892111811467257097624814635258985073512630492058901597 +404864371649851510984020115534028021934478375765870460883666269412675 +21'); my $denominator = Math::BigFloat->new('1000000000000000000000000000000 +000000000000000000000000000000000000000000000000000000000000000000000 +0000000000000000000000000000000000000000000000000000000000000'); printf("numerator: %s\n", $numerator + 1 ); printf("numerator: %s\n", $numerator ); printf("denominator: %s\n", $denominator ); my $x = $numerator->copy(); $x->accuracy(10); print scalar $x->bdiv($denominator), "\n"; $x = $numerator->copy(); $x->accuracy(1000); print scalar $x->bdiv($denominator), "\n"; [download] -- gam3 A picture is worth a thousand words, but takes 200K.	[reply] [d/l]
Re^2: BigFloat Precision by Willman023 (Scribe) on Apr 24, 2005 at 20:26 UTC
Ok I ran your suggestion, with the only change that the top line is my $numerator = Math:BigFloat->new('etc..'); The following is the output from that: numerator: 21870577150530061764992964086154170371291288972142382561136 +983444049083343369446319177410339000963705410329813528974794149028917 +796716909528358808344846743188027.33584241646881139793416728178209809 0127033581814271708569400601328616984193260382932770387188042899011550 +894245553293461838921611078461617347606684150327638383700337990394156 +8669190776279824783435928019360618178076885689415041227681 8446707649008533782606319999155922737975837211409433592932894129983821 +736527017232653521016258297747915878323970901479428587529263155523617 +8335246013656297975051588700684768403541023195036481891450 2590711655366289482609717574837339905541901751286406418328318293008678 +171684392671085100451990131991851475002542698170815681092942517642292 +6165775177862224171041495549862008650592604166929542174211 9523636180161632297653091372254124577566138184421784877430704105285817 +913486193962161495507937288053699039378006674483273529742432126440239 +9647634986599223736661814927327512892111811467257097624814 6352589850735126304920589015974048643716498515109840201155340280219344 +7837576587046088366626941267521 numerator: 21870577150530061764992964086154170371291288972142382561136 +983444049083343369446319177410339000963705410329813528974794149028917 +796716909528358808344846743188026.33584241646881139793416728178209809 0127033581814271708569400601328616984193260382932770387188042899011550 +894245553293461838921611078461617347606684150327638383700337990394156 +8669190776279824783435928019360618178076885689415041227681 8446707649008533782606319999155922737975837211409433592932894129983821 +736527017232653521016258297747915878323970901479428587529263155523617 +8335246013656297975051588700684768403541023195036481891450 2590711655366289482609717574837339905541901751286406418328318293008678 +171684392671085100451990131991851475002542698170815681092942517642292 +6165775177862224171041495549862008650592604166929542174211 9523636180161632297653091372254124577566138184421784877430704105285817 +913486193962161495507937288053699039378006674483273529742432126440239 +9647634986599223736661814927327512892111811467257097624814 6352589850735126304920589015974048643716498515109840201155340280219344 +7837576587046088366626941267521 denominator: 100000000000000000000000000000000000000000000000000000000 +000000000000000000000000000000000000000000000000000000000000000000000 +00000000000000000000000000000000000 2.187057715 2.18705771505300617649929640861541703712912889721423825611369834440490 +833433694463191774103390009637054103298135289747941490289177967169095 +283588083448467431880263358424164688113979341672817820980901270335818 +142717085694006013286169841932603829327703871880428990115508942455532 +934618389216110784616173476066841503276383837003379903941568669190776 +279824783435928019360618178076885689415041227681844670764900853378260 +631999915592273797583721140943359293289412998382173652701723265352101 +625829774791587832397090147942858752926315552361783352460136562979750 +515887006847684035410231950364818914502590711655366289482609717574837 +339905541901751286406418328318293008678171684392671085100451990131991 +851475002542698170815681092942517642292616577517786222417104149554986 +200865059260416692954217421195236361801616322976530913722541245775661 +381844217848774307041052858179134861939621614955079372880536990393780 +066744832735297424321264402399647634986599223736661814927327512892111 +8114672570976248146352589850735126 [download] The following is my code with my numerator and denominator whereby I attempt the same series of divisions and accuracy changes, (Also tried increasing the accuracy even more): my $numerator = Math::BigFloat->new($rho)->bpow(Math::BigFloat->new($c +[$i])); my $denominator = Math::BigFloat->new(Factorial($c[$i], \@fact))->bmul +($G); printf("numerator: %s\n", $numerator ); printf("denominator: %s\n", $denominator ); my $x = $numerator->copy(); $x->accuracy(10); print scalar $x->bdiv($denominator), "\n"; $x = $numerator->copy(); $x->accuracy(1000); print scalar $x->bdiv($denominator), "\n"; exit(); Output: numerator: 131254331066182.0671680983636563623093384645074900631038196 +916724632735210716611428572612616761979022584315623777773238891425361 +932674318752893801 denominator: 200000000000000 1.6 #<--- Way wrong should be around .65 [download] Thanks, bW	[reply] [d/l] [select]
Re: BigFloat Precision by belg4mit (Prior) on Apr 26, 2005 at 17:26 UTC
Take a look at my comments deeper down. Otherwise all I can say is good luck, I've had similar problems and was unable to solve them satisfactorily. I did notice though that unless you need precision exceeding 30+ decimal places that perl scalars (in 5.8 at least) will give the same results and are much much faster. PS> I wonder if the problem is not so much precision as accuracy ;-) `-- I'm not belgian but I play one on TV. On dit que je parle comme un belge aussi.`	[reply]
Re^2: BigFloat Precision by Anonymous Monk on Apr 26, 2005 at 17:48 UTC
I did notice though that unless you need precision exceeding 30+ decimal places that perl scalars (in 5.8 at least) will give the same results and are much much faster. You'll only get 30+ decimal digits if you're using some special hardware and a special version of perl to take advantage of it (Sparc-V9 hardware has 128bit quad floating point numbers). If you're using IEEE754 64bit doubles (like on Intel hardware) you'll only get about 15 significant decimal digits. `perl -e 'printf "%.30f\n",(1/3)'`	[reply] [d/l]