=head1 NAME

Catacomb::MP - Multiprecision integer arithmetic

=head1 SYNOPSIS

    use Catacomb qw(:const :mp :random :pgen);

    $x = Catacomb::MP->new($str, [$radix]);
    $x = Catacomb::MP->new($i);
    $x = Catacomb::MP->loadb($bytes);
    $x = Catacomb::MP->loadl($bytes);
    $x = Catacomb::MP->loadb2c($bytes);
    $x = Catacomb::MP->loadl2c($bytes);
    ($x, $rest) = Catacomb::MP->fromstring($str, [$radix]);
    $x = mp($str, [$radix]);
    $x = mp($i);
    ($x, $rest) = mp_fromstring($str, [$radix]);
    $a = $b + $c;
    $a = $b - $c;
    $a = $b * $c;
    $a = $b / $c;
    $a = $b % $c;
    $a = $b ** $n;
    $a = $b << $n;
    $a = $b >> $n;
    $a = $b & $c;
    $a = $b | $c;
    $a = $b ^ $c;
    $a = ~$b;
    $p = $b == $c;
    $p = $b != $c;
    $p = $b < $c;
    $p = $b > $c;
    $p = $b <= $c;
    $p = $b >= $c;
    $a = sqrt($b);
    $a = -$b;
    $a = $b->add($c);
    $a = $b->sub($c);
    $a = $b->mul($c);
    ($q, $r) = $a->div($b);
    $a = $b->sqr();
    $a = $b->exp($c);
    $a = $b->sqrt();
    $a = $b->neg();
    $a = $b->not();
    $a = $b->not2c();
    $a = $b->mod($c);
    $p = $b->eq($c);
    $cmp = $b->cmp($c);
    $a = $b->and($c);
    $a = $b->and2c($c);
    $a = $b->or($c);
    $a = $b->or2c($c);
    $a = $b->xor($c);
    $a = $b->xor2c($c);
    $a = $b->nand($c);
    $a = $b->nand2c($c);
    $a = $b->nor($c);
    $a = $b->nor2c($c);
    $a = $b->not();
    $a = $b->not2c();
    $a = $b->lsl($n);
    $a = $b->lsl2c($n);
    $a = $b->lsr($n);
    $a = $b->lsr2c($n);
    $a = $b->setbit($n);
    $a = $b->setbit2c($n);
    $a = $b->clearbit($n);
    $a = $b->clearbit2c($n);
    $p = $b->testbit($n);
    $p = $b->testbit2c($n);
    $x = $y->copy(); # largely useless
    $g = $a->gcd($b);
    ($g, $u, $v) = $a->gcd($b);
    ($s, $t) = $m->odd(); # m = 2^s t
    $a = $p->modexp($x, $e);
    $a = $p->modinv($b);
    $r = $p->modsqrt($x);
    $q = $n->jac($a);
    $p = $x->primep([$rng]);
    $nbits = $x->bits();
    $nbytes = $x->octets();
    $bytes = $x->storeb([$nbytes]);
    $bytes = $x->storel([$nbytes]);
    $bytes = $x->storeb2c([$nbytes]);
    $bytes = $x->storel2c([$nbytes]);
    $str = $x->tostring([$radix]);
    $n = $x->toint();

    $barrett = Catacomb::MP::Barrett->new($p);
    $barrett = $p->barrett();
    $p = $barrett->m();
    $x = $barrett->reduce($y);
    $a = $barrett->exp($x, $y);

    $mont = Catacomb::MP::Mont->new($p);
    $mont = $p->mont();
    $r = $mont->r();
    $r2 = $mont->r2();
    $p = $mont->m();
    $x = $mont->in($y);
    $a = $mont->mul($x, $y);
    $a = $mont->expr($x, $y);
    $a = $mont->mexpr($x0, $e0, $x1, $e1, ...);
    $x = $mont->reduce($y);
    $x = $mont->out($y);
    $a = $mont->exp($x, $y);
    $a = $mont->mexp($x0, $e0, $x1, $e1, ...);

    $reduce = Catacomb::MP::Reduce->new($p);
    $reduce = $p->mkreduce();
    $p = $reduce->m();
    $x = $reduce->reduce($y);
    $a = $barrett->exp($x, $y);

    $crt = Catacomb::MP::CRT->new(@n);
    $n = $crt->product();
    @n = $crt->moduli();
    $x = $crt->solve(@r);

