[perl #81400] Fix bmodinv() part of RT 63237

The following standard definition is used: z is the modular inverse of
x (mod y) if and only if x*z (mod y) = 1 (mod y).

- dist/Math-BigInt/lib/Math/BigInt.pm: Fix the code in bmodinv() so it can
  handle negative arguments. The code can be optimized further for speed,
  but correctnes first.

- dist/Math-BigInt/t/bigintpm.inc: Fix the test case for modinv(-3, -5).
  The output should be -3, since -3 * -3 (mod -5) = -9 (mod -5) = -4, and
  1 (mod -5) = -4.

- dist/Math-BigRat/t/bigratpm.inc: Fix same test case as above.
  Math::BigRat::bmodinv() only handles integers, and is essentially just a
  front-end to Math::BigInt::bmodinv().

3 files changed, 74 insertions, 27 deletions

diff --git a/dist/Math-BigInt/lib/Math/BigInt.pm b/dist/Math-BigInt/lib/Math/BigInt.pm
index 3e9196a985..52acc7a682 100644
--- a/dist/Math-BigInt/lib/Math/BigInt.pm
+++ b/dist/Math-BigInt/lib/Math/BigInt.pm
@@ -1785,10 +1785,12 @@ sub bmod
 
 sub bmodinv
   {
-  # Modular inverse.  given a number which is (hopefully) relatively
-  # prime to the modulus, calculate its inverse using Euclid's
-  # algorithm.  If the number is not relatively prime to the modulus
-  # (i.e. their gcd is not one) then NaN is returned.
+  # Return modular multiplicative inverse: z is the modular inverse of x (mod
+  # y) if and only if x*z (mod y) = 1 (mod y). If the modulus y is larger than
+  # one, x and z are relative primes (i.e., their greatest common divisor is
+  # one).
+  #
+  # If no modular multiplicative inverse exists, NaN is returned.
 
   # set up parameters
   my ($self,$x,$y,@r) = (undef,@_);
@@ -1800,22 +1802,58 @@ sub bmodinv
 
   return $x if $x->modify('bmodinv');
 
-  return $x->bnan()
-        if ($y->{sign} ne '+'                           # -, NaN, +inf, -inf
-         || $x->is_zero()                               # or num == 0
-         || $x->{sign} !~ /^[+-]$/                      # or num NaN, inf, -inf
-        );
-
-  # put least residue into $x if $x was negative, and thus make it positive
-  $x->bmod($y) if $x->{sign} eq '-';
-
-  my $sign;
-  ($x->{value},$sign) = $CALC->_modinv($x->{value},$y->{value});
-  return $x->bnan() if !defined $x->{value};		# in case no GCD found
-  return $x if !defined $sign;			# already real result
-  $x->{sign} = $sign;				# flip/flop see below
-  $x->bmod($y);					# calc real result
-  $x;
+  # Return NaN if one or both arguments is +inf, -inf, or nan.
+
+  return $x->bnan() if ($y->{sign} !~ /^[+-]$/ ||
+                        $x->{sign} !~ /^[+-]$/);
+
+  # Return NaN if $y is zero; 1 % 0 makes no sense.
+
+  return $x->bnan() if $y->is_zero();
+
+  # Return 0 in the trivial case. $x % 1 or $x % -1 is zero for all finite
+  # integers $x.
+
+  return $x->bzero() if ($y->is_one() ||
+                         $y->is_one('-'));
+
+  # Return NaN if $x = 0, or $x modulo $y is zero. The only valid case when
+  # $x = 0 is when $y = 1 or $y = -1, but that was covered above.
+  #
+  # Note that computing $x modulo $y here affects the value we'll feed to
+  # $CALC->_modinv() below when $x and $y have opposite signs. E.g., if $x =
+  # 5 and $y = 7, those two values are fed to _modinv(), but if $x = -5 and
+  # $y = 7, the values fed to _modinv() are $x = 2 (= -5 % 7) and $y = 7.
+  # The value if $x is affected only when $x and $y have opposite signs.
+
+  $x->bmod($y);
+  return $x->bnan() if $x->is_zero();
+
+  # Compute the modular multiplicative inverse of the absolute values. We'll
+  # correct for the signs of $x and $y later. Return NaN if no GCD is found.
+
+  ($x->{value}, $x->{sign}) = $CALC->_modinv($x->{value}, $y->{value});
+  return $x->bnan() if !defined $x->{value};
+
+  # When one or both arguments are negative, we have the following
+  # relations.  If x and y are positive:
+  #
+  #   modinv(-x, -y) = -modinv(x, y)
+  #   modinv(-x,  y) = y - modinv(x, y)  = -modinv(x, y) (mod y)
+  #   modinv( x, -y) = modinv(x, y) - y  =  modinv(x, y) (mod -y)
+
+  # We must swap the sign of the result if the original $x is negative.
+  # However, we must compensate for ignoring the signs when computing the
+  # inverse modulo. The net effect is that we must swap the sign of the
+  # result if $y is negative.
+
+  $x -> bneg() if $y->{sign} eq '-';
+
+  # Compute $x modulo $y again after correcting the sign.
+
+  $x -> bmod($y) if $x->{sign} ne $y->{sign};
+
+  return $x;
   }
 
 sub bmodpow
@@ -3175,7 +3213,7 @@ Math::BigInt - Arbitrary size integer/float math package
 
   $x->bmod($y);		   # modulus (x % y)
   $x->bmodpow($exp,$mod);  # modular exponentation (($num**$exp) % $mod))
-  $x->bmodinv($mod);	   # the inverse of $x in the given modulus $mod
+  $x->bmodinv($mod);	   # the multiplicative inverse of $x modulo $mod
 
   $x->bpow($y);		   # power of arguments (x ** y)
   $x->blsft($y);	   # left shift in base 2
@@ -3694,11 +3732,20 @@ This method was added in v1.87 of Math::BigInt (June 2007).
 
 =head2 bmodinv()
 
-	num->bmodinv($mod);		# modular inverse
+	$x->bmodinv($mod);		# modular multiplicative inverse
+
+Returns the multiplicative inverse of C<$x> modulo C<$mod>. If
+
+        $y = $x -> copy() -> bmodinv($mod)
+
+then C<$y> is the number closest to zero, and with the same sign as C<$mod>,
+satisfying
+
+        ($x * $y) % $mod = 1 % $mod
 
-Returns the inverse of C<$num> in the given modulus C<$mod>.  'C<NaN>' is
-returned unless C<$num> is relatively prime to C<$mod>, i.e. unless
-C<bgcd($num, $mod)==1>.
+If C<$x> and C<$y> are non-zero, they must be relative primes, i.e.,
+C<bgcd($y, $mod)==1>. 'C<NaN>' is returned when no modular multiplicative
+inverse exists.
 
 =head2 bmodpow()
 
diff --git a/dist/Math-BigInt/t/bigintpm.inc b/dist/Math-BigInt/t/bigintpm.inc
index e8460041a6..822b44610e 100644
--- a/dist/Math-BigInt/t/bigintpm.inc
+++ b/dist/Math-BigInt/t/bigintpm.inc
@@ -1704,13 +1704,13 @@ abc:5:NaN
 # bmodinv Expected Results from normal use
 1:5:1
 3:5:2
+3:-5:-3
 -2:5:2
 8:5033:4404
 1234567891:13:6
 -1234567891:13:7
 324958749843759385732954874325984357439658735983745:2348249874968739:1741662881064902
 ## bmodinv Error cases / useless use of function
-3:-5:NaN
 inf:5:NaN
 5:inf:NaN
 -inf:5:NaN
diff --git a/dist/Math-BigRat/t/bigratpm.inc b/dist/Math-BigRat/t/bigratpm.inc
index 13ec6240d8..de5e9e1687 100644
--- a/dist/Math-BigRat/t/bigratpm.inc
+++ b/dist/Math-BigRat/t/bigratpm.inc
@@ -185,13 +185,13 @@ abc:5:NaN
 # bmodinv Expected Results from normal use
 1:5:1
 3:5:2
+3:-5:-3
 -2:5:2
 8:5033:4404
 1234567891:13:6
 -1234567891:13:7
 324958749843759385732954874325984357439658735983745:2348249874968739:1741662881064902
 ## bmodinv Error cases / useless use of function
-3:-5:NaN
 inf:5:NaN
 5:inf:NaN
 -inf:5:NaN