Bug #43460 for Math-BigFloat: Bug in Math::BigFloat - gives wrong results for bpow and bexp

Fri Feb 20 07:44:21 2009 nir.nussbaum [...] gmail.com - Ticket created

Subject:

Bug in Math::BigFloat - gives wrong results for bpow and bexp

Hello, I'm working with Perl 5.10 and tried it on Ubuntu and Feodora. bexp gives wrong result for high numbers (negative) and takes ages to calculate, if it is a real negative number such as -394.84010945716266885 (or even -394.8). bpow gave also wrong result for power -394.8 of e. For numbers like -394 it would give the right result. The code: my $i=Math::BigFloat->new(-394.84010945715266885, 20)->bexp(); print ("exp: ".exp(-394.84010945715266885).", BigFloat exp: $i\n"); I'd expect it to give the same results, no? The output: exp: 3.33517962278651e-172, BigFloat exp: 19871330168979164802000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000 0000 Would be happy to hear from you. Best regards, Nir.

Wed Oct 30 18:01:24 2013 DANAJ [...] cpan.org - Correspondence added

Attached is a patch to fix. This is caused by needing far more precision in the intermediate calculations. The easy fix is shown in "Fast Algorithms for High-Precision Computations of Elementary Functions" by Brent, 2006. When x < 0, use exp(x) = 1/exp(-x). For example: perl -Mbignum -E 'say 45.6**-28.4' before: -295.151342321403826190692606470102755169 after: 0.000000000000000000000000000000000000000000000007680658721702948549500271947332317125487 I will point out that that paper also shows methods to speed up the calculation by substantial amounts. A simple demonstration is calculating -8000.8->bexp() with the current code, then with line 1201 in bexp changed from: $x->bpow($x_org, @params); to: $x->bpow(scalar $x_org->bdiv(16,$scale+5),$scale+5)->bpow(16, @params); which runs over 400x faster for me (and passes the test suite, and agrees with PARI/GP). A real patch would be dynamic and also consider things like expm1 or other methods.

Subject:

bigfloat-bpow-fix.patch

--- lib/Math/BigFloat.pm.orig 2013-10-30 13:55:42.426664527 -0700 +++ lib/Math/BigFloat.pm 2013-10-30 13:55:47.262664515 -0700 @@ -2342,6 +2342,15 @@ $HALF = $self->new($HALF) unless ref($HALF); return $x->bsqrt(@r,$y) if $y->bcmp($HALF) == 0; + # if y < 0, use: exp(x) = 1/exp(-x) + if ($y->{sign} eq '-') + { + my $neg_y = $y->copy()->bneg(); + my $xpow = $x->copy->_pow($neg_y, @r); + $x->bone->bdiv($xpow); + return $x; + } + # Using: # a ** x == e ** (x * ln a)

Wed Oct 30 18:01:25 2013 The RT System itself - Status changed from 'new' to 'open'

Wed Oct 30 19:07:47 2013 DANAJ [...] cpan.org - Correspondence added

Here is an alternate patch that fixes the original issue as well as the related issue: perl -Iblib/lib -Iblib/arch -Mbignum -E 'say .01**28.4' before: 648531815674.9002481119466764195627749488 after: 0.000000000000000000000000000000000000000000000000000000001584893192461113485202101373391507013269 .01 ** 23.5: before: -124.4365830026889488973983751515276488572 after: 0.00000000000000000000000000000000000000000000001 The previous patch looked at the sign of y for x^y. This one looks at the sign of u where u = y*ln(x). The should cover both cases.

Subject:

bigfloat-bpow-fix.patch

--- lib/Math/BigFloat.pm.orig 2013-10-30 13:55:42.426664527 -0700 +++ lib/Math/BigFloat.pm 2013-10-30 15:55:08.406646654 -0700 @@ -2389,6 +2389,13 @@ my ($limit,$v,$u,$below,$factor,$next,$over); $u = $x->copy()->blog(undef,$scale)->bmul($y); + # if u < 0, use: exp(x) = 1/exp(-x) + my $do_invert = 0; + if ($u->{sign} eq '-') + { + $do_invert = 1; + $u->bneg(); + } $v = $self->bone(); # 1 $factor = $self->new(2); # 2 $x->bone(); # first term: 1 @@ -2413,6 +2420,8 @@ #$steps++; } + + $x = $x->copy->bone->bdiv($x, $scale) if $do_invert; # shortcut to not run through _find_round_parameters again if (defined $params[0])

Wed Oct 30 20:45:16 2013 DANAJ [...] cpan.org - Correspondence added

It looks like the last line of the last patch needs to be modified to make sure $x is modified in place: if ($do_invert) { my $x_copy = $x->copy; $x->bone->bdiv($x_copy, $scale); } perl -Iblib/lib -Iblib/arch -Mbignum -E 'say 0+(-112.5)->bexp' should print something like: 0.0000000000000000000000000000000000000000000000001386343293641170635024487001056748670217 (perl 5.8.9 through 5.19.5 all print 1650.xxx)

Thu Oct 31 15:56:06 2013 nospam-abuse [...] bloodgate.com - Correspondence added

Subject:	Re: [rt.cpan.org #43460] Bug in Math::BigFloat - gives wrong results for bpow and bexp
Date:	Thu, 31 Oct 2013 15:55:49 -0400
To:	bug-Math-BigFloat [...] rt.cpan.org
From:	"Tels" <nospam-abuse [...] bloodgate.com>

On Wed, October 30, 2013 8:45 pm, Dana Jacobsen via RT wrote: Show quoted text

> Queue: Math-BigFloat > Ticket <URL: https://rt.cpan.org/Ticket/Display.html?id=43460 > > > It looks like the last line of the last patch needs to be modified to make > sure $x is modified in place: > > if ($do_invert) > { > my $x_copy = $x->copy; > $x->bone->bdiv($x_copy, $scale); > }

Might also make sense to add a method to compute 1/x in place like $x->reciprocal()? Regards, Tels PS: Please keep up the good work, it's good to see that the Math modules are still improved. (And I can't say mea culpa oten enough for the problems I left behind :) -- http://bloodgate.com/galleries/

Wed Jan 06 03:59:49 2016 peter.john.acklam [...] gmail.com - Correspondence added

RT-Send-CC:

nospam-abuse [...] bloodgate.com

I have no permission to close this ticket. MARKB or TELS are the only two with permission to do that. This issue fixed in Math-BigInt1.999715.