=head1 DESCRIPTION

This is a full-featured multiprecision maths library.  Integer objects
belong to the class C<Catacomb::MP>.  Most standard arithmetic operators
are overloaded and do the Right Thing.  See below for more details.

Note that multiprecision integers are I<immutable>.  No operation can
change the value of an integer.

=head2 Constructors

=over

=item B<Catacomb::MP-E<gt>new(>i<int>B<)>

=item B<mp(>I<int>B<)>

Return the multiprecision integer whose value is I<int>.

=item B<Catacomb::MP-E<gt>new(>I<string>, [I<radix>]B<)>

=item B<mp(>I<string>, [I<radix>]B<)>

Returns the multiprecision integer whose value is represented by
I<string>, in the given I<radix>.  If I<radix> is 0, as it is by
default, then I<string> may have a prefix giving its radix: it may be
B<0> for octal, B<0x> for hex, or I<r>B<_> for radix I<r>.  If the
argument is invalid, B<undef> is returned.

=item B<Catacomb::MP-E<gt>loadb(>I<bytes>B<)>

=item B<Catacomb::MP-E<gt>loadl(>I<bytes>B<)>

=item B<mp_loadb(>I<bytes>B<)>

=item B<mp_loadl(>I<bytes>B<)>

Returns the multiprecision integer whose value is represented in big- or
little-endian base-256 form by I<bytes>.  This will always be
nonnegative.

Returns the multiprecision integer whose value is represented in
little-endian base-256 form by I<bytes>.  This will always be
nonnegative.

=item B<Catacomb::MP-E<gt>loadb2c(>I<bytes>B<)>

=item B<Catacomb::MP-E<gt>loadl2c(>I<bytes>B<)>

=item B<mp_loadb2c(>I<bytes>B<)>

=item B<mp_loadl2c(>I<bytes>B<)>

Returns the multiprecision integer whose value is represented in
two's complement big- or little-endian base-256 form by I<bytes>.

=item B<Catacomb::MP-E<gt>fromstring(>I<string>, [I<radix>]B<)>

=item B<mp_fromstring(>I<string>, [I<radix>]B<)>

Returns the multiprecision integer whose value is represented by
I<string>, and (in a list context) the remainder of I<string>.  If
I<radix> is 0 then I<string> may have a prefix, as for B<new()>.  If
I<string> is invalid, an empty list is returned.

=back

=head2 Methods

The B<Catacomb::MP> package overloads the standard arithmetic
operators.  As long as one operand is a B<Catacomb::MP> object, the
other may be any argument acceptable to B<new()>.  As a special
exception, either argument to B<*> may be a B<Catacomb::EC::Pt> object,
to perform elliptic curve point multiplication.

=over

=item I<a>B<-E<gt>add(>I<b>B<)>

=item I<a>B<-E<gt>sub(>I<b>B<)>

=item I<a>B<-E<gt>mul(>I<b>B<)>

=item I<a>B<-E<gt>div(>I<b>B<)>

=item I<a>B<-E<gt>mod(>I<b>B<)>

The basic algorithms.  The argument I<b> may be any value acceptable to
the B<new()> constructor.

In an array context, the B<div()> method returns a two values: the
quotient and remainder.  The quotient returned is the I<floor> of the
true quotient: the remainder therefore has the same sign as the divisor.
This is usually the right result for number-theoretic purposes.

=item I<a>B<-E<gt>sqr()>

