From: Tels Date: Fri, 22 Jun 2007 19:02:22 +0000 (+0200) Subject: Math::BigInt v1.87 take 10 X-Git-Url: http://git.shadowcat.co.uk/gitweb/gitweb.cgi?a=commitdiff_plain;h=60a1aa196c6751722bae1e1ee83a99d0d965146d;p=p5sagit%2Fp5-mst-13.2.git Math::BigInt v1.87 take 10 Message-Id: <200706221902.22487@bloodgate.com> p4raw-id: //depot/perl@31449 --- diff --git a/lib/Math/BigFloat.pm b/lib/Math/BigFloat.pm index bc4ca90..3c3cf52 100644 --- a/lib/Math/BigFloat.pm +++ b/lib/Math/BigFloat.pm @@ -2244,13 +2244,13 @@ sub bfac sub _pow { - # Calculate a power where $y is a non-integer, like 2 ** 0.5 - my ($x,$y,$a,$p,$r) = @_; + # Calculate a power where $y is a non-integer, like 2 ** 0.3 + my ($x,$y,@r) = @_; my $self = ref($x); # if $y == 0.5, it is sqrt($x) $HALF = $self->new($HALF) unless ref($HALF); - return $x->bsqrt($a,$p,$r,$y) if $y->bcmp($HALF) == 0; + return $x->bsqrt(@r,$y) if $y->bcmp($HALF) == 0; # Using: # a ** x == e ** (x * ln a) @@ -2264,7 +2264,7 @@ sub _pow # we need to limit the accuracy to protect against overflow my $fallback = 0; my ($scale,@params); - ($x,@params) = $x->_find_round_parameters($a,$p,$r); + ($x,@params) = $x->_find_round_parameters(@r); return $x if $x->is_nan(); # error in _find_round_parameters? @@ -2275,7 +2275,7 @@ sub _pow $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; # round mode by caller or undef + $params[2] = $r[2]; # round mode by caller or undef $fallback = 1; # to clear a/p afterwards } else @@ -2380,7 +2380,7 @@ sub bpow } if ($x_is_zero) { - return $x->bone() if $y_is_zero; + #return $x->bone() if $y_is_zero; return $x if $y->{sign} eq '+'; # 0**y => 0 (if not y <= 0) # 0 ** -y => 1 / (0 ** y) => 1 / 0! (1 / 0 => +inf) return $x->binf(); @@ -2522,6 +2522,8 @@ sub bpi { # called like Math::BigFloat::bpi(10); $n = $self; $self = $class; + # called like Math::BigFloat->bpi(); + $n = undef if $n eq 'Math::BigFloat'; } $self = ref($self) if ref($self); $n = 40 if !defined $n || $n < 1; @@ -2570,6 +2572,310 @@ sub bpi $x; } +sub bcos + { + # Calculate a cosinus of x. + my ($self,$x,@r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_); + + # Taylor: x^2 x^4 x^6 x^8 + # cos = 1 - --- + --- - --- + --- ... + # 2! 4! 6! 8! + + # we need to limit the accuracy to protect against overflow + my $fallback = 0; + my ($scale,@params); + ($x,@params) = $x->_find_round_parameters(@r); + + # constant object or error in _find_round_parameters? + return $x if $x->modify('bcos') || $x->is_nan(); + + return $x->bone(@r) if $x->is_zero(); + + # 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 + } + + # when user set globals, they would interfere with our calculation, so + # disable them and later re-enable them + no strict 'refs'; + my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef; + my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef; + # we also need to disable any set A or P on $x (_find_round_parameters took + # them already into account), since these would interfere, too + delete $x->{_a}; delete $x->{_p}; + # need to disable $upgrade in BigInt, to avoid deep recursion + local $Math::BigInt::upgrade = undef; + + my $last = 0; + my $over = $x * $x; # X ^ 2 + my $x2 = $over->copy(); # X ^ 2; difference between terms + my $sign = 1; # start with -= + my $below = $self->new(2); my $factorial = $self->new(3); + $x->bone(); delete $x->{a}; delete $x->{p}; + + my $limit = $self->new("1E-". ($scale-1)); + #my $steps = 0; + while (3 < 5) + { + # we calculate the next term, and add it to the last + # when the next term is below our limit, it won't affect the outcome + # anymore, so we stop: + my $next = $over->copy()->bdiv($below,$scale); + last if $next->bacmp($limit) <= 0; + + if ($sign == 0) + { + $x->badd($next); + } + else + { + $x->bsub($next); + } + $sign = 1-$sign; # alternate + # calculate things for the next term + $over->bmul($x2); # $x*$x + $below->bmul($factorial); $factorial->binc(); # n*(n+1) + $below->bmul($factorial); $factorial->binc(); # n*(n+1) + } + + # shortcut to not run through _find_round_parameters again + if (defined $params[0]) + { + $x->bround($params[0],$params[2]); # then round accordingly + } + else + { + $x->bfround($params[1],$params[2]); # then round accordingly + } + if ($fallback) + { + # clear a/p after round, since user did not request it + delete $x->{_a}; delete $x->{_p}; + } + # restore globals + $$abr = $ab; $$pbr = $pb; + $x; + } + +sub bsin + { + # Calculate a sinus of x. + my ($self,$x,@r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_); + + # taylor: x^3 x^5 x^7 x^9 + # sin = x - --- + --- - --- + --- ... + # 3! 5! 7! 9! + + # we need to limit the accuracy to protect against overflow + my $fallback = 0; + my ($scale,@params); + ($x,@params) = $x->_find_round_parameters(@r); + + # constant object or error in _find_round_parameters? + return $x if $x->modify('bsin') || $x->is_nan(); + + return $x->bzero(@r) if $x->is_zero(); + + # 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 + } + + # when user set globals, they would interfere with our calculation, so + # disable them and later re-enable them + no strict 'refs'; + my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef; + my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef; + # we also need to disable any set A or P on $x (_find_round_parameters took + # them already into account), since these would interfere, too + delete $x->{_a}; delete $x->{_p}; + # need to disable $upgrade in BigInt, to avoid deep recursion + local $Math::BigInt::upgrade = undef; + + my $last = 0; + my $over = $x * $x; # X ^ 2 + my $x2 = $over->copy(); # X ^ 2; difference between terms + $over->bmul($x); # X ^ 3 as starting value + my $sign = 1; # start with -= + my $below = $self->new(6); my $factorial = $self->new(4); + delete $x->{a}; delete $x->{p}; + + my $limit = $self->new("1E-". ($scale-1)); + #my $steps = 0; + while (3 < 5) + { + # we calculate the next term, and add it to the last + # when the next term is below our limit, it won't affect the outcome + # anymore, so we stop: + my $next = $over->copy()->bdiv($below,$scale); + last if $next->bacmp($limit) <= 0; + + if ($sign == 0) + { + $x->badd($next); + } + else + { + $x->bsub($next); + } + $sign = 1-$sign; # alternate + # calculate things for the next term + $over->bmul($x2); # $x*$x + $below->bmul($factorial); $factorial->binc(); # n*(n+1) + $below->bmul($factorial); $factorial->binc(); # n*(n+1) + } + + # shortcut to not run through _find_round_parameters again + if (defined $params[0]) + { + $x->bround($params[0],$params[2]); # then round accordingly + } + else + { + $x->bfround($params[1],$params[2]); # then round accordingly + } + if ($fallback) + { + # clear a/p after round, since user did not request it + delete $x->{_a}; delete $x->{_p}; + } + # restore globals + $$abr = $ab; $$pbr = $pb; + $x; + } + +sub batan + { + # Calculate a arcus tangens of x. + my ($x,@r) = @_; + my $self = ref($x); + + # taylor: x^3 x^5 x^7 x^9 + # 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()!"); + } + + # we need to limit the accuracy to protect against overflow + my $fallback = 0; + my ($scale,@params); + ($x,@params) = $x->_find_round_parameters(@r); + + # constant object or error in _find_round_parameters? + return $x if $x->modify('batan') || $x->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 + } + + # when user set globals, they would interfere with our calculation, so + # disable them and later re-enable them + no strict 'refs'; + my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef; + my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef; + # we also need to disable any set A or P on $x (_find_round_parameters took + # them already into account), since these would interfere, too + delete $x->{_a}; delete $x->{_p}; + # need to disable $upgrade in BigInt, to avoid deep recursion + local $Math::BigInt::upgrade = undef; + + my $last = 0; + my $over = $x * $x; # X ^ 2 + my $x2 = $over->copy(); # X ^ 2; difference between terms + $over->bmul($x); # X ^ 3 as starting value + my $sign = 1; # start with -= + my $below = $self->new(3); + my $two = $self->new(2); + $x->bone(); delete $x->{a}; delete $x->{p}; + + my $limit = $self->new("1E-". ($scale-1)); + #my $steps = 0; + while (3 < 5) + { + # we calculate the next term, and add it to the last + # when the next term is below our limit, it won't affect the outcome + # anymore, so we stop: + my $next = $over->copy()->bdiv($below,$scale); + last if $next->bacmp($limit) <= 0; + + if ($sign == 0) + { + $x->badd($next); + } + else + { + $x->bsub($next); + } + $sign = 1-$sign; # alternate + # calculate things for the next term + $over->bmul($x2); # $x*$x + $below->badd($two); # n += 2 + } + + # shortcut to not run through _find_round_parameters again + if (defined $params[0]) + { + $x->bround($params[0],$params[2]); # then round accordingly + } + else + { + $x->bfround($params[1],$params[2]); # then round accordingly + } + if ($fallback) + { + # clear a/p after round, since user did not request it + delete $x->{_a}; delete $x->{_p}; + } + # restore globals + $$abr = $ab; $$pbr = $pb; + $x; + } + ############################################################################### # rounding functions @@ -3141,6 +3447,11 @@ Math::BigFloat - Arbitrary size floating point math package my $pi = Math::BigFloat->bpi(100); # PI to 100 digits + # the following examples compute their result to 100 digits accuracy: + my $cos = Math::BigFloat->new(1)->bcos(100); # cosinus(1) + my $sin = Math::BigFloat->new(1)->bsin(100); # sinus(1) + my $atan = Math::BigFloat->new(1)->batan(100); # arcus tangens(1) + # Testing $x->is_zero(); # true if arg is +0 $x->is_nan(); # true if arg is NaN @@ -3508,6 +3819,33 @@ Calculate PI to N digits (including the 3 before the dot). This method was added in v1.87 of Math::BigInt (June 2007). +=head2 bcos() + + my $x = Math::BigFloat->new(1); + print $x->bcos(100), "\n"; + +Calculate the cosinus of $x, modifying $x in place. + +This method was added in v1.87 of Math::BigInt (June 2007). + +=head2 bsin() + + my $x = Math::BigFloat->new(1); + print $x->bsin(100), "\n"; + +Calculate the sinus of $x, modifying $x in place. + +This method was added in v1.87 of Math::BigInt (June 2007). + +=head2 batan() + + my $x = Math::BigFloat->new(1); + print $x->batan(100), "\n"; + +Calculate the arcus tanges of $x, modifying $x in place. + +This method was added in v1.87 of Math::BigInt (June 2007). + =head2 bmuladd() $x->bmuladd($y,$z); diff --git a/lib/Math/BigInt.pm b/lib/Math/BigInt.pm index ea4876c..134a4b7 100644 --- a/lib/Math/BigInt.pm +++ b/lib/Math/BigInt.pm @@ -81,9 +81,8 @@ use overload "$_[1]" cmp $_[0]->bstr() : $_[0]->bstr() cmp "$_[1]" }, -# make cos()/sin()/atan2() "work" with BigInt's or subclasses -'cos' => sub { cos($_[0]->numify()) }, -'sin' => sub { sin($_[0]->numify()) }, +'cos' => sub { $_[0]->copy->bcos(); }, +'sin' => sub { $_[0]->copy->bsin(); }, 'atan2' => sub { $_[2] ? atan2($_[1],$_[0]->numify()) : atan2($_[0]->numify(),$_[1]) }, @@ -2927,6 +2926,65 @@ sub bpi $self->new(3); } +sub bcos + { + # Calculate cosinus(x) to N digits. Unless upgrading is in effect, returns the + # result truncated to an integer. + my ($self,$x,@r) = ref($_[0]) ? (undef,@_) : objectify(1,@_); + + return $x if $x->modify('bcos'); + + return $x->bnan() if $x->{sign} !~ /^[+-]\z/; # -inf +inf or NaN => NaN + + return $upgrade->new($x)->bcos(@r) if defined $upgrade; + + # calculate the result and truncate it to integer + my $t = Math::BigFloat->new($x)->bcos(@r)->as_int(); + + $x->bone() if $t->is_one(); + $x->bzero() if $t->is_zero(); + $x->round(@r); + } + +sub bsin + { + # Calculate sinus(x) to N digits. Unless upgrading is in effect, returns the + # result truncated to an integer. + my ($self,$x,@r) = ref($_[0]) ? (undef,@_) : objectify(1,@_); + + return $x if $x->modify('bsin'); + + return $x->bnan() if $x->{sign} !~ /^[+-]\z/; # -inf +inf or NaN => NaN + + return $upgrade->new($x)->bsin(@r) if defined $upgrade; + + # calculate the result and truncate it to integer + my $t = Math::BigFloat->new($x)->bsin(@r)->as_int(); + + $x->bone() if $t->is_one(); + $x->bzero() if $t->is_zero(); + $x->round(@r); + } + +sub batan + { + # Calculate arcus tangens of x to N digits. Unless upgrading is in effect, returns the + # result truncated to an integer. + my ($self,$x,@r) = ref($_[0]) ? (undef,@_) : objectify(1,@_); + + return $x if $x->modify('batan'); + + return $x->bnan() if $x->{sign} !~ /^[+-]\z/; # -inf +inf or NaN => NaN + + return $upgrade->new($x)->batan(@r) if defined $upgrade; + + # calculate the result and truncate it to integer + my $t = Math::BigFloat->new($x)->batan(@r); + + $x->{value} = $CALC->_new( $x->as_int()->bstr() ); + $x->round(@r); + } + ############################################################################### # this method returns 0 if the object can be modified, or 1 if not. # We use a fast constant sub() here, to avoid costly calls. Subclasses @@ -3624,6 +3682,33 @@ before the dot): This method was added in v1.87 of Math::BigInt (June 2007). +=head2 bcos() + + my $x = Math::BigFloat->new(1); + print $x->bcos(100), "\n"; + +Calculate the cosinus of $x, modifying $x in place. + +This method was added in v1.87 of Math::BigInt (June 2007). + +=head2 bsin() + + my $x = Math::BigFloat->new(1); + print $x->bsin(100), "\n"; + +Calculate the sinus of $x, modifying $x in place. + +This method was added in v1.87 of Math::BigInt (June 2007). + +=head2 batan() + + my $x = Math::BigFloat->new(0.5); + print $x->batan(100), "\n"; + +Calculate the arcus tanges of $x, modifying $x in place. + +This method was added in v1.87 of Math::BigInt (June 2007). + =head2 blsft() $x->blsft($y); # left shift in base 2 diff --git a/lib/Math/BigInt/t/bare_mbf.t b/lib/Math/BigInt/t/bare_mbf.t index 7428535..1ff4b52 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 => 2200; + plan tests => 2226; } use Math::BigFloat lib => 'BareCalc'; diff --git a/lib/Math/BigInt/t/bigfltpm.inc b/lib/Math/BigInt/t/bigfltpm.inc index c99cc70..873a4ca 100644 --- a/lib/Math/BigInt/t/bigfltpm.inc +++ b/lib/Math/BigInt/t/bigfltpm.inc @@ -119,6 +119,12 @@ while () $try .= '$x ** $y;'; } elsif ($f eq "bnok") { $try .= '$x->bnok($y);'; + } elsif ($f eq "bcos") { + $try .= '$x->bcos($y);'; + } elsif ($f eq "bsin") { + $try .= '$x->bsin($y);'; + } elsif ($f eq "batan") { + $try .= '$x->batan($y);'; } elsif ($f eq "froot") { $try .= "$setup; \$x->froot(\$y);"; } elsif ($f eq "fadd") { @@ -142,9 +148,9 @@ while () $try .= "\$z = $class->new(\"$args[2]\");"; if( $f eq "bmodpow") { - $try .= "\$x->bmodpow(\$y,\$z);"; + $try .= '$x->bmodpow($y,$z);'; } elsif ($f eq "bmuladd"){ - $try .= "\$x->bmuladd(\$y,\$z);"; + $try .= '$x->bmuladd($y,$z);'; } else { warn "Unknown op '$f'"; } } } @@ -410,6 +416,21 @@ abc:+0:NaN +27:+90:270 +1034:+804:415668 $div_scale = 40; +&bcos +1.2:10:0.3623577545 +2.4:12:-0.737393715541 +0:10:1 +0:20:1 +1:10:0.5403023059 +1:12:0.540302305868 +&bsin +1:10:0.8414709848 +0:10:0 +0:20:0 +2.1:12:0.863209366649 +1.2:13:0.9320390859672 +0.2:13:0.1986693307951 +3.2:12:-0.0583741434276 &bpi 77:3.1415926535897932384626433832795028841971693993751058209749445923078164062862 +0:3.141592653589793238462643383279502884197 diff --git a/lib/Math/BigInt/t/bigfltpm.t b/lib/Math/BigInt/t/bigfltpm.t index 1b9127e..bcaf6e1 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 => 2200 + plan tests => 2226 + 5; # own tests } diff --git a/lib/Math/BigInt/t/fallback.t b/lib/Math/BigInt/t/fallback.t index 010ff8a..c69a3ff 100644 --- a/lib/Math/BigInt/t/fallback.t +++ b/lib/Math/BigInt/t/fallback.t @@ -28,7 +28,7 @@ BEGIN } print "# INC = @INC\n"; - plan tests => 9; + plan tests => 7; } # The tests below test that cos(BigInt) = cos(Scalar) which is DWIM, but not @@ -41,18 +41,14 @@ use Math::BigFloat; my $bi = Math::BigInt->new(1); -ok (cos($bi), cos(1)); -ok (sin($bi), sin(1)); +ok (cos($bi), int(cos(1))); +ok (sin($bi), int(sin(1))); ok (atan2($bi,$bi), atan2(1,1)); my $bf = Math::BigInt->new(0); -ok (cos($bf), cos(0)); -ok (sin($bf), sin(0)); +ok (cos($bf), int(cos(0))); +ok (sin($bf), int(sin(0))); ok (atan2($bi,$bf), atan2(1,0)); ok (atan2($bf,$bi), atan2(0,1)); -my $bone = Math::BigInt->new(1); -ok (cos($bone), cos(1)); -ok (sin($bone), sin(1)); - diff --git a/lib/Math/BigInt/t/sub_mbf.t b/lib/Math/BigInt/t/sub_mbf.t index 08b1ae4..8c15747 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 => 2200 + plan tests => 2226 + 6; # + our own tests } diff --git a/lib/Math/BigInt/t/with_sub.t b/lib/Math/BigInt/t/with_sub.t index a8d3ec8..695fb08 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 => 2200 + plan tests => 2226 + 1; }