summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--lib/Math/BigFloat.pm244
-rw-r--r--lib/Math/BigInt.pm40
-rw-r--r--lib/Math/BigInt/t/bare_mbf.t2
-rw-r--r--lib/Math/BigInt/t/bare_mbi.t2
-rw-r--r--lib/Math/BigInt/t/bigfltpm.inc32
-rwxr-xr-xlib/Math/BigInt/t/bigfltpm.t2
-rw-r--r--lib/Math/BigInt/t/bigintpm.inc23
-rwxr-xr-xlib/Math/BigInt/t/bigintpm.t2
-rwxr-xr-xlib/Math/BigInt/t/sub_mbf.t2
-rwxr-xr-xlib/Math/BigInt/t/sub_mbi.t2
-rw-r--r--lib/Math/BigInt/t/with_sub.t2
11 files changed, 279 insertions, 74 deletions
diff --git a/lib/Math/BigFloat.pm b/lib/Math/BigFloat.pm
index 3762574e9d..7fb9863547 100644
--- a/lib/Math/BigFloat.pm
+++ b/lib/Math/BigFloat.pm
@@ -2488,11 +2488,11 @@ sub _atan_inv
# sign = 1 - sign
my $a = $MBI->_one();
- my $b = $MBI->_new($x);
+ my $b = $MBI->_copy($x);
- my $x2 = $MBI->_mul( $MBI->_new($x), $b); # x2 = x * x
+ my $x2 = $MBI->_mul( $MBI->_copy($x), $b); # x2 = x * x
my $d = $MBI->_new( 3 ); # d = 3
- my $c = $MBI->_mul( $MBI->_new($x), $x2); # c = x ^ 3
+ my $c = $MBI->_mul( $MBI->_copy($x), $x2); # c = x ^ 3
my $two = $MBI->_new( 2 );
# run the first step unconditionally
@@ -2567,7 +2567,7 @@ sub _atan_inv
}
-# print "Took $step steps for $x\n";
+# print "Took $step steps for ", $MBI->_str($x),"\n";
# print "a=",$MBI->_str($a),"\n"; print "b=",$MBI->_str($b),"\n";
# return a/b so that a/b approximates atan(1/x)
($a,$b);
@@ -2597,12 +2597,12 @@ sub bpi
# a few more to prevent rounding errors
$n += 4;
- my ($a,$b) = $self->_atan_inv(239,$n);
- my ($c,$d) = $self->_atan_inv(1023,$n);
- my ($e,$f) = $self->_atan_inv(5832,$n);
- my ($g,$h) = $self->_atan_inv(110443,$n);
- my ($i,$j) = $self->_atan_inv(4841182,$n);
- my ($k,$l) = $self->_atan_inv(6826318,$n);
+ my ($a,$b) = $self->_atan_inv( $MBI->_new(239),$n);
+ my ($c,$d) = $self->_atan_inv( $MBI->_new(1023),$n);
+ my ($e,$f) = $self->_atan_inv( $MBI->_new(5832),$n);
+ my ($g,$h) = $self->_atan_inv( $MBI->_new(110443),$n);
+ my ($i,$j) = $self->_atan_inv( $MBI->_new(4841182),$n);
+ my ($k,$l) = $self->_atan_inv( $MBI->_new(6826318),$n);
$MBI->_mul($a, $MBI->_new(732));
$MBI->_mul($c, $MBI->_new(128));
@@ -2830,38 +2830,139 @@ sub bsin
sub batan2
{
- # calculate arcus tangens of ($x/$y)
+ # calculate arcus tangens of ($y/$x)
# set up parameters
- my ($self,$x,$y,@r) = (ref($_[0]),@_);
+ my ($self,$y,$x,@r) = (ref($_[0]),@_);
# objectify is costly, so avoid it
if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
{
- ($self,$x,$y,@r) = objectify(2,@_);
+ ($self,$y,$x,@r) = objectify(2,@_);
}
- return $x if $x->modify('batan2');
+ return $y if $y->modify('batan2');
- return $x->bnan() if (($x->{sign} eq $nan) ||
- ($y->{sign} eq $nan) ||
- ($x->is_zero() && $y->is_zero()));
+ return $y->bnan() if ($y->{sign} eq $nan) || ($x->{sign} eq $nan);
- # inf handling
- if (($x->{sign} =~ /^[+-]inf$/) || ($y->{sign} =~ /^[+-]inf$/))
+ return $upgrade->new($y)->batan2($upgrade->new($x),@r) if defined $upgrade;
+
+ # Y X
+ # 0 0 result is 0
+ # 0 +x result is 0
+ return $y->bzero(@r) if $y->is_zero() && $x->{sign} eq '+';
+
+ # Y X
+ # 0 -x result is PI
+ if ($y->is_zero())
{
- # XXX TODO:
- return $x->bnan();
+ # calculate PI
+ my $pi = $self->bpi(@r);
+ # modify $x in place
+ $y->{_m} = $pi->{_m};
+ $y->{_e} = $pi->{_e};
+ $y->{_es} = $pi->{_es};
+ $y->{sign} = '+';
+ return $y;
+ }
+
+ # Y X
+ # +y 0 result is PI/2
+ # -y 0 result is -PI/2
+ if ($y->is_inf() || $x->is_zero())
+ {
+ # calculate PI/2
+ my $pi = $self->bpi(@r);
+ # modify $x in place
+ $y->{_m} = $pi->{_m};
+ $y->{_e} = $pi->{_e};
+ $y->{_es} = $pi->{_es};
+ # -y => -PI/2, +y => PI/2
+ $y->{sign} = substr($y->{sign},0,1); # +inf => +
+ $MBI->_div($y->{_m}, $MBI->_new(2));
+ return $y;
}
- return $upgrade->new($x)->batan2($upgrade->new($y),@r) if defined $upgrade;
+ # we need to limit the accuracy to protect against overflow
+ my $fallback = 0;
+ my ($scale,@params);
+ ($y,@params) = $y->_find_round_parameters(@r);
+
+ # error in _find_round_parameters?
+ return $y if $y->is_nan();
+
+ # no rounding at all, so must use fallback
+ if (scalar @params == 0)
+ {
+ # simulate old behaviour
+ $params[0] = $self->div_scale(); # and round to it as accuracy
+ $params[1] = undef; # disable P
+ $scale = $params[0]+4; # at least four more for proper round
+ $params[2] = $r[2]; # round mode by caller or undef
+ $fallback = 1; # to clear a/p afterwards
+ }
+ else
+ {
+ # the 4 below is empirical, and there might be cases where it is not
+ # enough...
+ $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
+ }
+
+ # inlined is_one() && is_one('-')
+ if ($MBI->_is_one($y->{_m}) && $MBI->_is_zero($y->{_e}))
+ {
+ # shortcut: 1 1 result is PI/4
+ # inlined is_one() && is_one('-')
+ if ($MBI->_is_one($x->{_m}) && $MBI->_is_zero($x->{_e}))
+ {
+ # 1,1 => PI/4
+ my $pi_4 = $self->bpi( $scale - 3);
+ # modify $x in place
+ $y->{_m} = $pi_4->{_m};
+ $y->{_e} = $pi_4->{_e};
+ $y->{_es} = $pi_4->{_es};
+ # 1 1 => +
+ # -1 1 => -
+ # 1 -1 => -
+ # -1 -1 => +
+ $y->{sign} = $x->{sign} eq $y->{sign} ? '+' : '-';
+ $MBI->_div($y->{_m}, $MBI->_new(4));
+ return $y;
+ }
+ # shortcut: 1 int(X) result is _atan_inv(X)
+
+ # is integer
+ if ($x->{_es} eq '+')
+ {
+ my $x1 = $MBI->_copy($x->{_m});
+ $MBI->_lsft($x1, $x->{_e},10) unless $MBI->_is_zero($x->{_e});
+
+ my ($a,$b) = $self->_atan_inv($x1, $scale);
+ my $y_sign = $y->{sign};
+ # calculate A/B
+ $y->bone(); $y->{_m} = $a; my $y_d = $self->bone(); $y_d->{_m} = $b;
+ $y->bdiv($y_d, @r);
+ $y->{sign} = $y_sign;
+ return $y;
+ }
+ }
+
+ # handle all other cases
+ # X Y
+ # +x +y 0 to PI/2
+ # -x +y PI/2 to PI
+ # +x -y 0 to -PI/2
+ # -x -y -PI/2 to -PI
+
+ my $y_sign = $y->{sign};
# divide $x by $y
- $x->bdiv($y)->batan(@r);
+ $y->bdiv($x, $scale) unless $x->is_one();
+ $y->batan(@r);
- # set the sign of $x depending on $y
- $x->{sign} = '-' if $y->{sign} eq '-';
+ # restore sign
+ $y->{sign} = $y_sign;
- $x;
+ $y;
}
sub batan
@@ -2874,17 +2975,6 @@ sub batan
# atan = x - --- + --- - --- + --- ...
# 3 5 7 9
- # XXX TODO:
- # This series is only valid if -1 < x < 1, so for other x we need to
- # find a different way.
-
- if ($x < -1 || $x > 1)
- {
- die("$x is out of range for batan()!");
- }
-
- return $x->bzero(@r) if $x->is_zero();
-
# we need to limit the accuracy to protect against overflow
my $fallback = 0;
my ($scale,@params);
@@ -2893,6 +2983,24 @@ sub batan
# constant object or error in _find_round_parameters?
return $x if $x->modify('batan') || $x->is_nan();
+ if ($x->{sign} =~ /^[+-]inf\z/)
+ {
+ # +inf result is PI/2
+ # -inf result is -PI/2
+ # calculate PI/2
+ my $pi = $self->bpi(@r);
+ # modify $x in place
+ $x->{_m} = $pi->{_m};
+ $x->{_e} = $pi->{_e};
+ $x->{_es} = $pi->{_es};
+ # -y => -PI/2, +y => PI/2
+ $x->{sign} = substr($x->{sign},0,1); # +inf => +
+ $MBI->_div($x->{_m}, $MBI->_new(2));
+ return $x;
+ }
+
+ return $x->bzero(@r) if $x->is_zero();
+
# no rounding at all, so must use fallback
if (scalar @params == 0)
{
@@ -2910,6 +3018,35 @@ sub batan
$scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
}
+ # 1 or -1 => PI/4
+ # inlined is_one() && is_one('-')
+ if ($MBI->_is_one($x->{_m}) && $MBI->_is_zero($x->{_e}))
+ {
+ my $pi = $self->bpi($scale - 3);
+ # modify $x in place
+ $x->{_m} = $pi->{_m};
+ $x->{_e} = $pi->{_e};
+ $x->{_es} = $pi->{_es};
+ # leave the sign of $x alone (+1 => +PI/4, -1 => -PI/4)
+ $MBI->_div($x->{_m}, $MBI->_new(4));
+ return $x;
+ }
+
+ # This series is only valid if -1 < x < 1, so for other x we need to
+ # to calculate PI/ - atan(1/x):
+ my $one = $MBI->_new(1);
+ my $pi = undef;
+ if ($x->{_es} eq '+' && ($MBI->_acmp($x->{_m},$one) >= 0))
+ {
+ # calculate PI/2
+ $pi = $self->bpi($scale - 3);
+ $MBI->_div($pi->{_m}, $MBI->_new(2));
+ # calculate 1/$x:
+ my $x_copy = $x->copy();
+ # modify $x in place
+ $x->bone(); $x->bdiv($x_copy,$scale);
+ }
+
# when user set globals, they would interfere with our calculation, so
# disable them and later re-enable them
no strict 'refs';
@@ -2954,6 +3091,17 @@ sub batan
$below->badd($two); # n += 2
}
+ if (defined $pi)
+ {
+ my $x_copy = $x->copy();
+ # modify $x in place
+ $x->{_m} = $pi->{_m};
+ $x->{_e} = $pi->{_e};
+ $x->{_es} = $pi->{_es};
+ # PI/2 - $x
+ $x->bsub($x_copy);
+ }
+
# shortcut to not run through _find_round_parameters again
if (defined $params[0])
{
@@ -3549,6 +3697,10 @@ Math::BigFloat - Arbitrary size floating point math package
my $sin = Math::BigFloat->new(1)->bsin(100); # sinus(1)
my $atan = Math::BigFloat->new(1)->batan(100); # arcus tangens(1)
+ my $atan2 = Math::BigFloat->new( 1 )->batan2( 1 ,100); # batan(1)
+ my $atan2 = Math::BigFloat->new( 1 )->batan2( 8 ,100); # batan(1/8)
+ my $atan2 = Math::BigFloat->new( -2 )->batan2( 1 ,100); # batan(-2)
+
# Testing
$x->is_zero(); # true if arg is +0
$x->is_nan(); # true if arg is NaN
@@ -3935,21 +4087,23 @@ Calculate the sinus of $x, modifying $x in place.
This method was added in v1.87 of Math::BigInt (June 2007).
-=head2 batan()
+=head2 batan2()
- my $x = Math::BigFloat->new(1);
- print $x->batan(100), "\n";
+ my $y = Math::BigFloat->new(2);
+ my $x = Math::BigFloat->new(3);
+ print $y->batan2($x), "\n";
-Calculate the arcus tanges of $x, modifying $x in place.
+Calculate the arcus tanges of C<$y> divided by C<$x>, modifying $y in place.
+See also L<batan()>.
This method was added in v1.87 of Math::BigInt (June 2007).
-=head2 batan2()
+=head2 batan()
- my $x = Math::BigFloat->new(0.5);
- print $x->batan2(0.5), "\n";
+ my $x = Math::BigFloat->new(1);
+ print $x->batan(100), "\n";
-Calculate the arcus tanges of $x / $y, modifying $x in place.
+Calculate the arcus tanges of $x, modifying $x in place. See also L<batan2()>.
This method was added in v1.87 of Math::BigInt (June 2007).
diff --git a/lib/Math/BigInt.pm b/lib/Math/BigInt.pm
index fcb6a786a9..a9bc6c7101 100644
--- a/lib/Math/BigInt.pm
+++ b/lib/Math/BigInt.pm
@@ -2970,33 +2970,30 @@ sub bsin
sub batan2
{
- # calculate arcus tangens of ($x/$y)
+ # calculate arcus tangens of ($y/$x)
# set up parameters
- my ($self,$x,$y,@r) = (ref($_[0]),@_);
+ my ($self,$y,$x,@r) = (ref($_[0]),@_);
# objectify is costly, so avoid it
if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
{
- ($self,$x,$y,@r) = objectify(2,@_);
+ ($self,$y,$x,@r) = objectify(2,@_);
}
- return $x if $x->modify('batan2');
+ return $y if $y->modify('batan2');
- return $x->bnan() if (($x->{sign} eq $nan) ||
- ($y->{sign} eq $nan) ||
- ($x->is_zero() && $y->is_zero()));
+ return $y->bnan() if ($y->{sign} eq $nan) || ($x->{sign} eq $nan);
+
+ return $y->bzero() if $y->is_zero() && $x->{sign} eq '+'; # x >= 0
# inf handling
- if (($x->{sign} =~ /^[+-]inf$/) || ($y->{sign} =~ /^[+-]inf$/))
- {
- # XXX TODO:
- return $x->bnan();
- }
+ # +-inf => --PI/2 => +-1
+ return $y->bone( substr($y->{sign},0,1) ) if $y->{sign} =~ /^[+-]inf$/;
- return $upgrade->new($x)->batan2($upgrade->new($y),@r) if defined $upgrade;
+ return $upgrade->new($y)->batan2($upgrade->new($x),@r) if defined $upgrade;
require Math::BigFloat;
- my $r = Math::BigFloat->new($x)->batan2(Math::BigFloat->new($y),@r)->as_int();
+ my $r = Math::BigFloat->new($y)->batan2(Math::BigFloat->new($x),@r)->as_int();
$x->{value} = $r->{value};
$x->{sign} = $r->{sign};
@@ -3744,24 +3741,25 @@ integer.
This method was added in v1.87 of Math::BigInt (June 2007).
-=head2 batan()
+=head2 batan2()
my $x = Math::BigInt->new(1);
- print $x->batan(100), "\n";
+ my $y = Math::BigInt->new(1);
+ print $y->batan2($x), "\n";
-Calculate the arcus tanges of $x, modifying $x in place.
+Calculate the arcus tangens of C<$y> divided by C<$x>, modifying $y in place.
In BigInt, unless upgrading is in effect, the result is truncated to an
integer.
This method was added in v1.87 of Math::BigInt (June 2007).
-=head2 batan2()
+=head2 batan()
- my $x = Math::BigInt->new(1);
- print $x->batan2(5), "\n";
+ my $x = Math::BigFloat->new(0.5);
+ print $x->batan(100), "\n";
-Calculate the arcus tanges of $x / $y, modifying $x in place.
+Calculate the arcus tangens of $x, modifying $x in place.
In BigInt, unless upgrading is in effect, the result is truncated to an
integer.
diff --git a/lib/Math/BigInt/t/bare_mbf.t b/lib/Math/BigInt/t/bare_mbf.t
index b501695687..53b9fb4ba0 100644
--- a/lib/Math/BigInt/t/bare_mbf.t
+++ b/lib/Math/BigInt/t/bare_mbf.t
@@ -27,7 +27,7 @@ BEGIN
}
print "# INC = @INC\n";
- plan tests => 2236;
+ plan tests => 2292;
}
use Math::BigFloat lib => 'BareCalc';
diff --git a/lib/Math/BigInt/t/bare_mbi.t b/lib/Math/BigInt/t/bare_mbi.t
index 793fd6956c..6b572eee57 100644
--- a/lib/Math/BigInt/t/bare_mbi.t
+++ b/lib/Math/BigInt/t/bare_mbi.t
@@ -26,7 +26,7 @@ BEGIN
}
print "# INC = @INC\n";
- plan tests => 3217;
+ plan tests => 3257;
}
use Math::BigInt lib => 'BareCalc';
diff --git a/lib/Math/BigInt/t/bigfltpm.inc b/lib/Math/BigInt/t/bigfltpm.inc
index 6225146d11..60b5cd2c16 100644
--- a/lib/Math/BigInt/t/bigfltpm.inc
+++ b/lib/Math/BigInt/t/bigfltpm.inc
@@ -434,11 +434,41 @@ $div_scale = 40;
0.2:13:0.1986693307951
3.2:12:-0.0583741434276
&batan
-0.2:14:0.19739555984988
+NaN:10:NaN
+inf:14:1.5707963267949
+-inf:14:-1.5707963267949
0.2:13:0.1973955598499
+0.2:14:0.19739555984988
+0:10:0
+1:14:0.78539816339744
+-1:14:-0.78539816339744
+# test an argument X > 1
+2:14:1.1071487177941
&batan2
+NaN:1:10:NaN
+NaN:NaN:10:NaN
+1:NaN:10:NaN
+inf:1:14:1.5707963267949
+-inf:1:14:-1.5707963267949
1:5:13:0.1973955598499
1:5:14:0.19739555984988
+0:0:10:0
+0:1:14:0
+0:2:14:0
+1:0:14:1.5707963267949
+5:0:14:1.5707963267949
+-1:0:11:-1.5707963268
+-2:0:77:-1.5707963267948966192313216916397514420985846996875529104874722961539082031431
+2:0:77:1.5707963267948966192313216916397514420985846996875529104874722961539082031431
+-1:5:14:-0.19739555984988
+1:5:14:0.19739555984988
+-1:8:14:-0.12435499454676
+1:8:14:0.12435499454676
+-1:1:14:-0.78539816339744
+# test an argument X > 1 and one X < 1
+1:2:24:0.463647609000806116214256
+2:1:14:1.1071487177941
+-2:1:14:-1.1071487177941
&bpi
150:3.14159265358979323846264338327950288419716939937510582097494459230781640628620899862803482534211706798214808651328230664709384460955058223172535940813
77:3.1415926535897932384626433832795028841971693993751058209749445923078164062862
diff --git a/lib/Math/BigInt/t/bigfltpm.t b/lib/Math/BigInt/t/bigfltpm.t
index 83a5da37c6..0956883573 100755
--- a/lib/Math/BigInt/t/bigfltpm.t
+++ b/lib/Math/BigInt/t/bigfltpm.t
@@ -26,7 +26,7 @@ BEGIN
}
print "# INC = @INC\n";
- plan tests => 2236
+ plan tests => 2292
+ 5; # own tests
}
diff --git a/lib/Math/BigInt/t/bigintpm.inc b/lib/Math/BigInt/t/bigintpm.inc
index dbcb469e47..f9481563f6 100644
--- a/lib/Math/BigInt/t/bigintpm.inc
+++ b/lib/Math/BigInt/t/bigintpm.inc
@@ -176,6 +176,8 @@ while (<DATA>)
$try .= "\$x->bmodinv(\$y);";
}elsif ($f eq "digit"){
$try .= "\$x->digit(\$y);";
+ }elsif ($f eq "batan2"){
+ $try .= "\$x->batan2(\$y);";
} else {
# Functions with three arguments
$try .= "\$z = $class->new(\"$args[2]\");";
@@ -2296,6 +2298,27 @@ NaN:NaN
inf:inf
1:2
2:7
+&batan2
+NaN:1:10:NaN
+NaN:NaN:10:NaN
+1:NaN:10:NaN
+inf:1:14:1
+-inf:1:14:-1
+1:5:13:0
+1:5:14:0
+0:0:10:0
+0:1:14:0
+0:2:14:0
+1:0:14:1
+5:0:14:1
+-1:0:11:-1
+-2:0:77:-1
+2:0:77:1
+-1:5:14:0
+1:5:14:0
+-1:8:14:0
+1:8:14:0
+-1:1:14:0
&bpi
77:3
+0:3
diff --git a/lib/Math/BigInt/t/bigintpm.t b/lib/Math/BigInt/t/bigintpm.t
index 86a41c58d6..948f05b097 100755
--- a/lib/Math/BigInt/t/bigintpm.t
+++ b/lib/Math/BigInt/t/bigintpm.t
@@ -10,7 +10,7 @@ BEGIN
my $location = $0; $location =~ s/bigintpm.t//;
unshift @INC, $location; # to locate the testing files
chdir 't' if -d 't';
- plan tests => 3217;
+ plan tests => 3257;
}
use Math::BigInt lib => 'Calc';
diff --git a/lib/Math/BigInt/t/sub_mbf.t b/lib/Math/BigInt/t/sub_mbf.t
index 09a0b29948..6d64cf76f4 100755
--- a/lib/Math/BigInt/t/sub_mbf.t
+++ b/lib/Math/BigInt/t/sub_mbf.t
@@ -26,7 +26,7 @@ BEGIN
}
print "# INC = @INC\n";
- plan tests => 2236
+ plan tests => 2292
+ 6; # + our own tests
}
diff --git a/lib/Math/BigInt/t/sub_mbi.t b/lib/Math/BigInt/t/sub_mbi.t
index b9c598acd0..ea38bffcb8 100755
--- a/lib/Math/BigInt/t/sub_mbi.t
+++ b/lib/Math/BigInt/t/sub_mbi.t
@@ -26,7 +26,7 @@ BEGIN
}
print "# INC = @INC\n";
- plan tests => 3217
+ plan tests => 3257
+ 5; # +5 own tests
}
diff --git a/lib/Math/BigInt/t/with_sub.t b/lib/Math/BigInt/t/with_sub.t
index e6e4815f1c..3e829b1519 100644
--- a/lib/Math/BigInt/t/with_sub.t
+++ b/lib/Math/BigInt/t/with_sub.t
@@ -28,7 +28,7 @@ BEGIN
}
print "# INC = @INC\n";
- plan tests => 2236
+ plan tests => 2292
+ 1;
}
View Lorry Controller Status