Returns the square of its argument.  Usually faster than multiplying it
by itself.

=item I<a>B<-E<gt>exp(>I<b>B<)>

Returns the result I<a>^I<b>.  Theoretically, I<b> can be a
multiprecision integer, but the result will be impossibly huge if this
is the case.  The exponent I<b> may not be negative.

=item I<a>B<-E<gt>sqrt()>

Returns the integer square-root of I<a> -- i.e., the largest integer
less than or equal to the true square root.

=item I<a>B<-E<gt>eq(>I<b>B<)>

Returns true if I<a> and I<b> are equal; or false otherwise.

=item I<a>B<-E<gt>cmp(>I<b>B<)>

Returns -1, 0 or +1 according to whether I<a> is less than, equal to or
greater than I<b>.

=item I<a>B<-E<gt>and(>I<b>B<)>

=item I<a>B<-E<gt>and2c(>I<b>B<)>

=item I<a>B<-E<gt>or(>I<b>B<)>

=item I<a>B<-E<gt>or2c(>I<b>B<)>

=item I<a>B<-E<gt>xor(>I<b>B<)>

=item I<a>B<-E<gt>xor2c(>I<b>B<)>

=item I<a>B<-E<gt>nand(>I<b>B<)>

=item I<a>B<-E<gt>nand2c(>I<b>B<)>

=item I<a>B<-E<gt>nor(>I<b>B<)>

=item I<a>B<-E<gt>nor2c(>I<b>B<)>

=item I<a>B<-E<gt>not()>

=item I<a>B<-E<gt>not2c()>

Standard bitwise operators.  The unadorned methods ignore the sign of
their operands and use only the absolute values.  The B<2c> versions use
a two's complement representation, with the sign bit repeated
infinitely.

=item I<a>B<-E<gt>lsl(>I<n>B<)>

=item I<a>B<-E<gt>lsl2c(>I<n>B<)>

=item I<a>B<-E<gt>lsr(>I<n>B<)>

=item I<a>B<-E<gt>lsr2c(>I<n>B<)>

Standard shifting operators.  The unadorned methods treat only absolute
values, while the B<2c> versions use a two's complement representation,
with the sign bit repeated infinitely.

=item I<a>B<-E<gt>setbit(>I<n>B<)>

=item I<a>B<-E<gt>setbit2c(>I<n>B<)>

=item I<a>B<-E<gt>clearbit(>I<n>B<)>

=item I<a>B<-E<gt>clearbit2c(>I<n>B<)>

Returns a copy of the argument I<a> with bit I<n> set or clear.  The
unadorned methods treat only absolute values, while the B<2c> variants
use a two's complement representation with the sign bit repeated
infinitely.

=item I<a>B<-E<gt>testbit(>I<n>B<)>

=item I<a>B<-E<gt>testbit2c(>I<n>B<)>

Returns true or false according to whether I<a> has bit I<n> set or
clear.  The unadorned methods treat only absolute values, while the
B<2c> variant uses a two's complement representation with the sign bit
repeated infinitely.

=item I<a>B<-E<gt>copy()>

Returns a copy of I<a>.  Not useful: use assignment instead.

=item I<a>B<-E<gt>gcd(>I<b>B<)>

Returns the greatest common divisor of I<a> and I<b>.  This will never
be negative.  In a list context, it also returns two further results,
I<u> and I<v>: if the GCD is I<g> then we have I<a> I<u> + I<b> I<v> =
I<g>.  Furthermore, I<v> will have the same sign as I<b>.

=item I<a>B<-E<gt>odd()>

Returns a pair of numbers I<s> and I<t>, so that I<a> = I<t> << I<s>,
with I<t> odd and I<s> as large as possible.  I<t> will be
multiprecision; I<s> will be a standard Perl integer.

=item I<a>B<-E<gt>modexp(>I<b>B<,> I<n>B<)>

Returns I<b>^I<n> mod I<a>.

