3 Catacomb::MP - Multiprecision integer arithmetic
7 use Catacomb qw(:const :mp :random :pgen);
9 $x = Catacomb::MP->new($str, [$radix]);
10 $x = Catacomb::MP->new($i);
11 $x = Catacomb::MP->loadb($bytes);
12 $x = Catacomb::MP->loadl($bytes);
13 $x = Catacomb::MP->loadb2c($bytes);
14 $x = Catacomb::MP->loadl2c($bytes);
15 ($x, $rest) = Catacomb::MP->fromstring($str, [$radix]);
16 $x = mp($str, [$radix]);
18 ($x, $rest) = mp_fromstring($str, [$radix]);
42 ($q, $r) = $a->div($b);
69 $a = $b->setbit2c($n);
70 $a = $b->clearbit($n);
71 $a = $b->clearbit2c($n);
73 $p = $b->testbit2c($n);
74 $x = $y->copy(); # largely useless
76 ($g, $u, $v) = $a->gcd($b);
77 ($s, $t) = $m->odd(); # m = 2^s t
78 $a = $p->modexp($x, $e);
82 $p = $x->primep([$rng]);
84 $nbytes = $x->octets();
85 $bytes = $x->storeb([$nbytes]);
86 $bytes = $x->storel([$nbytes]);
87 $bytes = $x->storeb2c([$nbytes]);
88 $bytes = $x->storel2c([$nbytes]);
89 $str = $x->tostring([$radix]);
92 $barrett = Catacomb::MP::Barrett->new($p);
93 $barrett = $p->barrett();
95 $x = $barrett->reduce($y);
96 $a = $barrett->exp($x, $y);
98 $mont = Catacomb::MP::Mont->new($p);
104 $a = $mont->mul($x, $y);
105 $a = $mont->expr($x, $y);
106 $a = $mont->mexpr($x0, $e0, $x1, $e1, ...);
107 $x = $mont->reduce($y);
109 $a = $mont->exp($x, $y);
110 $a = $mont->mexp($x0, $e0, $x1, $e1, ...);
112 $reduce = Catacomb::MP::Reduce->new($p);
113 $reduce = $p->mkreduce();
115 $x = $reduce->reduce($y);
116 $a = $barrett->exp($x, $y);
118 $crt = Catacomb::MP::CRT->new(@n);
119 $n = $crt->product();
121 $x = $crt->solve(@r);
125 This is a full-featured multiprecision maths library. Integer objects
126 belong to the class C<Catacomb::MP>. Most standard arithmetic operators
127 are overloaded and do the Right Thing. See below for more details.
129 Note that multiprecision integers are I<immutable>. No operation can
130 change the value of an integer.
136 =item B<Catacomb::MP-E<gt>new(>i<int>B<)>
138 =item B<mp(>I<int>B<)>
140 Return the multiprecision integer whose value is I<int>.
142 =item B<Catacomb::MP-E<gt>new(>I<string>, [I<radix>]B<)>
144 =item B<mp(>I<string>, [I<radix>]B<)>
146 Returns the multiprecision integer whose value is represented by
147 I<string>, in the given I<radix>. If I<radix> is 0, as it is by
148 default, then I<string> may have a prefix giving its radix: it may be
149 B<0> for octal, B<0x> for hex, or I<r>B<_> for radix I<r>. If the
150 argument is invalid, B<undef> is returned.
152 =item B<Catacomb::MP-E<gt>loadb(>I<bytes>B<)>
154 =item B<Catacomb::MP-E<gt>loadl(>I<bytes>B<)>
156 =item B<mp_loadb(>I<bytes>B<)>
158 =item B<mp_loadl(>I<bytes>B<)>
160 Returns the multiprecision integer whose value is represented in big- or
161 little-endian base-256 form by I<bytes>. This will always be
164 Returns the multiprecision integer whose value is represented in
165 little-endian base-256 form by I<bytes>. This will always be
168 =item B<Catacomb::MP-E<gt>loadb2c(>I<bytes>B<)>
170 =item B<Catacomb::MP-E<gt>loadl2c(>I<bytes>B<)>
172 =item B<mp_loadb2c(>I<bytes>B<)>
174 =item B<mp_loadl2c(>I<bytes>B<)>
176 Returns the multiprecision integer whose value is represented in
177 two's complement big- or little-endian base-256 form by I<bytes>.
179 =item B<Catacomb::MP-E<gt>fromstring(>I<string>, [I<radix>]B<)>
181 =item B<mp_fromstring(>I<string>, [I<radix>]B<)>
183 Returns the multiprecision integer whose value is represented by
184 I<string>, and (in a list context) the remainder of I<string>. If
185 I<radix> is 0 then I<string> may have a prefix, as for B<new()>. If
186 I<string> is invalid, an empty list is returned.
192 The B<Catacomb::MP> package overloads the standard arithmetic
193 operators. As long as one operand is a B<Catacomb::MP> object, the
194 other may be any argument acceptable to B<new()>. As a special
195 exception, either argument to B<*> may be a B<Catacomb::EC::Pt> object,
196 to perform elliptic curve point multiplication.
200 =item I<a>B<-E<gt>add(>I<b>B<)>
202 =item I<a>B<-E<gt>sub(>I<b>B<)>
204 =item I<a>B<-E<gt>mul(>I<b>B<)>
206 =item I<a>B<-E<gt>div(>I<b>B<)>
208 =item I<a>B<-E<gt>mod(>I<b>B<)>
210 The basic algorithms. The argument I<b> may be any value acceptable to
211 the B<new()> constructor.
213 In an array context, the B<div()> method returns a two values: the
214 quotient and remainder. The quotient returned is the I<floor> of the
215 true quotient: the remainder therefore has the same sign as the divisor.
216 This is usually the right result for number-theoretic purposes.
218 =item I<a>B<-E<gt>sqr()>
220 Returns the square of its argument. Usually faster than multiplying it
223 =item I<a>B<-E<gt>exp(>I<b>B<)>
225 Returns the result I<a>^I<b>. Theoretically, I<b> can be a
226 multiprecision integer, but the result will be impossibly huge if this
227 is the case. The exponent I<b> may not be negative.
229 =item I<a>B<-E<gt>sqrt()>
231 Returns the integer square-root of I<a> -- i.e., the largest integer
232 less than or equal to the true square root.
234 =item I<a>B<-E<gt>eq(>I<b>B<)>
236 Returns true if I<a> and I<b> are equal; or false otherwise.
238 =item I<a>B<-E<gt>cmp(>I<b>B<)>
240 Returns -1, 0 or +1 according to whether I<a> is less than, equal to or
243 =item I<a>B<-E<gt>and(>I<b>B<)>
245 =item I<a>B<-E<gt>and2c(>I<b>B<)>
247 =item I<a>B<-E<gt>or(>I<b>B<)>
249 =item I<a>B<-E<gt>or2c(>I<b>B<)>
251 =item I<a>B<-E<gt>xor(>I<b>B<)>
253 =item I<a>B<-E<gt>xor2c(>I<b>B<)>
255 =item I<a>B<-E<gt>nand(>I<b>B<)>
257 =item I<a>B<-E<gt>nand2c(>I<b>B<)>
259 =item I<a>B<-E<gt>nor(>I<b>B<)>
261 =item I<a>B<-E<gt>nor2c(>I<b>B<)>
263 =item I<a>B<-E<gt>not()>
265 =item I<a>B<-E<gt>not2c()>
267 Standard bitwise operators. The unadorned methods ignore the sign of
268 their operands and use only the absolute values. The B<2c> versions use
269 a two's complement representation, with the sign bit repeated
272 =item I<a>B<-E<gt>lsl(>I<n>B<)>
274 =item I<a>B<-E<gt>lsl2c(>I<n>B<)>
276 =item I<a>B<-E<gt>lsr(>I<n>B<)>
278 =item I<a>B<-E<gt>lsr2c(>I<n>B<)>
280 Standard shifting operators. The unadorned methods treat only absolute
281 values, while the B<2c> versions use a two's complement representation,
282 with the sign bit repeated infinitely.
284 =item I<a>B<-E<gt>setbit(>I<n>B<)>
286 =item I<a>B<-E<gt>setbit2c(>I<n>B<)>
288 =item I<a>B<-E<gt>clearbit(>I<n>B<)>
290 =item I<a>B<-E<gt>clearbit2c(>I<n>B<)>
292 Returns a copy of the argument I<a> with bit I<n> set or clear. The
293 unadorned methods treat only absolute values, while the B<2c> variants
294 use a two's complement representation with the sign bit repeated
297 =item I<a>B<-E<gt>testbit(>I<n>B<)>
299 =item I<a>B<-E<gt>testbit2c(>I<n>B<)>
301 Returns true or false according to whether I<a> has bit I<n> set or
302 clear. The unadorned methods treat only absolute values, while the
303 B<2c> variant uses a two's complement representation with the sign bit
306 =item I<a>B<-E<gt>copy()>
308 Returns a copy of I<a>. Not useful: use assignment instead.
310 =item I<a>B<-E<gt>gcd(>I<b>B<)>
312 Returns the greatest common divisor of I<a> and I<b>. This will never
313 be negative. In a list context, it also returns two further results,
314 I<u> and I<v>: if the GCD is I<g> then we have I<a> I<u> + I<b> I<v> =
315 I<g>. Furthermore, I<v> will have the same sign as I<b>.
317 =item I<a>B<-E<gt>odd()>
319 Returns a pair of numbers I<s> and I<t>, so that I<a> = I<t> << I<s>,
320 with I<t> odd and I<s> as large as possible. I<t> will be
321 multiprecision; I<s> will be a standard Perl integer.
323 =item I<a>B<-E<gt>modexp(>I<b>B<,> I<n>B<)>
325 Returns I<b>^I<n> mod I<a>.
327 =item I<a>B<-E<gt>modinv(>I<b>B<)>
329 Returns the inverse of I<b> modulo I<a>. If no inverse exists then
332 =item I<a>B<-E<gt>modsqrt(>I<b>B<)>
334 If I<a> is prime, and I<b> is a quadratic residue mod I<a>, then return
335 a square root of I<b> mod I<a>. If I<a> is not prime, does something
336 arbitrary; if I<b> is a nonresidue then report an error.
338 =item I<a>B<-E<gt>jac(>I<b>B<)>
340 Returns the Jacobi symbol (I<b>/I<a>). If I<a> is prime, then this is
341 I<b>^((I<a> - 1)/2), i.e., 0 if I<b> is a multiple of I<a>, +1 if I<b>
342 is a quadratic residue modulo I<a>, or -1 if I<b> is a quadratic
345 =item I<a>B<-E<gt>primep(>[I<rng>]I<)>
347 Returns true or false, depending on whether I<a> is (probably) prime.
348 If I<rng> is specified, it must be a random number generator -- see
351 =item I<a>B<-E<gt>bits()>
353 Returns the smallest number of bits required to represent the integer
356 =item I<a>B<-E<gt>octets()>
358 Returns the smallest number of bytes required to represent the integer
361 =item I<a>B<-E<gt>octets2c()>
363 Returns the smallest number of bytes required to represent the integer
364 I<a> as two's complement. This may be one greater than the result
365 returned by B<octets()> in the case where the top bit of I<a> is the top
368 =item I<a>B<-E<gt>storeb(>[I<nbytes>]B<)>
370 =item I<a>B<-E<gt>storel(>[I<nbytes>]B<)>
372 Return (the absolute value of) I<a> as a raw base-256 string in big- or
373 little-endian order. The returned string has length I<nbytes>, which
374 defaults to the smallest nonzero length necessary to represent the
375 argument. If I<nbytes> is too small, more significant bits are silently
376 discarded. The sign is ignored.
378 =item I<a>B<-E<gt>storeb2c(>[I<nbytes>]B<)>
380 =item I<a>B<-E<gt>storel2c(>[I<nbytes>]B<)>
382 Return I<a> as a raw two's complement base-256 string in big- or
383 little-endian order. The returned string has length I<nbytes>, which
384 defaults to the smallest nonzero length necessary to represent the
385 argument. If I<nbytes> is too small, more significant bits are silently
386 discarded. An alternative way to look at this operation is that the
387 returned value is congruent to I<a> modulo 2^(8 I<nbytes>).
389 =item I<a>B<-E<gt>tostring(>[I<radix>]B<)>
391 Return I<a> expressed as a base-I<radix> string. The default I<radix>
392 is 10. The result may be preceded by a minus sign if I<a> is negative.
393 No leading zeroes or other radix prefix are prepended.
395 =item I<a>B<-E<gt>toint()>
397 Returns I<a> as a Perl integer. If I<a> cannot be represented as a Perl
398 integer then an undefined result is returned.
402 =head2 Barrett reduction
404 Barrett reduction is an efficient way of doing modular reduction on
405 numbers which aren't too much bigger than the modulus. Barrett
406 reduction isn't as efficient as Montgomery reduction (see below) but is
407 simpler and works on even moduli.
409 The precomputed information used to perform Barrett reduction are stored
410 in an object of type B<Catacomb::MP::Barrett>.
414 =item B<Catacomb::MP::Barrett-E<gt>new(>I<m>B<)>
416 =item I<m>B<-E<gt>barrett()>
418 Create a new Barrett reduction context, set up for reducing modulo I<m>.
420 =item I<b>B<-E<gt>m()>
422 Returns the modulus stored in the reduction context I<b>.
424 =item I<b>B<-E<gt>reduce(>I<x>B<)>
426 Let I<m> be the modulus stored in the reduction context I<b>. Then, if
427 I<x> is nonnegative and less then I<m>^2, returns I<x> mod I<m>.
428 Otherwise the return value is undefined.
430 =item I<b>B<-E<gt>exp(>I<x>B<,> I<n>B<)>
432 Let I<m> be the modulus stored in the reduction context I<b>. Then
433 return I<x>^I<n> mod I<m>. If I<n> is negative then I<x> must have an
434 inverse modulo I<m>; otherwise you'll get a big mess.
438 =head2 Montgomery reduction
440 Montgomery reduction is a clever and efficient way of doing modular
441 reduction on numbers which aren't too much bigger than the modulus. It
442 only works if the modulus is positive and odd, and it's a little
445 Let I<m> be the modulus in question. There is a value I<R> which is
446 computed in some way from I<m>. The Montgomery reduction operation,
447 given a number I<x>, actually returns I<x> I<R>^(-1) mod I<m>, which
448 doesn't sound useful. However, if you multiply all your numbers by a
449 factor of I<R> then the cleverness becomes clear: if you reduce I<x>
450 I<y> I<R>^2, you get I<x> I<y> I<R> mod I<m>. Indeed, there's an
451 efficient multiply-and-reduce operation which takes two operands, each
452 with a factor of I<R> and returns a result with the factor of I<R>. One
453 more reduction step drops off the final factor of I<R> to give you the
454 real answer. Getting numbers into the required form involves
455 multiplying by I<R>^2 mod I<m> and reducing: conveniently, I<R>^2 mod
456 I<m> is precomputed. Addition and subtraction leave the factor of I<R>
457 where it was, so you can just test-and-subtract to reduce. One final
458 tip, if you're just multiplying two numbers, multiply them together
459 first, and then multiply the result by I<R>^2 to fix the answer.
461 If that all sounded too difficult, then you can use the B<in()> method to
462 convert to internal form, B<mul()> to multiply numbers in internal form,
463 B<expr()> and B<mexpr()> to do modular exponentiation, and B<out()> to
464 convert the result to external form.
466 The precomputed information used to perform Montgomery reduction are
467 stored in an object of type B<Catacomb::MP::Mont>.
471 =item B<Catacomb::MP::Mont-E<gt>new(>I<m>B<)>
473 =item I<m>B<-E<gt>mont()>
475 Constructs a Montgomery reduction context for reduction modulo I<m>.
476 Returns B<undef> unless I<m> is an odd positive integer.
478 =item I<mm>B<-E<gt>r()>
480 Returns the Montgomery reduction factor I<R>.
482 =item I<mm>B<-E<gt>r2()>
484 Returns the Montgomerization factor I<R>^2.
486 =item I<mm>B<-E<gt>m()>
488 Returns the modulus I<m>.
490 =item I<mm>B<-E<gt>in(>I<x>B<)>
492 Returns the Montgomerized form of I<x>, i.e., I<x> I<R> mod I<m>. I<x>
495 =item I<mm>B<-E<gt>mul(>I<x>B<,> I<y>B<)>
497 Montgomery multiplication. If I<x> and I<y> are Montgomerized, then
498 returns their Montgomerized product; i.e., I<x> I<y> I<R>^(-1) mod I<m>.
500 =item I<mm>B<-E<gt>expr(>I<x>B<,> I<n>B<)>
502 Montgomery exponentiation. If I<x> is Montgomerized, then returns the
503 Montgomerized value of I<x>^I<n> mod I<m>; i.e., if I<x> = I<x'> I<R>,
504 it returns I<x'>^I<n> I<R> mod I<m>. Here, if I<n> is negative, then
505 I<x> must have an inverse modulo I<m>; otherwise there'll be a big mess.
507 =item I<mm>B<-E<gt>mexpr(>I<x0>B<,> I<n0>B<,> I<x1>B<,> I<n1>B<,> ...B<)>
509 Simultaneous Montgomery exponentiation. If the I<xi> are Montgomerized,
510 then returns the Montgomerized value of I<x0>^I<n0> I<x1>^I<n1> ... mod
511 I<m>; i.e., if I<xi> = I<xi'> I<R>, it returns I<x0'>^I<n0> I<x1'>^I<n1>
512 ... I<R> mod I<m>. Here, if some I<ni> is negative, then the
513 corresponding I<xi> must have an inverse modulo I<m>; otherwise there'll
516 =item I<mm>B<-E<gt>reduce(>I<x>B<)>
518 =item I<mm>B<-E<gt>out(>I<x>B<)>
520 Returns the unMontgomerized form of I<x>; i.e., I<x> I<R>^(-1) mod I<m>.
522 =item I<mm>B<-E<gt>expr(>I<x>B<,> I<n>B<)>
524 Modular exponentiation. Returns I<x>^I<n> mod I<m>. Here, if I<n> is
525 negative, then I<x> must have an inverse modulo I<m>; otherwise there'll
528 =item I<mm>B<-E<gt>mexpr(>I<x0>B<,> I<n0>B<,> I<x1>B<,> I<n1>B<,> ...B<)>
530 Simultaneous modular exponentiation. Returns I<x0>^I<n0> I<x1>^I<n1>
531 ... mod I<m>. Here, if some I<ni> is negative, then the corresponding
532 I<xi> must have an inverse modulo I<m>; otherwise there'll be a big
537 =head2 Nice prime reduction
539 For some numbers, we can do modular reduction quite rapidly.
540 Specifically, these are numbers of the form
544 2^I<n0> - 2^I<n1> +/- 2^I<n2> +/- ... +/- 2^I<nk>
548 where I<n0> > I<n1> > ... I<nk>, for small values of I<k>. We call such
549 numbers `nice' (that's specific to Catacomb, not a term in general use).
550 Usually, numbers of interest to us will be prime; hence, we shall speak
553 Information for efficient reduction modulo a nice number is stored in an
554 object of type B<Catacomb::MP::Reduce>.
558 =item B<Catacomb::MP::Reduce-E<gt>new(>I<m>B<)>
560 =item I<m>B<-E<gt>mkreduce()>
562 Create a new nice-modulus reduction context, set up for reducing modulo
565 =item I<r>B<-E<gt>m()>
567 Returns the modulus stored in the reduction context I<r>.
569 =item I<r>B<-E<gt>reduce(>I<x>B<)>
571 Let I<m> be the modulus stored in the reduction context I<r>. Then, if
572 I<x>, returns I<x> mod I<m>. Otherwise the return value is undefined.
574 =item I<r>B<-E<gt>exp(>I<x>B<,> I<n>B<)>
576 Let I<m> be the modulus stored in the reduction context I<r>. Then
577 return I<x>^I<n> mod I<m>. If I<n> is negative then I<x> must have an
578 inverse modulo I<m>; otherwise you'll get a big mess.
582 =head2 Chinese Remainder Theorem solution
584 Catacomb can solve simultaneous congruence problems, i.e., of the form
585 I<x> = I<ri> (mod I<ni>) using the Chinese Remainder Theorem (CRT). It
586 is required that the I<ni> be positive and pairwise relatively prime;
587 otherwise you'll get an unpleasant mess.
591 =item B<Catacomb::MP::CRT-E<gt>new(>I<n0>B<,> I<n1>B<,> ...B<)>
593 Construct a new CRT solving context, for solving congruences mod the
594 I<ni>, and return it.
596 =item I<crt>B<-E<gt>product()>
598 Returns the product of the moduli I<ni>. In a scalar context, return
599 the number of moduli.
601 =item I<crt>B<-E<gt>moduli()>
603 Returns a list of the I<ni>, as passed to B<new()>.
605 =item I<crt>B<-E<gt>solve(>I<r0>B<,> I<r1>B<,>, ...B<)>
607 Returns the unique answer I<x> (modulo the product of the I<ni>) to the
608 congruences I<x> = I<ri> (mod I<ni>). You will get an unhelpful answer
609 if any I<ri> >= I<ni>, and an error if you pass the wrong number of
620 Mark Wooding, <mdw@nsict.org>