--- /dev/null
+=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>
+
~mdw / catacomb-perl / blobdiff