=item I<a>B<-E<gt>modinv(>I<b>B<)>

Returns the inverse of I<b> modulo I<a>.  If no inverse exists then
report an error.

=item I<a>B<-E<gt>modsqrt(>I<b>B<)>

If I<a> is prime, and I<b> is a quadratic residue mod I<a>, then return
a square root of I<b> mod I<a>.  If I<a> is not prime, does something
arbitrary; if I<b> is a nonresidue then report an error.

=item I<a>B<-E<gt>jac(>I<b>B<)>

Returns the Jacobi symbol (I<b>/I<a>).  If I<a> is prime, then this is
I<b>^((I<a> - 1)/2), i.e., 0 if I<b> is a multiple of I<a>, +1 if I<b>
is a quadratic residue modulo I<a>, or -1 if I<b> is a quadratic
nonresidue mod I<a>.

=item I<a>B<-E<gt>primep(>[I<rng>]I<)>

Returns true or false, depending on whether I<a> is (probably) prime.
If I<rng> is specified, it must be a random number generator -- see
Catacomb::Crypto(3).

=item I<a>B<-E<gt>bits()>

Returns the smallest number of bits required to represent the integer
I<a>.

=item I<a>B<-E<gt>octets()>

Returns the smallest number of bytes required to represent the integer
I<a>, ignoring sign.

=item I<a>B<-E<gt>octets2c()>

Returns the smallest number of bytes required to represent the integer
I<a> as two's complement.  This may be one greater than the result
returned by B<octets()> in the case where the top bit of I<a> is the top
bit of an octet.

=item I<a>B<-E<gt>storeb(>[I<nbytes>]B<)>

=item I<a>B<-E<gt>storel(>[I<nbytes>]B<)>

Return (the absolute value of) I<a> as a raw base-256 string in big- or
little-endian order.  The returned string has length I<nbytes>, which
defaults to the smallest nonzero length necessary to represent the
argument.  If I<nbytes> is too small, more significant bits are silently
discarded.  The sign is ignored.

=item I<a>B<-E<gt>storeb2c(>[I<nbytes>]B<)>

=item I<a>B<-E<gt>storel2c(>[I<nbytes>]B<)>

Return I<a> as a raw two's complement base-256 string in big- or
little-endian order.    The returned string has length I<nbytes>, which
defaults to the smallest nonzero length necessary to represent the
argument.  If I<nbytes> is too small, more significant bits are silently
discarded.  An alternative way to look at this operation is that the
returned value is congruent to I<a> modulo 2^(8 I<nbytes>).

=item I<a>B<-E<gt>tostring(>[I<radix>]B<)>

Return I<a> expressed as a base-I<radix> string.  The default I<radix>
is 10.  The result may be preceded by a minus sign if I<a> is negative.
No leading zeroes or other radix prefix are prepended.

=item I<a>B<-E<gt>toint()>

Returns I<a> as a Perl integer.  If I<a> cannot be represented as a Perl
integer then an undefined result is returned.

=back

=head2 Barrett reduction

Barrett reduction is an efficient way of doing modular reduction on
numbers which aren't too much bigger than the modulus.  Barrett
reduction isn't as efficient as Montgomery reduction (see below) but is
simpler and works on even moduli.

The precomputed information used to perform Barrett reduction are stored
in an object of type B<Catacomb::MP::Barrett>.

=over

=item B<Catacomb::MP::Barrett-E<gt>new(>I<m>B<)>

=item I<m>B<-E<gt>barrett()>

Create a new Barrett reduction context, set up for reducing modulo I<m>.

=item I<b>B<-E<gt>m()>

Returns the modulus stored in the reduction context I<b>.

=item I<b>B<-E<gt>reduce(>I<x>B<)>

Let I<m> be the modulus stored in the reduction context I<b>.  Then, if
I<x> is nonnegative and less then I<m>^2, returns I<x> mod I<m>.
Otherwise the return value is undefined.

=item I<b>B<-E<gt>exp(>I<x>B<,> I<n>B<)>

Let I<m> be the modulus stored in the reduction context I<b>.  Then
return I<x>^I<n> mod I<m>.  If I<n> is negative then I<x> must have an
inverse modulo I<m>; otherwise you'll get a big mess.

=back

=head2 Montgomery reduction

Montgomery reduction is a clever and efficient way of doing modular
reduction on numbers which aren't too much bigger than the modulus.  It
only works if the modulus is positive and odd, and it's a little
complicated.

Let I<m> be the modulus in question.  There is a value I<R> which is
computed in some way from I<m>.  The Montgomery reduction operation,
given a number I<x>, actually returns I<x> I<R>^(-1) mod I<m>, which
doesn't sound useful.  However, if you multiply all your numbers by a
factor of I<R> then the cleverness becomes clear: if you reduce I<x>
I<y> I<R>^2, you get I<x> I<y> I<R> mod I<m>.  Indeed, there's an
efficient multiply-and-reduce operation which takes two operands, each
with a factor of I<R> and returns a result with the factor of I<R>.  One
more reduction step drops off the final factor of I<R> to give you the
real answer.  Getting numbers into the required form involves
multiplying by I<R>^2 mod I<m> and reducing: conveniently, I<R>^2 mod
I<m> is precomputed.  Addition and subtraction leave the factor of I<R>
where it was, so you can just test-and-subtract to reduce.  One final
tip, if you're just multiplying two numbers, multiply them together
first, and then multiply the result by I<R>^2 to fix the answer.

If that all sounded too difficult, then you can use the B<in()> method to
convert to internal form, B<mul()> to multiply numbers in internal form,
B<expr()> and B<mexpr()> to do modular exponentiation, and B<out()> to
convert the result to external form.

The precomputed information used to perform Montgomery reduction are
stored in an object of type B<Catacomb::MP::Mont>.

=over

=item B<Catacomb::MP::Mont-E<gt>new(>I<m>B<)>

=item I<m>B<-E<gt>mont()>

Constructs a Montgomery reduction context for reduction modulo I<m>.
Returns B<undef> unless I<m> is an odd positive integer.

=item I<mm>B<-E<gt>r()>

Returns the Montgomery reduction factor I<R>.

=item I<mm>B<-E<gt>r2()>

Returns the Montgomerization factor I<R>^2.

=item I<mm>B<-E<gt>m()>

Returns the modulus I<m>.

=item I<mm>B<-E<gt>in(>I<x>B<)>

Returns the Montgomerized form of I<x>, i.e., I<x> I<R> mod I<m>.  I<x>
can be any integer.

=item I<mm>B<-E<gt>mul(>I<x>B<,> I<y>B<)>

Montgomery multiplication.  If I<x> and I<y> are Montgomerized, then
returns their Montgomerized product; i.e., I<x> I<y> I<R>^(-1) mod I<m>.

=item I<mm>B<-E<gt>expr(>I<x>B<,> I<n>B<)>

Montgomery exponentiation.  If I<x> is Montgomerized, then returns the
Montgomerized value of I<x>^I<n> mod I<m>; i.e., if I<x> = I<x'> I<R>,
it returns I<x'>^I<n> I<R> mod I<m>.  Here, if I<n> is negative, then
I<x> must have an inverse modulo I<m>; otherwise there'll be a big mess.

=item I<mm>B<-E<gt>mexpr(>I<x0>B<,> I<n0>B<,> I<x1>B<,> I<n1>B<,> ...B<)>

Simultaneous Montgomery exponentiation.  If the I<xi> are Montgomerized,
then returns the Montgomerized value of I<x0>^I<n0> I<x1>^I<n1> ... mod
I<m>; i.e., if I<xi> = I<xi'> I<R>, it returns I<x0'>^I<n0> I<x1'>^I<n1>
... I<R> mod I<m>.  Here, if some I<ni> is negative, then the
corresponding I<xi> must have an inverse modulo I<m>; otherwise there'll
be a big mess.

=item I<mm>B<-E<gt>reduce(>I<x>B<)>

=item I<mm>B<-E<gt>out(>I<x>B<)>

Returns the unMontgomerized form of I<x>; i.e., I<x> I<R>^(-1) mod I<m>.

=item I<mm>B<-E<gt>expr(>I<x>B<,> I<n>B<)>

Modular exponentiation.  Returns I<x>^I<n> mod I<m>.  Here, if I<n> is
negative, then I<x> must have an inverse modulo I<m>; otherwise there'll
be a big mess.

=item I<mm>B<-E<gt>mexpr(>I<x0>B<,> I<n0>B<,> I<x1>B<,> I<n1>B<,> ...B<)>

Simultaneous modular exponentiation.  Returns I<x0>^I<n0> I<x1>^I<n1>
... mod I<m>.  Here, if some I<ni> is negative, then the corresponding
I<xi> must have an inverse modulo I<m>; otherwise there'll be a big
mess.

=back

=head2 Nice prime reduction

For some numbers, we can do modular reduction quite rapidly.
Specifically, these are numbers of the form

=over

2^I<n0> - 2^I<n1> +/- 2^I<n2> +/- ... +/- 2^I<nk>

=back

where I<n0> > I<n1> > ... I<nk>, for small values of I<k>.  We call such
numbers `nice' (that's specific to Catacomb, not a term in general use).
Usually, numbers of interest to us will be prime; hence, we shall speak
of `nice primes'.

Information for efficient reduction modulo a nice number is stored in an
object of type B<Catacomb::MP::Reduce>.

=over

=item B<Catacomb::MP::Reduce-E<gt>new(>I<m>B<)>

=item I<m>B<-E<gt>mkreduce()>

Create a new nice-modulus reduction context, set up for reducing modulo
I<m>.

=item I<r>B<-E<gt>m()>

Returns the modulus stored in the reduction context I<r>.

=item I<r>B<-E<gt>reduce(>I<x>B<)>

Let I<m> be the modulus stored in the reduction context I<r>.  Then, if
I<x>, returns I<x> mod I<m>.  Otherwise the return value is undefined.

=item I<r>B<-E<gt>exp(>I<x>B<,> I<n>B<)>

Let I<m> be the modulus stored in the reduction context I<r>.  Then
return I<x>^I<n> mod I<m>.  If I<n> is negative then I<x> must have an
inverse modulo I<m>; otherwise you'll get a big mess.

=back

=head2 Chinese Remainder Theorem solution

Catacomb can solve simultaneous congruence problems, i.e., of the form
I<x> = I<ri> (mod I<ni>) using the Chinese Remainder Theorem (CRT).  It
is required that the I<ni> be positive and pairwise relatively prime;
otherwise you'll get an unpleasant mess.

=over

=item B<Catacomb::MP::CRT-E<gt>new(>I<n0>B<,> I<n1>B<,> ...B<)>

Construct a new CRT solving context, for solving congruences mod the
I<ni>, and return it.

=item I<crt>B<-E<gt>product()>

Returns the product of the moduli I<ni>.  In a scalar context, return
the number of moduli.

=item I<crt>B<-E<gt>moduli()>

Returns a list of the I<ni>, as passed to B<new()>.

=item I<crt>B<-E<gt>solve(>I<r0>B<,> I<r1>B<,>, ...B<)>

Returns the unique answer I<x> (modulo the product of the I<ni>) to the
congruences I<x> = I<ri> (mod I<ni>).  You will get an unhelpful answer
if any I<ri> >= I<ni>, and an error if you pass the wrong number of
I<ri> arguments.

=back

=head1 SEE ALSO

Catacomb(3).

=head1 AUTHOR

Mark Wooding, <mdw@nsict.org>