2 /* vim: set expandtab tabstop=4 shiftwidth=4 softtabstop=4: */
5 * Pure-PHP arbitrary precision integer arithmetic library.
7 * Supports base-2, base-10, base-16, and base-256 numbers. Uses the GMP or BCMath extensions, if available,
8 * and an internal implementation, otherwise.
10 * PHP versions 4 and 5
12 * {@internal (all DocBlock comments regarding implementation - such as the one that follows - refer to the
13 * {@link MATH_BIGINTEGER_MODE_INTERNAL MATH_BIGINTEGER_MODE_INTERNAL} mode)
15 * Math_BigInteger uses base-2**26 to perform operations such as multiplication and division and
16 * base-2**52 (ie. two base 2**26 digits) to perform addition and subtraction. Because the largest possible
17 * value when multiplying two base-2**26 numbers together is a base-2**52 number, double precision floating
18 * point numbers - numbers that should be supported on most hardware and whose significand is 53 bits - are
19 * used. As a consequence, bitwise operators such as >> and << cannot be used, nor can the modulo operator %,
20 * which only supports integers. Although this fact will slow this library down, the fact that such a high
21 * base is being used should more than compensate.
23 * When PHP version 6 is officially released, we'll be able to use 64-bit integers. This should, once again,
24 * allow bitwise operators, and will increase the maximum possible base to 2**31 (or 2**62 for addition /
27 * Numbers are stored in {@link http://en.wikipedia.org/wiki/Endianness little endian} format. ie.
28 * (new Math_BigInteger(pow(2, 26)))->value = array(0, 1)
30 * Useful resources are as follows:
32 * - {@link http://www.cacr.math.uwaterloo.ca/hac/about/chap14.pdf Handbook of Applied Cryptography (HAC)}
33 * - {@link http://math.libtomcrypt.com/files/tommath.pdf Multi-Precision Math (MPM)}
34 * - Java's BigInteger classes. See /j2se/src/share/classes/java/math in jdk-1_5_0-src-jrl.zip
36 * Here's an example of how to use this library:
39 * include('Math/BigInteger.php');
41 * $a = new Math_BigInteger(2);
42 * $b = new Math_BigInteger(3);
46 * echo $c->toString(); // outputs 5
50 * LICENSE: This library is free software; you can redistribute it and/or
51 * modify it under the terms of the GNU Lesser General Public
52 * License as published by the Free Software Foundation; either
53 * version 2.1 of the License, or (at your option) any later version.
55 * This library is distributed in the hope that it will be useful,
56 * but WITHOUT ANY WARRANTY; without even the implied warranty of
57 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
58 * Lesser General Public License for more details.
60 * You should have received a copy of the GNU Lesser General Public
61 * License along with this library; if not, write to the Free Software
62 * Foundation, Inc., 59 Temple Place, Suite 330, Boston,
66 * @package Math_BigInteger
67 * @author Jim Wigginton <terrafrost@php.net>
68 * @copyright MMVI Jim Wigginton
69 * @license http://www.gnu.org/licenses/lgpl.txt
70 * @version $Id: BigInteger.php,v 1.33 2010/03/22 22:32:03 terrafrost Exp $
71 * @link http://pear.php.net/package/Math_BigInteger
78 * @see Math_BigInteger::_reduce()
81 * @see Math_BigInteger::_montgomery()
82 * @see Math_BigInteger::_prepMontgomery()
84 define('MATH_BIGINTEGER_MONTGOMERY', 0);
86 * @see Math_BigInteger::_barrett()
88 define('MATH_BIGINTEGER_BARRETT', 1);
90 * @see Math_BigInteger::_mod2()
92 define('MATH_BIGINTEGER_POWEROF2', 2);
94 * @see Math_BigInteger::_remainder()
96 define('MATH_BIGINTEGER_CLASSIC', 3);
98 * @see Math_BigInteger::__clone()
100 define('MATH_BIGINTEGER_NONE', 4);
106 * Rather than create a thousands and thousands of new Math_BigInteger objects in repeated function calls to add() and
107 * multiply() or whatever, we'll just work directly on arrays, taking them in as parameters and returning them.
112 * $result[MATH_BIGINTEGER_VALUE] contains the value.
114 define('MATH_BIGINTEGER_VALUE', 0);
116 * $result[MATH_BIGINTEGER_SIGN] contains the sign.
118 define('MATH_BIGINTEGER_SIGN', 1);
123 * @see Math_BigInteger::_montgomery()
124 * @see Math_BigInteger::_barrett()
129 * $cache[MATH_BIGINTEGER_VARIABLE] tells us whether or not the cached data is still valid.
131 define('MATH_BIGINTEGER_VARIABLE', 0);
133 * $cache[MATH_BIGINTEGER_DATA] contains the cached data.
135 define('MATH_BIGINTEGER_DATA', 1);
142 * @see Math_BigInteger::Math_BigInteger()
145 * To use the pure-PHP implementation
147 define('MATH_BIGINTEGER_MODE_INTERNAL', 1);
149 * To use the BCMath library
151 * (if enabled; otherwise, the internal implementation will be used)
153 define('MATH_BIGINTEGER_MODE_BCMATH', 2);
155 * To use the GMP library
157 * (if present; otherwise, either the BCMath or the internal implementation will be used)
159 define('MATH_BIGINTEGER_MODE_GMP', 3);
163 * The largest digit that may be used in addition / subtraction
165 * (we do pow(2, 52) instead of using 4503599627370496, directly, because some PHP installations
166 * will truncate 4503599627370496)
170 define('MATH_BIGINTEGER_MAX_DIGIT52', pow(2, 52));
175 * At what point do we switch between Karatsuba multiplication and schoolbook long multiplication?
179 define('MATH_BIGINTEGER_KARATSUBA_CUTOFF', 25);
182 * Pure-PHP arbitrary precision integer arithmetic library. Supports base-2, base-10, base-16, and base-256
185 * @author Jim Wigginton <terrafrost@php.net>
188 * @package Math_BigInteger
190 class Math_BigInteger {
192 * Holds the BigInteger's value.
200 * Holds the BigInteger's magnitude.
205 var $is_negative = false;
208 * Random number generator function
210 * @see setRandomGenerator()
213 var $generator = 'mt_rand';
218 * @see setPrecision()
226 * @see setPrecision()
229 var $bitmask = false;
232 * Mode independant value used for serialization.
234 * If the bcmath or gmp extensions are installed $this->value will be a non-serializable resource, hence the need for
235 * a variable that'll be serializable regardless of whether or not extensions are being used. Unlike $this->value,
236 * however, $this->hex is only calculated when $this->__sleep() is called.
246 * Converts base-2, base-10, base-16, and binary strings (eg. base-256) to BigIntegers.
248 * If the second parameter - $base - is negative, then it will be assumed that the number's are encoded using
249 * two's compliment. The sole exception to this is -10, which is treated the same as 10 is.
254 * include('Math/BigInteger.php');
256 * $a = new Math_BigInteger('0x32', 16); // 50 in base-16
258 * echo $a->toString(); // outputs 50
262 * @param optional $x base-10 number or base-$base number if $base set.
263 * @param optional integer $base
264 * @return Math_BigInteger
267 function Math_BigInteger($x = 0, $base = 10)
269 if ( !defined('MATH_BIGINTEGER_MODE') ) {
271 case extension_loaded('gmp'):
272 define('MATH_BIGINTEGER_MODE', MATH_BIGINTEGER_MODE_GMP);
274 case extension_loaded('bcmath'):
275 define('MATH_BIGINTEGER_MODE', MATH_BIGINTEGER_MODE_BCMATH);
278 define('MATH_BIGINTEGER_MODE', MATH_BIGINTEGER_MODE_INTERNAL);
282 switch ( MATH_BIGINTEGER_MODE ) {
283 case MATH_BIGINTEGER_MODE_GMP:
284 if (is_resource($x) && get_resource_type($x) == 'GMP integer') {
288 $this->value = gmp_init(0);
290 case MATH_BIGINTEGER_MODE_BCMATH:
294 $this->value = array();
303 if (ord($x[0]) & 0x80) {
305 $this->is_negative = true;
308 switch ( MATH_BIGINTEGER_MODE ) {
309 case MATH_BIGINTEGER_MODE_GMP:
310 $sign = $this->is_negative ? '-' : '';
311 $this->value = gmp_init($sign . '0x' . bin2hex($x));
313 case MATH_BIGINTEGER_MODE_BCMATH:
314 // round $len to the nearest 4 (thanks, DavidMJ!)
315 $len = (strlen($x) + 3) & 0xFFFFFFFC;
317 $x = str_pad($x, $len, chr(0), STR_PAD_LEFT);
319 for ($i = 0; $i < $len; $i+= 4) {
320 $this->value = bcmul($this->value, '4294967296', 0); // 4294967296 == 2**32
321 $this->value = bcadd($this->value, 0x1000000 * ord($x[$i]) + ((ord($x[$i + 1]) << 16) | (ord($x[$i + 2]) << 8) | ord($x[$i + 3])), 0);
324 if ($this->is_negative) {
325 $this->value = '-' . $this->value;
329 // converts a base-2**8 (big endian / msb) number to base-2**26 (little endian / lsb)
332 $this->value[] = $this->_bytes2int($this->_base256_rshift($x, 26));
336 if ($this->is_negative) {
337 if (MATH_BIGINTEGER_MODE != MATH_BIGINTEGER_MODE_INTERNAL) {
338 $this->is_negative = false;
340 $temp = $this->add(new Math_BigInteger('-1'));
341 $this->value = $temp->value;
346 if ($base > 0 && $x[0] == '-') {
347 $this->is_negative = true;
351 $x = preg_replace('#^(?:0x)?([A-Fa-f0-9]*).*#', '$1', $x);
353 $is_negative = false;
354 if ($base < 0 && hexdec($x[0]) >= 8) {
355 $this->is_negative = $is_negative = true;
356 $x = bin2hex(~pack('H*', $x));
359 switch ( MATH_BIGINTEGER_MODE ) {
360 case MATH_BIGINTEGER_MODE_GMP:
361 $temp = $this->is_negative ? '-0x' . $x : '0x' . $x;
362 $this->value = gmp_init($temp);
363 $this->is_negative = false;
365 case MATH_BIGINTEGER_MODE_BCMATH:
366 $x = ( strlen($x) & 1 ) ? '0' . $x : $x;
367 $temp = new Math_BigInteger(pack('H*', $x), 256);
368 $this->value = $this->is_negative ? '-' . $temp->value : $temp->value;
369 $this->is_negative = false;
372 $x = ( strlen($x) & 1 ) ? '0' . $x : $x;
373 $temp = new Math_BigInteger(pack('H*', $x), 256);
374 $this->value = $temp->value;
378 $temp = $this->add(new Math_BigInteger('-1'));
379 $this->value = $temp->value;
384 $x = preg_replace('#^(-?[0-9]*).*#', '$1', $x);
386 switch ( MATH_BIGINTEGER_MODE ) {
387 case MATH_BIGINTEGER_MODE_GMP:
388 $this->value = gmp_init($x);
390 case MATH_BIGINTEGER_MODE_BCMATH:
391 // explicitly casting $x to a string is necessary, here, since doing $x[0] on -1 yields different
392 // results then doing it on '-1' does (modInverse does $x[0])
393 $this->value = (string) $x;
396 $temp = new Math_BigInteger();
398 // array(10000000) is 10**7 in base-2**26. 10**7 is the closest to 2**26 we can get without passing it.
399 $multiplier = new Math_BigInteger();
400 $multiplier->value = array(10000000);
403 $this->is_negative = true;
407 $x = str_pad($x, strlen($x) + (6 * strlen($x)) % 7, 0, STR_PAD_LEFT);
410 $temp = $temp->multiply($multiplier);
411 $temp = $temp->add(new Math_BigInteger($this->_int2bytes(substr($x, 0, 7)), 256));
415 $this->value = $temp->value;
418 case 2: // base-2 support originally implemented by Lluis Pamies - thanks!
420 if ($base > 0 && $x[0] == '-') {
421 $this->is_negative = true;
425 $x = preg_replace('#^([01]*).*#', '$1', $x);
426 $x = str_pad($x, strlen($x) + (3 * strlen($x)) % 4, 0, STR_PAD_LEFT);
430 $part = substr($x, 0, 4);
431 $str.= dechex(bindec($part));
435 if ($this->is_negative) {
439 $temp = new Math_BigInteger($str, 8 * $base); // ie. either -16 or +16
440 $this->value = $temp->value;
441 $this->is_negative = $temp->is_negative;
445 // base not supported, so we'll let $this == 0
450 * Converts a BigInteger to a byte string (eg. base-256).
452 * Negative numbers are saved as positive numbers, unless $twos_compliment is set to true, at which point, they're
453 * saved as two's compliment.
458 * include('Math/BigInteger.php');
460 * $a = new Math_BigInteger('65');
462 * echo $a->toBytes(); // outputs chr(65)
466 * @param Boolean $twos_compliment
469 * @internal Converts a base-2**26 number to base-2**8
471 function toBytes($twos_compliment = false)
473 if ($twos_compliment) {
474 $comparison = $this->compare(new Math_BigInteger());
475 if ($comparison == 0) {
476 return $this->precision > 0 ? str_repeat(chr(0), ($this->precision + 1) >> 3) : '';
479 $temp = $comparison < 0 ? $this->add(new Math_BigInteger(1)) : $this->copy();
480 $bytes = $temp->toBytes();
482 if (empty($bytes)) { // eg. if the number we're trying to convert is -1
486 if (ord($bytes[0]) & 0x80) {
487 $bytes = chr(0) . $bytes;
490 return $comparison < 0 ? ~$bytes : $bytes;
493 switch ( MATH_BIGINTEGER_MODE ) {
494 case MATH_BIGINTEGER_MODE_GMP:
495 if (gmp_cmp($this->value, gmp_init(0)) == 0) {
496 return $this->precision > 0 ? str_repeat(chr(0), ($this->precision + 1) >> 3) : '';
499 $temp = gmp_strval(gmp_abs($this->value), 16);
500 $temp = ( strlen($temp) & 1 ) ? '0' . $temp : $temp;
501 $temp = pack('H*', $temp);
503 return $this->precision > 0 ?
504 substr(str_pad($temp, $this->precision >> 3, chr(0), STR_PAD_LEFT), -($this->precision >> 3)) :
505 ltrim($temp, chr(0));
506 case MATH_BIGINTEGER_MODE_BCMATH:
507 if ($this->value === '0') {
508 return $this->precision > 0 ? str_repeat(chr(0), ($this->precision + 1) >> 3) : '';
512 $current = $this->value;
514 if ($current[0] == '-') {
515 $current = substr($current, 1);
518 while (bccomp($current, '0', 0) > 0) {
519 $temp = bcmod($current, '16777216');
520 $value = chr($temp >> 16) . chr($temp >> 8) . chr($temp) . $value;
521 $current = bcdiv($current, '16777216', 0);
524 return $this->precision > 0 ?
525 substr(str_pad($value, $this->precision >> 3, chr(0), STR_PAD_LEFT), -($this->precision >> 3)) :
526 ltrim($value, chr(0));
529 if (!count($this->value)) {
530 return $this->precision > 0 ? str_repeat(chr(0), ($this->precision + 1) >> 3) : '';
532 $result = $this->_int2bytes($this->value[count($this->value) - 1]);
534 $temp = $this->copy();
536 for ($i = count($temp->value) - 2; $i >= 0; --$i) {
537 $temp->_base256_lshift($result, 26);
538 $result = $result | str_pad($temp->_int2bytes($temp->value[$i]), strlen($result), chr(0), STR_PAD_LEFT);
541 return $this->precision > 0 ?
542 str_pad(substr($result, -(($this->precision + 7) >> 3)), ($this->precision + 7) >> 3, chr(0), STR_PAD_LEFT) :
547 * Converts a BigInteger to a hex string (eg. base-16)).
549 * Negative numbers are saved as positive numbers, unless $twos_compliment is set to true, at which point, they're
550 * saved as two's compliment.
555 * include('Math/BigInteger.php');
557 * $a = new Math_BigInteger('65');
559 * echo $a->toHex(); // outputs '41'
563 * @param Boolean $twos_compliment
566 * @internal Converts a base-2**26 number to base-2**8
568 function toHex($twos_compliment = false)
570 return bin2hex($this->toBytes($twos_compliment));
574 * Converts a BigInteger to a bit string (eg. base-2).
576 * Negative numbers are saved as positive numbers, unless $twos_compliment is set to true, at which point, they're
577 * saved as two's compliment.
582 * include('Math/BigInteger.php');
584 * $a = new Math_BigInteger('65');
586 * echo $a->toBits(); // outputs '1000001'
590 * @param Boolean $twos_compliment
593 * @internal Converts a base-2**26 number to base-2**2
595 function toBits($twos_compliment = false)
597 $hex = $this->toHex($twos_compliment);
599 for ($i = 0; $i < strlen($hex); $i+=8) {
600 $bits.= str_pad(decbin(hexdec(substr($hex, $i, 8))), 32, '0', STR_PAD_LEFT);
602 return $this->precision > 0 ? substr($bits, -$this->precision) : ltrim($bits, '0');
606 * Converts a BigInteger to a base-10 number.
611 * include('Math/BigInteger.php');
613 * $a = new Math_BigInteger('50');
615 * echo $a->toString(); // outputs 50
621 * @internal Converts a base-2**26 number to base-10**7 (which is pretty much base-10)
625 switch ( MATH_BIGINTEGER_MODE ) {
626 case MATH_BIGINTEGER_MODE_GMP:
627 return gmp_strval($this->value);
628 case MATH_BIGINTEGER_MODE_BCMATH:
629 if ($this->value === '0') {
633 return ltrim($this->value, '0');
636 if (!count($this->value)) {
640 $temp = $this->copy();
641 $temp->is_negative = false;
643 $divisor = new Math_BigInteger();
644 $divisor->value = array(10000000); // eg. 10**7
646 while (count($temp->value)) {
647 list($temp, $mod) = $temp->divide($divisor);
648 $result = str_pad(isset($mod->value[0]) ? $mod->value[0] : '', 7, '0', STR_PAD_LEFT) . $result;
650 $result = ltrim($result, '0');
651 if (empty($result)) {
655 if ($this->is_negative) {
656 $result = '-' . $result;
665 * PHP5 passes objects by reference while PHP4 passes by value. As such, we need a function to guarantee
666 * that all objects are passed by value, when appropriate. More information can be found here:
668 * {@link http://php.net/language.oop5.basic#51624}
672 * @return Math_BigInteger
676 $temp = new Math_BigInteger();
677 $temp->value = $this->value;
678 $temp->is_negative = $this->is_negative;
679 $temp->generator = $this->generator;
680 $temp->precision = $this->precision;
681 $temp->bitmask = $this->bitmask;
686 * __toString() magic method
688 * Will be called, automatically, if you're supporting just PHP5. If you're supporting PHP4, you'll need to call
692 * @internal Implemented per a suggestion by Techie-Michael - thanks!
694 function __toString()
696 return $this->toString();
700 * __clone() magic method
702 * Although you can call Math_BigInteger::__toString() directly in PHP5, you cannot call Math_BigInteger::__clone()
703 * directly in PHP5. You can in PHP4 since it's not a magic method, but in PHP5, you have to call it by using the PHP5
704 * only syntax of $y = clone $x. As such, if you're trying to write an application that works on both PHP4 and PHP5,
705 * call Math_BigInteger::copy(), instead.
709 * @return Math_BigInteger
713 return $this->copy();
717 * __sleep() magic method
719 * Will be called, automatically, when serialize() is called on a Math_BigInteger object.
726 $this->hex = $this->toHex(true);
727 $vars = array('hex');
728 if ($this->generator != 'mt_rand') {
729 $vars[] = 'generator';
731 if ($this->precision > 0) {
732 $vars[] = 'precision';
739 * __wakeup() magic method
741 * Will be called, automatically, when unserialize() is called on a Math_BigInteger object.
748 $temp = new Math_BigInteger($this->hex, -16);
749 $this->value = $temp->value;
750 $this->is_negative = $temp->is_negative;
751 $this->setRandomGenerator($this->generator);
752 if ($this->precision > 0) {
753 // recalculate $this->bitmask
754 $this->setPrecision($this->precision);
759 * Adds two BigIntegers.
764 * include('Math/BigInteger.php');
766 * $a = new Math_BigInteger('10');
767 * $b = new Math_BigInteger('20');
771 * echo $c->toString(); // outputs 30
775 * @param Math_BigInteger $y
776 * @return Math_BigInteger
778 * @internal Performs base-2**52 addition
782 switch ( MATH_BIGINTEGER_MODE ) {
783 case MATH_BIGINTEGER_MODE_GMP:
784 $temp = new Math_BigInteger();
785 $temp->value = gmp_add($this->value, $y->value);
787 return $this->_normalize($temp);
788 case MATH_BIGINTEGER_MODE_BCMATH:
789 $temp = new Math_BigInteger();
790 $temp->value = bcadd($this->value, $y->value, 0);
792 return $this->_normalize($temp);
795 $temp = $this->_add($this->value, $this->is_negative, $y->value, $y->is_negative);
797 $result = new Math_BigInteger();
798 $result->value = $temp[MATH_BIGINTEGER_VALUE];
799 $result->is_negative = $temp[MATH_BIGINTEGER_SIGN];
801 return $this->_normalize($result);
807 * @param Array $x_value
808 * @param Boolean $x_negative
809 * @param Array $y_value
810 * @param Boolean $y_negative
814 function _add($x_value, $x_negative, $y_value, $y_negative)
816 $x_size = count($x_value);
817 $y_size = count($y_value);
821 MATH_BIGINTEGER_VALUE => $y_value,
822 MATH_BIGINTEGER_SIGN => $y_negative
824 } else if ($y_size == 0) {
826 MATH_BIGINTEGER_VALUE => $x_value,
827 MATH_BIGINTEGER_SIGN => $x_negative
831 // subtract, if appropriate
832 if ( $x_negative != $y_negative ) {
833 if ( $x_value == $y_value ) {
835 MATH_BIGINTEGER_VALUE => array(),
836 MATH_BIGINTEGER_SIGN => false
840 $temp = $this->_subtract($x_value, false, $y_value, false);
841 $temp[MATH_BIGINTEGER_SIGN] = $this->_compare($x_value, false, $y_value, false) > 0 ?
842 $x_negative : $y_negative;
847 if ($x_size < $y_size) {
855 $value[] = 0; // just in case the carry adds an extra digit
858 for ($i = 0, $j = 1; $j < $size; $i+=2, $j+=2) {
859 $sum = $x_value[$j] * 0x4000000 + $x_value[$i] + $y_value[$j] * 0x4000000 + $y_value[$i] + $carry;
860 $carry = $sum >= MATH_BIGINTEGER_MAX_DIGIT52; // eg. floor($sum / 2**52); only possible values (in any base) are 0 and 1
861 $sum = $carry ? $sum - MATH_BIGINTEGER_MAX_DIGIT52 : $sum;
863 $temp = (int) ($sum / 0x4000000);
865 $value[$i] = (int) ($sum - 0x4000000 * $temp); // eg. a faster alternative to fmod($sum, 0x4000000)
869 if ($j == $size) { // ie. if $y_size is odd
870 $sum = $x_value[$i] + $y_value[$i] + $carry;
871 $carry = $sum >= 0x4000000;
872 $value[$i] = $carry ? $sum - 0x4000000 : $sum;
873 ++$i; // ie. let $i = $j since we've just done $value[$i]
877 for (; $value[$i] == 0x3FFFFFF; ++$i) {
884 MATH_BIGINTEGER_VALUE => $this->_trim($value),
885 MATH_BIGINTEGER_SIGN => $x_negative
890 * Subtracts two BigIntegers.
895 * include('Math/BigInteger.php');
897 * $a = new Math_BigInteger('10');
898 * $b = new Math_BigInteger('20');
900 * $c = $a->subtract($b);
902 * echo $c->toString(); // outputs -10
906 * @param Math_BigInteger $y
907 * @return Math_BigInteger
909 * @internal Performs base-2**52 subtraction
911 function subtract($y)
913 switch ( MATH_BIGINTEGER_MODE ) {
914 case MATH_BIGINTEGER_MODE_GMP:
915 $temp = new Math_BigInteger();
916 $temp->value = gmp_sub($this->value, $y->value);
918 return $this->_normalize($temp);
919 case MATH_BIGINTEGER_MODE_BCMATH:
920 $temp = new Math_BigInteger();
921 $temp->value = bcsub($this->value, $y->value, 0);
923 return $this->_normalize($temp);
926 $temp = $this->_subtract($this->value, $this->is_negative, $y->value, $y->is_negative);
928 $result = new Math_BigInteger();
929 $result->value = $temp[MATH_BIGINTEGER_VALUE];
930 $result->is_negative = $temp[MATH_BIGINTEGER_SIGN];
932 return $this->_normalize($result);
936 * Performs subtraction.
938 * @param Array $x_value
939 * @param Boolean $x_negative
940 * @param Array $y_value
941 * @param Boolean $y_negative
945 function _subtract($x_value, $x_negative, $y_value, $y_negative)
947 $x_size = count($x_value);
948 $y_size = count($y_value);
952 MATH_BIGINTEGER_VALUE => $y_value,
953 MATH_BIGINTEGER_SIGN => !$y_negative
955 } else if ($y_size == 0) {
957 MATH_BIGINTEGER_VALUE => $x_value,
958 MATH_BIGINTEGER_SIGN => $x_negative
962 // add, if appropriate (ie. -$x - +$y or +$x - -$y)
963 if ( $x_negative != $y_negative ) {
964 $temp = $this->_add($x_value, false, $y_value, false);
965 $temp[MATH_BIGINTEGER_SIGN] = $x_negative;
970 $diff = $this->_compare($x_value, $x_negative, $y_value, $y_negative);
974 MATH_BIGINTEGER_VALUE => array(),
975 MATH_BIGINTEGER_SIGN => false
979 // switch $x and $y around, if appropriate.
980 if ( (!$x_negative && $diff < 0) || ($x_negative && $diff > 0) ) {
985 $x_negative = !$x_negative;
987 $x_size = count($x_value);
988 $y_size = count($y_value);
991 // at this point, $x_value should be at least as big as - if not bigger than - $y_value
994 for ($i = 0, $j = 1; $j < $y_size; $i+=2, $j+=2) {
995 $sum = $x_value[$j] * 0x4000000 + $x_value[$i] - $y_value[$j] * 0x4000000 - $y_value[$i] - $carry;
996 $carry = $sum < 0; // eg. floor($sum / 2**52); only possible values (in any base) are 0 and 1
997 $sum = $carry ? $sum + MATH_BIGINTEGER_MAX_DIGIT52 : $sum;
999 $temp = (int) ($sum / 0x4000000);
1001 $x_value[$i] = (int) ($sum - 0x4000000 * $temp);
1002 $x_value[$j] = $temp;
1005 if ($j == $y_size) { // ie. if $y_size is odd
1006 $sum = $x_value[$i] - $y_value[$i] - $carry;
1008 $x_value[$i] = $carry ? $sum + 0x4000000 : $sum;
1013 for (; !$x_value[$i]; ++$i) {
1014 $x_value[$i] = 0x3FFFFFF;
1020 MATH_BIGINTEGER_VALUE => $this->_trim($x_value),
1021 MATH_BIGINTEGER_SIGN => $x_negative
1026 * Multiplies two BigIntegers
1028 * Here's an example:
1031 * include('Math/BigInteger.php');
1033 * $a = new Math_BigInteger('10');
1034 * $b = new Math_BigInteger('20');
1036 * $c = $a->multiply($b);
1038 * echo $c->toString(); // outputs 200
1042 * @param Math_BigInteger $x
1043 * @return Math_BigInteger
1046 function multiply($x)
1048 switch ( MATH_BIGINTEGER_MODE ) {
1049 case MATH_BIGINTEGER_MODE_GMP:
1050 $temp = new Math_BigInteger();
1051 $temp->value = gmp_mul($this->value, $x->value);
1053 return $this->_normalize($temp);
1054 case MATH_BIGINTEGER_MODE_BCMATH:
1055 $temp = new Math_BigInteger();
1056 $temp->value = bcmul($this->value, $x->value, 0);
1058 return $this->_normalize($temp);
1061 $temp = $this->_multiply($this->value, $this->is_negative, $x->value, $x->is_negative);
1063 $product = new Math_BigInteger();
1064 $product->value = $temp[MATH_BIGINTEGER_VALUE];
1065 $product->is_negative = $temp[MATH_BIGINTEGER_SIGN];
1067 return $this->_normalize($product);
1071 * Performs multiplication.
1073 * @param Array $x_value
1074 * @param Boolean $x_negative
1075 * @param Array $y_value
1076 * @param Boolean $y_negative
1080 function _multiply($x_value, $x_negative, $y_value, $y_negative)
1082 //if ( $x_value == $y_value ) {
1084 // MATH_BIGINTEGER_VALUE => $this->_square($x_value),
1085 // MATH_BIGINTEGER_SIGN => $x_sign != $y_value
1089 $x_length = count($x_value);
1090 $y_length = count($y_value);
1092 if ( !$x_length || !$y_length ) { // a 0 is being multiplied
1094 MATH_BIGINTEGER_VALUE => array(),
1095 MATH_BIGINTEGER_SIGN => false
1100 MATH_BIGINTEGER_VALUE => min($x_length, $y_length) < 2 * MATH_BIGINTEGER_KARATSUBA_CUTOFF ?
1101 $this->_trim($this->_regularMultiply($x_value, $y_value)) :
1102 $this->_trim($this->_karatsuba($x_value, $y_value)),
1103 MATH_BIGINTEGER_SIGN => $x_negative != $y_negative
1108 * Performs long multiplication on two BigIntegers
1110 * Modeled after 'multiply' in MutableBigInteger.java.
1112 * @param Array $x_value
1113 * @param Array $y_value
1117 function _regularMultiply($x_value, $y_value)
1119 $x_length = count($x_value);
1120 $y_length = count($y_value);
1122 if ( !$x_length || !$y_length ) { // a 0 is being multiplied
1126 if ( $x_length < $y_length ) {
1128 $x_value = $y_value;
1131 $x_length = count($x_value);
1132 $y_length = count($y_value);
1135 $product_value = $this->_array_repeat(0, $x_length + $y_length);
1137 // the following for loop could be removed if the for loop following it
1138 // (the one with nested for loops) initially set $i to 0, but
1139 // doing so would also make the result in one set of unnecessary adds,
1140 // since on the outermost loops first pass, $product->value[$k] is going
1145 for ($j = 0; $j < $x_length; ++$j) { // ie. $i = 0
1146 $temp = $x_value[$j] * $y_value[0] + $carry; // $product_value[$k] == 0
1147 $carry = (int) ($temp / 0x4000000);
1148 $product_value[$j] = (int) ($temp - 0x4000000 * $carry);
1151 $product_value[$j] = $carry;
1153 // the above for loop is what the previous comment was talking about. the
1154 // following for loop is the "one with nested for loops"
1155 for ($i = 1; $i < $y_length; ++$i) {
1158 for ($j = 0, $k = $i; $j < $x_length; ++$j, ++$k) {
1159 $temp = $product_value[$k] + $x_value[$j] * $y_value[$i] + $carry;
1160 $carry = (int) ($temp / 0x4000000);
1161 $product_value[$k] = (int) ($temp - 0x4000000 * $carry);
1164 $product_value[$k] = $carry;
1167 return $product_value;
1171 * Performs Karatsuba multiplication on two BigIntegers
1173 * See {@link http://en.wikipedia.org/wiki/Karatsuba_algorithm Karatsuba algorithm} and
1174 * {@link http://math.libtomcrypt.com/files/tommath.pdf#page=120 MPM 5.2.3}.
1176 * @param Array $x_value
1177 * @param Array $y_value
1181 function _karatsuba($x_value, $y_value)
1183 $m = min(count($x_value) >> 1, count($y_value) >> 1);
1185 if ($m < MATH_BIGINTEGER_KARATSUBA_CUTOFF) {
1186 return $this->_regularMultiply($x_value, $y_value);
1189 $x1 = array_slice($x_value, $m);
1190 $x0 = array_slice($x_value, 0, $m);
1191 $y1 = array_slice($y_value, $m);
1192 $y0 = array_slice($y_value, 0, $m);
1194 $z2 = $this->_karatsuba($x1, $y1);
1195 $z0 = $this->_karatsuba($x0, $y0);
1197 $z1 = $this->_add($x1, false, $x0, false);
1198 $temp = $this->_add($y1, false, $y0, false);
1199 $z1 = $this->_karatsuba($z1[MATH_BIGINTEGER_VALUE], $temp[MATH_BIGINTEGER_VALUE]);
1200 $temp = $this->_add($z2, false, $z0, false);
1201 $z1 = $this->_subtract($z1, false, $temp[MATH_BIGINTEGER_VALUE], false);
1203 $z2 = array_merge(array_fill(0, 2 * $m, 0), $z2);
1204 $z1[MATH_BIGINTEGER_VALUE] = array_merge(array_fill(0, $m, 0), $z1[MATH_BIGINTEGER_VALUE]);
1206 $xy = $this->_add($z2, false, $z1[MATH_BIGINTEGER_VALUE], $z1[MATH_BIGINTEGER_SIGN]);
1207 $xy = $this->_add($xy[MATH_BIGINTEGER_VALUE], $xy[MATH_BIGINTEGER_SIGN], $z0, false);
1209 return $xy[MATH_BIGINTEGER_VALUE];
1219 function _square($x = false)
1221 return count($x) < 2 * MATH_BIGINTEGER_KARATSUBA_CUTOFF ?
1222 $this->_trim($this->_baseSquare($x)) :
1223 $this->_trim($this->_karatsubaSquare($x));
1227 * Performs traditional squaring on two BigIntegers
1229 * Squaring can be done faster than multiplying a number by itself can be. See
1230 * {@link http://www.cacr.math.uwaterloo.ca/hac/about/chap14.pdf#page=7 HAC 14.2.4} /
1231 * {@link http://math.libtomcrypt.com/files/tommath.pdf#page=141 MPM 5.3} for more information.
1233 * @param Array $value
1237 function _baseSquare($value)
1239 if ( empty($value) ) {
1242 $square_value = $this->_array_repeat(0, 2 * count($value));
1244 for ($i = 0, $max_index = count($value) - 1; $i <= $max_index; ++$i) {
1247 $temp = $square_value[$i2] + $value[$i] * $value[$i];
1248 $carry = (int) ($temp / 0x4000000);
1249 $square_value[$i2] = (int) ($temp - 0x4000000 * $carry);
1251 // note how we start from $i+1 instead of 0 as we do in multiplication.
1252 for ($j = $i + 1, $k = $i2 + 1; $j <= $max_index; ++$j, ++$k) {
1253 $temp = $square_value[$k] + 2 * $value[$j] * $value[$i] + $carry;
1254 $carry = (int) ($temp / 0x4000000);
1255 $square_value[$k] = (int) ($temp - 0x4000000 * $carry);
1258 // the following line can yield values larger 2**15. at this point, PHP should switch
1260 $square_value[$i + $max_index + 1] = $carry;
1263 return $square_value;
1267 * Performs Karatsuba "squaring" on two BigIntegers
1269 * See {@link http://en.wikipedia.org/wiki/Karatsuba_algorithm Karatsuba algorithm} and
1270 * {@link http://math.libtomcrypt.com/files/tommath.pdf#page=151 MPM 5.3.4}.
1272 * @param Array $value
1276 function _karatsubaSquare($value)
1278 $m = count($value) >> 1;
1280 if ($m < MATH_BIGINTEGER_KARATSUBA_CUTOFF) {
1281 return $this->_baseSquare($value);
1284 $x1 = array_slice($value, $m);
1285 $x0 = array_slice($value, 0, $m);
1287 $z2 = $this->_karatsubaSquare($x1);
1288 $z0 = $this->_karatsubaSquare($x0);
1290 $z1 = $this->_add($x1, false, $x0, false);
1291 $z1 = $this->_karatsubaSquare($z1[MATH_BIGINTEGER_VALUE]);
1292 $temp = $this->_add($z2, false, $z0, false);
1293 $z1 = $this->_subtract($z1, false, $temp[MATH_BIGINTEGER_VALUE], false);
1295 $z2 = array_merge(array_fill(0, 2 * $m, 0), $z2);
1296 $z1[MATH_BIGINTEGER_VALUE] = array_merge(array_fill(0, $m, 0), $z1[MATH_BIGINTEGER_VALUE]);
1298 $xx = $this->_add($z2, false, $z1[MATH_BIGINTEGER_VALUE], $z1[MATH_BIGINTEGER_SIGN]);
1299 $xx = $this->_add($xx[MATH_BIGINTEGER_VALUE], $xx[MATH_BIGINTEGER_SIGN], $z0, false);
1301 return $xx[MATH_BIGINTEGER_VALUE];
1305 * Divides two BigIntegers.
1307 * Returns an array whose first element contains the quotient and whose second element contains the
1308 * "common residue". If the remainder would be positive, the "common residue" and the remainder are the
1309 * same. If the remainder would be negative, the "common residue" is equal to the sum of the remainder
1310 * and the divisor (basically, the "common residue" is the first positive modulo).
1312 * Here's an example:
1315 * include('Math/BigInteger.php');
1317 * $a = new Math_BigInteger('10');
1318 * $b = new Math_BigInteger('20');
1320 * list($quotient, $remainder) = $a->divide($b);
1322 * echo $quotient->toString(); // outputs 0
1324 * echo $remainder->toString(); // outputs 10
1328 * @param Math_BigInteger $y
1331 * @internal This function is based off of {@link http://www.cacr.math.uwaterloo.ca/hac/about/chap14.pdf#page=9 HAC 14.20}.
1335 switch ( MATH_BIGINTEGER_MODE ) {
1336 case MATH_BIGINTEGER_MODE_GMP:
1337 $quotient = new Math_BigInteger();
1338 $remainder = new Math_BigInteger();
1340 list($quotient->value, $remainder->value) = gmp_div_qr($this->value, $y->value);
1342 if (gmp_sign($remainder->value) < 0) {
1343 $remainder->value = gmp_add($remainder->value, gmp_abs($y->value));
1346 return array($this->_normalize($quotient), $this->_normalize($remainder));
1347 case MATH_BIGINTEGER_MODE_BCMATH:
1348 $quotient = new Math_BigInteger();
1349 $remainder = new Math_BigInteger();
1351 $quotient->value = bcdiv($this->value, $y->value, 0);
1352 $remainder->value = bcmod($this->value, $y->value);
1354 if ($remainder->value[0] == '-') {
1355 $remainder->value = bcadd($remainder->value, $y->value[0] == '-' ? substr($y->value, 1) : $y->value, 0);
1358 return array($this->_normalize($quotient), $this->_normalize($remainder));
1361 if (count($y->value) == 1) {
1362 list($q, $r) = $this->_divide_digit($this->value, $y->value[0]);
1363 $quotient = new Math_BigInteger();
1364 $remainder = new Math_BigInteger();
1365 $quotient->value = $q;
1366 $remainder->value = array($r);
1367 $quotient->is_negative = $this->is_negative != $y->is_negative;
1368 return array($this->_normalize($quotient), $this->_normalize($remainder));
1372 if ( !isset($zero) ) {
1373 $zero = new Math_BigInteger();
1379 $x_sign = $x->is_negative;
1380 $y_sign = $y->is_negative;
1382 $x->is_negative = $y->is_negative = false;
1384 $diff = $x->compare($y);
1387 $temp = new Math_BigInteger();
1388 $temp->value = array(1);
1389 $temp->is_negative = $x_sign != $y_sign;
1390 return array($this->_normalize($temp), $this->_normalize(new Math_BigInteger()));
1394 // if $x is negative, "add" $y.
1396 $x = $y->subtract($x);
1398 return array($this->_normalize(new Math_BigInteger()), $this->_normalize($x));
1401 // normalize $x and $y as described in HAC 14.23 / 14.24
1402 $msb = $y->value[count($y->value) - 1];
1403 for ($shift = 0; !($msb & 0x2000000); ++$shift) {
1406 $x->_lshift($shift);
1407 $y->_lshift($shift);
1408 $y_value = &$y->value;
1410 $x_max = count($x->value) - 1;
1411 $y_max = count($y->value) - 1;
1413 $quotient = new Math_BigInteger();
1414 $quotient_value = &$quotient->value;
1415 $quotient_value = $this->_array_repeat(0, $x_max - $y_max + 1);
1417 static $temp, $lhs, $rhs;
1418 if (!isset($temp)) {
1419 $temp = new Math_BigInteger();
1420 $lhs = new Math_BigInteger();
1421 $rhs = new Math_BigInteger();
1423 $temp_value = &$temp->value;
1424 $rhs_value = &$rhs->value;
1426 // $temp = $y << ($x_max - $y_max-1) in base 2**26
1427 $temp_value = array_merge($this->_array_repeat(0, $x_max - $y_max), $y_value);
1429 while ( $x->compare($temp) >= 0 ) {
1430 // calculate the "common residue"
1431 ++$quotient_value[$x_max - $y_max];
1432 $x = $x->subtract($temp);
1433 $x_max = count($x->value) - 1;
1436 for ($i = $x_max; $i >= $y_max + 1; --$i) {
1437 $x_value = &$x->value;
1439 isset($x_value[$i]) ? $x_value[$i] : 0,
1440 isset($x_value[$i - 1]) ? $x_value[$i - 1] : 0,
1441 isset($x_value[$i - 2]) ? $x_value[$i - 2] : 0
1445 ( $y_max > 0 ) ? $y_value[$y_max - 1] : 0
1448 $q_index = $i - $y_max - 1;
1449 if ($x_window[0] == $y_window[0]) {
1450 $quotient_value[$q_index] = 0x3FFFFFF;
1452 $quotient_value[$q_index] = (int) (
1453 ($x_window[0] * 0x4000000 + $x_window[1])
1459 $temp_value = array($y_window[1], $y_window[0]);
1461 $lhs->value = array($quotient_value[$q_index]);
1462 $lhs = $lhs->multiply($temp);
1464 $rhs_value = array($x_window[2], $x_window[1], $x_window[0]);
1466 while ( $lhs->compare($rhs) > 0 ) {
1467 --$quotient_value[$q_index];
1469 $lhs->value = array($quotient_value[$q_index]);
1470 $lhs = $lhs->multiply($temp);
1473 $adjust = $this->_array_repeat(0, $q_index);
1474 $temp_value = array($quotient_value[$q_index]);
1475 $temp = $temp->multiply($y);
1476 $temp_value = &$temp->value;
1477 $temp_value = array_merge($adjust, $temp_value);
1479 $x = $x->subtract($temp);
1481 if ($x->compare($zero) < 0) {
1482 $temp_value = array_merge($adjust, $y_value);
1483 $x = $x->add($temp);
1485 --$quotient_value[$q_index];
1488 $x_max = count($x_value) - 1;
1491 // unnormalize the remainder
1492 $x->_rshift($shift);
1494 $quotient->is_negative = $x_sign != $y_sign;
1496 // calculate the "common residue", if appropriate
1498 $y->_rshift($shift);
1499 $x = $y->subtract($x);
1502 return array($this->_normalize($quotient), $this->_normalize($x));
1506 * Divides a BigInteger by a regular integer
1508 * abc / x = a00 / x + b0 / x + c / x
1510 * @param Array $dividend
1511 * @param Array $divisor
1515 function _divide_digit($dividend, $divisor)
1520 for ($i = count($dividend) - 1; $i >= 0; --$i) {
1521 $temp = 0x4000000 * $carry + $dividend[$i];
1522 $result[$i] = (int) ($temp / $divisor);
1523 $carry = (int) ($temp - $divisor * $result[$i]);
1526 return array($result, $carry);
1530 * Performs modular exponentiation.
1532 * Here's an example:
1535 * include('Math/BigInteger.php');
1537 * $a = new Math_BigInteger('10');
1538 * $b = new Math_BigInteger('20');
1539 * $c = new Math_BigInteger('30');
1541 * $c = $a->modPow($b, $c);
1543 * echo $c->toString(); // outputs 10
1547 * @param Math_BigInteger $e
1548 * @param Math_BigInteger $n
1549 * @return Math_BigInteger
1551 * @internal The most naive approach to modular exponentiation has very unreasonable requirements, and
1552 * and although the approach involving repeated squaring does vastly better, it, too, is impractical
1553 * for our purposes. The reason being that division - by far the most complicated and time-consuming
1554 * of the basic operations (eg. +,-,*,/) - occurs multiple times within it.
1556 * Modular reductions resolve this issue. Although an individual modular reduction takes more time
1557 * then an individual division, when performed in succession (with the same modulo), they're a lot faster.
1559 * The two most commonly used modular reductions are Barrett and Montgomery reduction. Montgomery reduction,
1560 * although faster, only works when the gcd of the modulo and of the base being used is 1. In RSA, when the
1561 * base is a power of two, the modulo - a product of two primes - is always going to have a gcd of 1 (because
1562 * the product of two odd numbers is odd), but what about when RSA isn't used?
1564 * In contrast, Barrett reduction has no such constraint. As such, some bigint implementations perform a
1565 * Barrett reduction after every operation in the modpow function. Others perform Barrett reductions when the
1566 * modulo is even and Montgomery reductions when the modulo is odd. BigInteger.java's modPow method, however,
1567 * uses a trick involving the Chinese Remainder Theorem to factor the even modulo into two numbers - one odd and
1568 * the other, a power of two - and recombine them, later. This is the method that this modPow function uses.
1569 * {@link http://islab.oregonstate.edu/papers/j34monex.pdf Montgomery Reduction with Even Modulus} elaborates.
1571 function modPow($e, $n)
1573 $n = $this->bitmask !== false && $this->bitmask->compare($n) < 0 ? $this->bitmask : $n->abs();
1575 if ($e->compare(new Math_BigInteger()) < 0) {
1578 $temp = $this->modInverse($n);
1579 if ($temp === false) {
1583 return $this->_normalize($temp->modPow($e, $n));
1586 switch ( MATH_BIGINTEGER_MODE ) {
1587 case MATH_BIGINTEGER_MODE_GMP:
1588 $temp = new Math_BigInteger();
1589 $temp->value = gmp_powm($this->value, $e->value, $n->value);
1591 return $this->_normalize($temp);
1592 case MATH_BIGINTEGER_MODE_BCMATH:
1593 $temp = new Math_BigInteger();
1594 $temp->value = bcpowmod($this->value, $e->value, $n->value, 0);
1596 return $this->_normalize($temp);
1599 if ( empty($e->value) ) {
1600 $temp = new Math_BigInteger();
1601 $temp->value = array(1);
1602 return $this->_normalize($temp);
1605 if ( $e->value == array(1) ) {
1606 list(, $temp) = $this->divide($n);
1607 return $this->_normalize($temp);
1610 if ( $e->value == array(2) ) {
1611 $temp = new Math_BigInteger();
1612 $temp->value = $this->_square($this->value);
1613 list(, $temp) = $temp->divide($n);
1614 return $this->_normalize($temp);
1617 return $this->_normalize($this->_slidingWindow($e, $n, MATH_BIGINTEGER_BARRETT));
1619 // is the modulo odd?
1620 if ( $n->value[0] & 1 ) {
1621 return $this->_normalize($this->_slidingWindow($e, $n, MATH_BIGINTEGER_MONTGOMERY));
1623 // if it's not, it's even
1625 // find the lowest set bit (eg. the max pow of 2 that divides $n)
1626 for ($i = 0; $i < count($n->value); ++$i) {
1627 if ( $n->value[$i] ) {
1628 $temp = decbin($n->value[$i]);
1629 $j = strlen($temp) - strrpos($temp, '1') - 1;
1634 // at this point, 2^$j * $n/(2^$j) == $n
1638 $mod2 = new Math_BigInteger();
1639 $mod2->value = array(1);
1642 $part1 = ( $mod1->value != array(1) ) ? $this->_slidingWindow($e, $mod1, MATH_BIGINTEGER_MONTGOMERY) : new Math_BigInteger();
1643 $part2 = $this->_slidingWindow($e, $mod2, MATH_BIGINTEGER_POWEROF2);
1645 $y1 = $mod2->modInverse($mod1);
1646 $y2 = $mod1->modInverse($mod2);
1648 $result = $part1->multiply($mod2);
1649 $result = $result->multiply($y1);
1651 $temp = $part2->multiply($mod1);
1652 $temp = $temp->multiply($y2);
1654 $result = $result->add($temp);
1655 list(, $result) = $result->divide($n);
1657 return $this->_normalize($result);
1661 * Performs modular exponentiation.
1663 * Alias for Math_BigInteger::modPow()
1665 * @param Math_BigInteger $e
1666 * @param Math_BigInteger $n
1667 * @return Math_BigInteger
1670 function powMod($e, $n)
1672 return $this->modPow($e, $n);
1676 * Sliding Window k-ary Modular Exponentiation
1678 * Based on {@link http://www.cacr.math.uwaterloo.ca/hac/about/chap14.pdf#page=27 HAC 14.85} /
1679 * {@link http://math.libtomcrypt.com/files/tommath.pdf#page=210 MPM 7.7}. In a departure from those algorithims,
1680 * however, this function performs a modular reduction after every multiplication and squaring operation.
1681 * As such, this function has the same preconditions that the reductions being used do.
1683 * @param Math_BigInteger $e
1684 * @param Math_BigInteger $n
1685 * @param Integer $mode
1686 * @return Math_BigInteger
1689 function _slidingWindow($e, $n, $mode)
1691 static $window_ranges = array(7, 25, 81, 241, 673, 1793); // from BigInteger.java's oddModPow function
1692 //static $window_ranges = array(0, 7, 36, 140, 450, 1303, 3529); // from MPM 7.3.1
1694 $e_value = $e->value;
1695 $e_length = count($e_value) - 1;
1696 $e_bits = decbin($e_value[$e_length]);
1697 for ($i = $e_length - 1; $i >= 0; --$i) {
1698 $e_bits.= str_pad(decbin($e_value[$i]), 26, '0', STR_PAD_LEFT);
1701 $e_length = strlen($e_bits);
1703 // calculate the appropriate window size.
1704 // $window_size == 3 if $window_ranges is between 25 and 81, for example.
1705 for ($i = 0, $window_size = 1; $e_length > $window_ranges[$i] && $i < count($window_ranges); ++$window_size, ++$i);
1707 $n_value = $n->value;
1709 // precompute $this^0 through $this^$window_size
1711 $powers[1] = $this->_prepareReduce($this->value, $n_value, $mode);
1712 $powers[2] = $this->_squareReduce($powers[1], $n_value, $mode);
1714 // we do every other number since substr($e_bits, $i, $j+1) (see below) is supposed to end
1715 // in a 1. ie. it's supposed to be odd.
1716 $temp = 1 << ($window_size - 1);
1717 for ($i = 1; $i < $temp; ++$i) {
1719 $powers[$i2 + 1] = $this->_multiplyReduce($powers[$i2 - 1], $powers[2], $n_value, $mode);
1723 $result = $this->_prepareReduce($result, $n_value, $mode);
1725 for ($i = 0; $i < $e_length; ) {
1726 if ( !$e_bits[$i] ) {
1727 $result = $this->_squareReduce($result, $n_value, $mode);
1730 for ($j = $window_size - 1; $j > 0; --$j) {
1731 if ( !empty($e_bits[$i + $j]) ) {
1736 for ($k = 0; $k <= $j; ++$k) {// eg. the length of substr($e_bits, $i, $j+1)
1737 $result = $this->_squareReduce($result, $n_value, $mode);
1740 $result = $this->_multiplyReduce($result, $powers[bindec(substr($e_bits, $i, $j + 1))], $n_value, $mode);
1746 $temp = new Math_BigInteger();
1747 $temp->value = $this->_reduce($result, $n_value, $mode);
1755 * For most $modes this will return the remainder.
1757 * @see _slidingWindow()
1761 * @param Integer $mode
1764 function _reduce($x, $n, $mode)
1767 case MATH_BIGINTEGER_MONTGOMERY:
1768 return $this->_montgomery($x, $n);
1769 case MATH_BIGINTEGER_BARRETT:
1770 return $this->_barrett($x, $n);
1771 case MATH_BIGINTEGER_POWEROF2:
1772 $lhs = new Math_BigInteger();
1774 $rhs = new Math_BigInteger();
1776 return $x->_mod2($n);
1777 case MATH_BIGINTEGER_CLASSIC:
1778 $lhs = new Math_BigInteger();
1780 $rhs = new Math_BigInteger();
1782 list(, $temp) = $lhs->divide($rhs);
1783 return $temp->value;
1784 case MATH_BIGINTEGER_NONE:
1787 // an invalid $mode was provided
1792 * Modular reduction preperation
1794 * @see _slidingWindow()
1798 * @param Integer $mode
1801 function _prepareReduce($x, $n, $mode)
1803 if ($mode == MATH_BIGINTEGER_MONTGOMERY) {
1804 return $this->_prepMontgomery($x, $n);
1806 return $this->_reduce($x, $n, $mode);
1812 * @see _slidingWindow()
1817 * @param Integer $mode
1820 function _multiplyReduce($x, $y, $n, $mode)
1822 if ($mode == MATH_BIGINTEGER_MONTGOMERY) {
1823 return $this->_montgomeryMultiply($x, $y, $n);
1825 $temp = $this->_multiply($x, false, $y, false);
1826 return $this->_reduce($temp[MATH_BIGINTEGER_VALUE], $n, $mode);
1832 * @see _slidingWindow()
1836 * @param Integer $mode
1839 function _squareReduce($x, $n, $mode)
1841 if ($mode == MATH_BIGINTEGER_MONTGOMERY) {
1842 return $this->_montgomeryMultiply($x, $x, $n);
1844 return $this->_reduce($this->_square($x), $n, $mode);
1848 * Modulos for Powers of Two
1850 * Calculates $x%$n, where $n = 2**$e, for some $e. Since this is basically the same as doing $x & ($n-1),
1851 * we'll just use this function as a wrapper for doing that.
1853 * @see _slidingWindow()
1855 * @param Math_BigInteger
1856 * @return Math_BigInteger
1860 $temp = new Math_BigInteger();
1861 $temp->value = array(1);
1862 return $this->bitwise_and($n->subtract($temp));
1866 * Barrett Modular Reduction
1868 * See {@link http://www.cacr.math.uwaterloo.ca/hac/about/chap14.pdf#page=14 HAC 14.3.3} /
1869 * {@link http://math.libtomcrypt.com/files/tommath.pdf#page=165 MPM 6.2.5} for more information. Modified slightly,
1870 * so as not to require negative numbers (initially, this script didn't support negative numbers).
1872 * Employs "folding", as described at
1873 * {@link http://www.cosic.esat.kuleuven.be/publications/thesis-149.pdf#page=66 thesis-149.pdf#page=66}. To quote from
1874 * it, "the idea [behind folding] is to find a value x' such that x (mod m) = x' (mod m), with x' being smaller than x."
1876 * Unfortunately, the "Barrett Reduction with Folding" algorithm described in thesis-149.pdf is not, as written, all that
1877 * usable on account of (1) its not using reasonable radix points as discussed in
1878 * {@link http://math.libtomcrypt.com/files/tommath.pdf#page=162 MPM 6.2.2} and (2) the fact that, even with reasonable
1879 * radix points, it only works when there are an even number of digits in the denominator. The reason for (2) is that
1880 * (x >> 1) + (x >> 1) != x / 2 + x / 2. If x is even, they're the same, but if x is odd, they're not. See the in-line
1881 * comments for details.
1883 * @see _slidingWindow()
1889 function _barrett($n, $m)
1891 static $cache = array(
1892 MATH_BIGINTEGER_VARIABLE => array(),
1893 MATH_BIGINTEGER_DATA => array()
1896 $m_length = count($m);
1898 // if ($this->_compare($n, $this->_square($m)) >= 0) {
1899 if (count($n) > 2 * $m_length) {
1900 $lhs = new Math_BigInteger();
1901 $rhs = new Math_BigInteger();
1904 list(, $temp) = $lhs->divide($rhs);
1905 return $temp->value;
1908 // if (m.length >> 1) + 2 <= m.length then m is too small and n can't be reduced
1909 if ($m_length < 5) {
1910 return $this->_regularBarrett($n, $m);
1915 if ( ($key = array_search($m, $cache[MATH_BIGINTEGER_VARIABLE])) === false ) {
1916 $key = count($cache[MATH_BIGINTEGER_VARIABLE]);
1917 $cache[MATH_BIGINTEGER_VARIABLE][] = $m;
1919 $lhs = new Math_BigInteger();
1920 $lhs_value = &$lhs->value;
1921 $lhs_value = $this->_array_repeat(0, $m_length + ($m_length >> 1));
1923 $rhs = new Math_BigInteger();
1926 list($u, $m1) = $lhs->divide($rhs);
1930 $cache[MATH_BIGINTEGER_DATA][] = array(
1931 'u' => $u, // m.length >> 1 (technically (m.length >> 1) + 1)
1932 'm1'=> $m1 // m.length
1935 extract($cache[MATH_BIGINTEGER_DATA][$key]);
1938 $cutoff = $m_length + ($m_length >> 1);
1939 $lsd = array_slice($n, 0, $cutoff); // m.length + (m.length >> 1)
1940 $msd = array_slice($n, $cutoff); // m.length >> 1
1941 $lsd = $this->_trim($lsd);
1942 $temp = $this->_multiply($msd, false, $m1, false);
1943 $n = $this->_add($lsd, false, $temp[MATH_BIGINTEGER_VALUE], false); // m.length + (m.length >> 1) + 1
1945 if ($m_length & 1) {
1946 return $this->_regularBarrett($n[MATH_BIGINTEGER_VALUE], $m);
1949 // (m.length + (m.length >> 1) + 1) - (m.length - 1) == (m.length >> 1) + 2
1950 $temp = array_slice($n[MATH_BIGINTEGER_VALUE], $m_length - 1);
1951 // if even: ((m.length >> 1) + 2) + (m.length >> 1) == m.length + 2
1952 // if odd: ((m.length >> 1) + 2) + (m.length >> 1) == (m.length - 1) + 2 == m.length + 1
1953 $temp = $this->_multiply($temp, false, $u, false);
1954 // if even: (m.length + 2) - ((m.length >> 1) + 1) = m.length - (m.length >> 1) + 1
1955 // if odd: (m.length + 1) - ((m.length >> 1) + 1) = m.length - (m.length >> 1)
1956 $temp = array_slice($temp[MATH_BIGINTEGER_VALUE], ($m_length >> 1) + 1);
1957 // if even: (m.length - (m.length >> 1) + 1) + m.length = 2 * m.length - (m.length >> 1) + 1
1958 // if odd: (m.length - (m.length >> 1)) + m.length = 2 * m.length - (m.length >> 1)
1959 $temp = $this->_multiply($temp, false, $m, false);
1961 // at this point, if m had an odd number of digits, we'd be subtracting a 2 * m.length - (m.length >> 1) digit
1962 // number from a m.length + (m.length >> 1) + 1 digit number. ie. there'd be an extra digit and the while loop
1963 // following this comment would loop a lot (hence our calling _regularBarrett() in that situation).
1965 $result = $this->_subtract($n[MATH_BIGINTEGER_VALUE], false, $temp[MATH_BIGINTEGER_VALUE], false);
1967 while ($this->_compare($result[MATH_BIGINTEGER_VALUE], $result[MATH_BIGINTEGER_SIGN], $m, false) >= 0) {
1968 $result = $this->_subtract($result[MATH_BIGINTEGER_VALUE], $result[MATH_BIGINTEGER_SIGN], $m, false);
1971 return $result[MATH_BIGINTEGER_VALUE];
1975 * (Regular) Barrett Modular Reduction
1977 * For numbers with more than four digits Math_BigInteger::_barrett() is faster. The difference between that and this
1978 * is that this function does not fold the denominator into a smaller form.
1980 * @see _slidingWindow()
1986 function _regularBarrett($x, $n)
1988 static $cache = array(
1989 MATH_BIGINTEGER_VARIABLE => array(),
1990 MATH_BIGINTEGER_DATA => array()
1993 $n_length = count($n);
1995 if (count($x) > 2 * $n_length) {
1996 $lhs = new Math_BigInteger();
1997 $rhs = new Math_BigInteger();
2000 list(, $temp) = $lhs->divide($rhs);
2001 return $temp->value;
2004 if ( ($key = array_search($n, $cache[MATH_BIGINTEGER_VARIABLE])) === false ) {
2005 $key = count($cache[MATH_BIGINTEGER_VARIABLE]);
2006 $cache[MATH_BIGINTEGER_VARIABLE][] = $n;
2007 $lhs = new Math_BigInteger();
2008 $lhs_value = &$lhs->value;
2009 $lhs_value = $this->_array_repeat(0, 2 * $n_length);
2011 $rhs = new Math_BigInteger();
2013 list($temp, ) = $lhs->divide($rhs); // m.length
2014 $cache[MATH_BIGINTEGER_DATA][] = $temp->value;
2017 // 2 * m.length - (m.length - 1) = m.length + 1
2018 $temp = array_slice($x, $n_length - 1);
2019 // (m.length + 1) + m.length = 2 * m.length + 1
2020 $temp = $this->_multiply($temp, false, $cache[MATH_BIGINTEGER_DATA][$key], false);
2021 // (2 * m.length + 1) - (m.length - 1) = m.length + 2
2022 $temp = array_slice($temp[MATH_BIGINTEGER_VALUE], $n_length + 1);
2025 $result = array_slice($x, 0, $n_length + 1);
2027 $temp = $this->_multiplyLower($temp, false, $n, false, $n_length + 1);
2028 // $temp == array_slice($temp->_multiply($temp, false, $n, false)->value, 0, $n_length + 1)
2030 if ($this->_compare($result, false, $temp[MATH_BIGINTEGER_VALUE], $temp[MATH_BIGINTEGER_SIGN]) < 0) {
2031 $corrector_value = $this->_array_repeat(0, $n_length + 1);
2032 $corrector_value[] = 1;
2033 $result = $this->_add($result, false, $corrector, false);
2034 $result = $result[MATH_BIGINTEGER_VALUE];
2037 // at this point, we're subtracting a number with m.length + 1 digits from another number with m.length + 1 digits
2038 $result = $this->_subtract($result, false, $temp[MATH_BIGINTEGER_VALUE], $temp[MATH_BIGINTEGER_SIGN]);
2039 while ($this->_compare($result[MATH_BIGINTEGER_VALUE], $result[MATH_BIGINTEGER_SIGN], $n, false) > 0) {
2040 $result = $this->_subtract($result[MATH_BIGINTEGER_VALUE], $result[MATH_BIGINTEGER_SIGN], $n, false);
2043 return $result[MATH_BIGINTEGER_VALUE];
2047 * Performs long multiplication up to $stop digits
2049 * If you're going to be doing array_slice($product->value, 0, $stop), some cycles can be saved.
2051 * @see _regularBarrett()
2052 * @param Array $x_value
2053 * @param Boolean $x_negative
2054 * @param Array $y_value
2055 * @param Boolean $y_negative
2059 function _multiplyLower($x_value, $x_negative, $y_value, $y_negative, $stop)
2061 $x_length = count($x_value);
2062 $y_length = count($y_value);
2064 if ( !$x_length || !$y_length ) { // a 0 is being multiplied
2066 MATH_BIGINTEGER_VALUE => array(),
2067 MATH_BIGINTEGER_SIGN => false
2071 if ( $x_length < $y_length ) {
2073 $x_value = $y_value;
2076 $x_length = count($x_value);
2077 $y_length = count($y_value);
2080 $product_value = $this->_array_repeat(0, $x_length + $y_length);
2082 // the following for loop could be removed if the for loop following it
2083 // (the one with nested for loops) initially set $i to 0, but
2084 // doing so would also make the result in one set of unnecessary adds,
2085 // since on the outermost loops first pass, $product->value[$k] is going
2090 for ($j = 0; $j < $x_length; ++$j) { // ie. $i = 0, $k = $i
2091 $temp = $x_value[$j] * $y_value[0] + $carry; // $product_value[$k] == 0
2092 $carry = (int) ($temp / 0x4000000);
2093 $product_value[$j] = (int) ($temp - 0x4000000 * $carry);
2097 $product_value[$j] = $carry;
2100 // the above for loop is what the previous comment was talking about. the
2101 // following for loop is the "one with nested for loops"
2103 for ($i = 1; $i < $y_length; ++$i) {
2106 for ($j = 0, $k = $i; $j < $x_length && $k < $stop; ++$j, ++$k) {
2107 $temp = $product_value[$k] + $x_value[$j] * $y_value[$i] + $carry;
2108 $carry = (int) ($temp / 0x4000000);
2109 $product_value[$k] = (int) ($temp - 0x4000000 * $carry);
2113 $product_value[$k] = $carry;
2118 MATH_BIGINTEGER_VALUE => $this->_trim($product_value),
2119 MATH_BIGINTEGER_SIGN => $x_negative != $y_negative
2124 * Montgomery Modular Reduction
2126 * ($x->_prepMontgomery($n))->_montgomery($n) yields $x % $n.
2127 * {@link http://math.libtomcrypt.com/files/tommath.pdf#page=170 MPM 6.3} provides insights on how this can be
2128 * improved upon (basically, by using the comba method). gcd($n, 2) must be equal to one for this function
2129 * to work correctly.
2131 * @see _prepMontgomery()
2132 * @see _slidingWindow()
2138 function _montgomery($x, $n)
2140 static $cache = array(
2141 MATH_BIGINTEGER_VARIABLE => array(),
2142 MATH_BIGINTEGER_DATA => array()
2145 if ( ($key = array_search($n, $cache[MATH_BIGINTEGER_VARIABLE])) === false ) {
2146 $key = count($cache[MATH_BIGINTEGER_VARIABLE]);
2147 $cache[MATH_BIGINTEGER_VARIABLE][] = $x;
2148 $cache[MATH_BIGINTEGER_DATA][] = $this->_modInverse67108864($n);
2153 $result = array(MATH_BIGINTEGER_VALUE => $x);
2155 for ($i = 0; $i < $k; ++$i) {
2156 $temp = $result[MATH_BIGINTEGER_VALUE][$i] * $cache[MATH_BIGINTEGER_DATA][$key];
2157 $temp = (int) ($temp - 0x4000000 * ((int) ($temp / 0x4000000)));
2158 $temp = $this->_regularMultiply(array($temp), $n);
2159 $temp = array_merge($this->_array_repeat(0, $i), $temp);
2160 $result = $this->_add($result[MATH_BIGINTEGER_VALUE], false, $temp, false);
2163 $result[MATH_BIGINTEGER_VALUE] = array_slice($result[MATH_BIGINTEGER_VALUE], $k);
2165 if ($this->_compare($result, false, $n, false) >= 0) {
2166 $result = $this->_subtract($result[MATH_BIGINTEGER_VALUE], false, $n, false);
2169 return $result[MATH_BIGINTEGER_VALUE];
2173 * Montgomery Multiply
2175 * Interleaves the montgomery reduction and long multiplication algorithms together as described in
2176 * {@link http://www.cacr.math.uwaterloo.ca/hac/about/chap14.pdf#page=13 HAC 14.36}
2178 * @see _prepMontgomery()
2179 * @see _montgomery()
2186 function _montgomeryMultiply($x, $y, $m)
2188 $temp = $this->_multiply($x, false, $y, false);
2189 return $this->_montgomery($temp[MATH_BIGINTEGER_VALUE], $m);
2191 static $cache = array(
2192 MATH_BIGINTEGER_VARIABLE => array(),
2193 MATH_BIGINTEGER_DATA => array()
2196 if ( ($key = array_search($m, $cache[MATH_BIGINTEGER_VARIABLE])) === false ) {
2197 $key = count($cache[MATH_BIGINTEGER_VARIABLE]);
2198 $cache[MATH_BIGINTEGER_VARIABLE][] = $m;
2199 $cache[MATH_BIGINTEGER_DATA][] = $this->_modInverse67108864($m);
2202 $n = max(count($x), count($y), count($m));
2203 $x = array_pad($x, $n, 0);
2204 $y = array_pad($y, $n, 0);
2205 $m = array_pad($m, $n, 0);
2206 $a = array(MATH_BIGINTEGER_VALUE => $this->_array_repeat(0, $n + 1));
2207 for ($i = 0; $i < $n; ++$i) {
2208 $temp = $a[MATH_BIGINTEGER_VALUE][0] + $x[$i] * $y[0];
2209 $temp = (int) ($temp - 0x4000000 * ((int) ($temp / 0x4000000)));
2210 $temp = $temp * $cache[MATH_BIGINTEGER_DATA][$key];
2211 $temp = (int) ($temp - 0x4000000 * ((int) ($temp / 0x4000000)));
2212 $temp = $this->_add($this->_regularMultiply(array($x[$i]), $y), false, $this->_regularMultiply(array($temp), $m), false);
2213 $a = $this->_add($a[MATH_BIGINTEGER_VALUE], false, $temp[MATH_BIGINTEGER_VALUE], false);
2214 $a[MATH_BIGINTEGER_VALUE] = array_slice($a[MATH_BIGINTEGER_VALUE], 1);
2216 if ($this->_compare($a[MATH_BIGINTEGER_VALUE], false, $m, false) >= 0) {
2217 $a = $this->_subtract($a[MATH_BIGINTEGER_VALUE], false, $m, false);
2219 return $a[MATH_BIGINTEGER_VALUE];
2223 * Prepare a number for use in Montgomery Modular Reductions
2225 * @see _montgomery()
2226 * @see _slidingWindow()
2232 function _prepMontgomery($x, $n)
2234 $lhs = new Math_BigInteger();
2235 $lhs->value = array_merge($this->_array_repeat(0, count($n)), $x);
2236 $rhs = new Math_BigInteger();
2239 list(, $temp) = $lhs->divide($rhs);
2240 return $temp->value;
2244 * Modular Inverse of a number mod 2**26 (eg. 67108864)
2246 * Based off of the bnpInvDigit function implemented and justified in the following URL:
2248 * {@link http://www-cs-students.stanford.edu/~tjw/jsbn/jsbn.js}
2250 * The following URL provides more info:
2252 * {@link http://groups.google.com/group/sci.crypt/msg/7a137205c1be7d85}
2254 * As for why we do all the bitmasking... strange things can happen when converting from floats to ints. For
2255 * instance, on some computers, var_dump((int) -4294967297) yields int(-1) and on others, it yields
2256 * int(-2147483648). To avoid problems stemming from this, we use bitmasks to guarantee that ints aren't
2257 * auto-converted to floats. The outermost bitmask is present because without it, there's no guarantee that
2258 * the "residue" returned would be the so-called "common residue". We use fmod, in the last step, because the
2259 * maximum possible $x is 26 bits and the maximum $result is 16 bits. Thus, we have to be able to handle up to
2260 * 40 bits, which only 64-bit floating points will support.
2262 * Thanks to Pedro Gimeno Fortea for input!
2264 * @see _montgomery()
2269 function _modInverse67108864($x) // 2**26 == 67108864
2272 $result = $x & 0x3; // x**-1 mod 2**2
2273 $result = ($result * (2 - $x * $result)) & 0xF; // x**-1 mod 2**4
2274 $result = ($result * (2 - ($x & 0xFF) * $result)) & 0xFF; // x**-1 mod 2**8
2275 $result = ($result * ((2 - ($x & 0xFFFF) * $result) & 0xFFFF)) & 0xFFFF; // x**-1 mod 2**16
2276 $result = fmod($result * (2 - fmod($x * $result, 0x4000000)), 0x4000000); // x**-1 mod 2**26
2277 return $result & 0x3FFFFFF;
2281 * Calculates modular inverses.
2283 * Say you have (30 mod 17 * x mod 17) mod 17 == 1. x can be found using modular inverses.
2285 * Here's an example:
2288 * include('Math/BigInteger.php');
2290 * $a = new Math_BigInteger(30);
2291 * $b = new Math_BigInteger(17);
2293 * $c = $a->modInverse($b);
2294 * echo $c->toString(); // outputs 4
2298 * $d = $a->multiply($c);
2299 * list(, $d) = $d->divide($b);
2300 * echo $d; // outputs 1 (as per the definition of modular inverse)
2304 * @param Math_BigInteger $n
2305 * @return mixed false, if no modular inverse exists, Math_BigInteger, otherwise.
2307 * @internal See {@link http://www.cacr.math.uwaterloo.ca/hac/about/chap14.pdf#page=21 HAC 14.64} for more information.
2309 function modInverse($n)
2311 switch ( MATH_BIGINTEGER_MODE ) {
2312 case MATH_BIGINTEGER_MODE_GMP:
2313 $temp = new Math_BigInteger();
2314 $temp->value = gmp_invert($this->value, $n->value);
2316 return ( $temp->value === false ) ? false : $this->_normalize($temp);
2320 if (!isset($zero)) {
2321 $zero = new Math_BigInteger();
2322 $one = new Math_BigInteger(1);
2325 // $x mod $n == $x mod -$n.
2328 if ($this->compare($zero) < 0) {
2329 $temp = $this->abs();
2330 $temp = $temp->modInverse($n);
2331 return $negated === false ? false : $this->_normalize($n->subtract($temp));
2334 extract($this->extendedGCD($n));
2336 if (!$gcd->equals($one)) {
2340 $x = $x->compare($zero) < 0 ? $x->add($n) : $x;
2342 return $this->compare($zero) < 0 ? $this->_normalize($n->subtract($x)) : $this->_normalize($x);
2346 * Calculates the greatest common divisor and Bézout's identity.
2348 * Say you have 693 and 609. The GCD is 21. Bézout's identity states that there exist integers x and y such that
2349 * 693*x + 609*y == 21. In point of fact, there are actually an infinite number of x and y combinations and which
2350 * combination is returned is dependant upon which mode is in use. See
2351 * {@link http://en.wikipedia.org/wiki/B%C3%A9zout%27s_identity Bézout's identity - Wikipedia} for more information.
2353 * Here's an example:
2356 * include('Math/BigInteger.php');
2358 * $a = new Math_BigInteger(693);
2359 * $b = new Math_BigInteger(609);
2361 * extract($a->extendedGCD($b));
2363 * echo $gcd->toString() . "\r\n"; // outputs 21
2364 * echo $a->toString() * $x->toString() + $b->toString() * $y->toString(); // outputs 21
2368 * @param Math_BigInteger $n
2369 * @return Math_BigInteger
2371 * @internal Calculates the GCD using the binary xGCD algorithim described in
2372 * {@link http://www.cacr.math.uwaterloo.ca/hac/about/chap14.pdf#page=19 HAC 14.61}. As the text above 14.61 notes,
2373 * the more traditional algorithim requires "relatively costly multiple-precision divisions".
2375 function extendedGCD($n)
2377 switch ( MATH_BIGINTEGER_MODE ) {
2378 case MATH_BIGINTEGER_MODE_GMP:
2379 extract(gmp_gcdext($this->value, $n->value));
2382 'gcd' => $this->_normalize(new Math_BigInteger($g)),
2383 'x' => $this->_normalize(new Math_BigInteger($s)),
2384 'y' => $this->_normalize(new Math_BigInteger($t))
2386 case MATH_BIGINTEGER_MODE_BCMATH:
2387 // it might be faster to use the binary xGCD algorithim here, as well, but (1) that algorithim works
2388 // best when the base is a power of 2 and (2) i don't think it'd make much difference, anyway. as is,
2389 // the basic extended euclidean algorithim is what we're using.
2399 while (bccomp($v, '0', 0) != 0) {
2400 $q = bcdiv($u, $v, 0);
2404 $v = bcsub($temp, bcmul($v, $q, 0), 0);
2408 $c = bcsub($temp, bcmul($a, $q, 0), 0);
2412 $d = bcsub($temp, bcmul($b, $q, 0), 0);
2416 'gcd' => $this->_normalize(new Math_BigInteger($u)),
2417 'x' => $this->_normalize(new Math_BigInteger($a)),
2418 'y' => $this->_normalize(new Math_BigInteger($b))
2424 $g = new Math_BigInteger();
2425 $g->value = array(1);
2427 while ( !(($x->value[0] & 1)|| ($y->value[0] & 1)) ) {
2436 $a = new Math_BigInteger();
2437 $b = new Math_BigInteger();
2438 $c = new Math_BigInteger();
2439 $d = new Math_BigInteger();
2441 $a->value = $d->value = $g->value = array(1);
2442 $b->value = $c->value = array();
2444 while ( !empty($u->value) ) {
2445 while ( !($u->value[0] & 1) ) {
2447 if ( (!empty($a->value) && ($a->value[0] & 1)) || (!empty($b->value) && ($b->value[0] & 1)) ) {
2449 $b = $b->subtract($x);
2455 while ( !($v->value[0] & 1) ) {
2457 if ( (!empty($d->value) && ($d->value[0] & 1)) || (!empty($c->value) && ($c->value[0] & 1)) ) {
2459 $d = $d->subtract($x);
2465 if ($u->compare($v) >= 0) {
2466 $u = $u->subtract($v);
2467 $a = $a->subtract($c);
2468 $b = $b->subtract($d);
2470 $v = $v->subtract($u);
2471 $c = $c->subtract($a);
2472 $d = $d->subtract($b);
2477 'gcd' => $this->_normalize($g->multiply($v)),
2478 'x' => $this->_normalize($c),
2479 'y' => $this->_normalize($d)
2484 * Calculates the greatest common divisor
2486 * Say you have 693 and 609. The GCD is 21.
2488 * Here's an example:
2491 * include('Math/BigInteger.php');
2493 * $a = new Math_BigInteger(693);
2494 * $b = new Math_BigInteger(609);
2496 * $gcd = a->extendedGCD($b);
2498 * echo $gcd->toString() . "\r\n"; // outputs 21
2502 * @param Math_BigInteger $n
2503 * @return Math_BigInteger
2508 extract($this->extendedGCD($n));
2515 * @return Math_BigInteger
2520 $temp = new Math_BigInteger();
2522 switch ( MATH_BIGINTEGER_MODE ) {
2523 case MATH_BIGINTEGER_MODE_GMP:
2524 $temp->value = gmp_abs($this->value);
2526 case MATH_BIGINTEGER_MODE_BCMATH:
2527 $temp->value = (bccomp($this->value, '0', 0) < 0) ? substr($this->value, 1) : $this->value;
2530 $temp->value = $this->value;
2537 * Compares two numbers.
2539 * Although one might think !$x->compare($y) means $x != $y, it, in fact, means the opposite. The reason for this is
2540 * demonstrated thusly:
2542 * $x > $y: $x->compare($y) > 0
2543 * $x < $y: $x->compare($y) < 0
2544 * $x == $y: $x->compare($y) == 0
2546 * Note how the same comparison operator is used. If you want to test for equality, use $x->equals($y).
2548 * @param Math_BigInteger $x
2549 * @return Integer < 0 if $this is less than $x; > 0 if $this is greater than $x, and 0 if they are equal.
2552 * @internal Could return $this->subtract($x), but that's not as fast as what we do do.
2554 function compare($y)
2556 switch ( MATH_BIGINTEGER_MODE ) {
2557 case MATH_BIGINTEGER_MODE_GMP:
2558 return gmp_cmp($this->value, $y->value);
2559 case MATH_BIGINTEGER_MODE_BCMATH:
2560 return bccomp($this->value, $y->value, 0);
2563 return $this->_compare($this->value, $this->is_negative, $y->value, $y->is_negative);
2567 * Compares two numbers.
2569 * @param Array $x_value
2570 * @param Boolean $x_negative
2571 * @param Array $y_value
2572 * @param Boolean $y_negative
2577 function _compare($x_value, $x_negative, $y_value, $y_negative)
2579 if ( $x_negative != $y_negative ) {
2580 return ( !$x_negative && $y_negative ) ? 1 : -1;
2583 $result = $x_negative ? -1 : 1;
2585 if ( count($x_value) != count($y_value) ) {
2586 return ( count($x_value) > count($y_value) ) ? $result : -$result;
2588 $size = max(count($x_value), count($y_value));
2590 $x_value = array_pad($x_value, $size, 0);
2591 $y_value = array_pad($y_value, $size, 0);
2593 for ($i = count($x_value) - 1; $i >= 0; --$i) {
2594 if ($x_value[$i] != $y_value[$i]) {
2595 return ( $x_value[$i] > $y_value[$i] ) ? $result : -$result;
2603 * Tests the equality of two numbers.
2605 * If you need to see if one number is greater than or less than another number, use Math_BigInteger::compare()
2607 * @param Math_BigInteger $x
2614 switch ( MATH_BIGINTEGER_MODE ) {
2615 case MATH_BIGINTEGER_MODE_GMP:
2616 return gmp_cmp($this->value, $x->value) == 0;
2618 return $this->value === $x->value && $this->is_negative == $x->is_negative;
2625 * Some bitwise operations give different results depending on the precision being used. Examples include left
2626 * shift, not, and rotates.
2628 * @param Math_BigInteger $x
2630 * @return Math_BigInteger
2632 function setPrecision($bits)
2634 $this->precision = $bits;
2635 if ( MATH_BIGINTEGER_MODE != MATH_BIGINTEGER_MODE_BCMATH ) {
2636 $this->bitmask = new Math_BigInteger(chr((1 << ($bits & 0x7)) - 1) . str_repeat(chr(0xFF), $bits >> 3), 256);
2638 $this->bitmask = new Math_BigInteger(bcpow('2', $bits, 0));
2641 $temp = $this->_normalize($this);
2642 $this->value = $temp->value;
2648 * @param Math_BigInteger $x
2650 * @internal Implemented per a request by Lluis Pamies i Juarez <lluis _a_ pamies.cat>
2651 * @return Math_BigInteger
2653 function bitwise_and($x)
2655 switch ( MATH_BIGINTEGER_MODE ) {
2656 case MATH_BIGINTEGER_MODE_GMP:
2657 $temp = new Math_BigInteger();
2658 $temp->value = gmp_and($this->value, $x->value);
2660 return $this->_normalize($temp);
2661 case MATH_BIGINTEGER_MODE_BCMATH:
2662 $left = $this->toBytes();
2663 $right = $x->toBytes();
2665 $length = max(strlen($left), strlen($right));
2667 $left = str_pad($left, $length, chr(0), STR_PAD_LEFT);
2668 $right = str_pad($right, $length, chr(0), STR_PAD_LEFT);
2670 return $this->_normalize(new Math_BigInteger($left & $right, 256));
2673 $result = $this->copy();
2675 $length = min(count($x->value), count($this->value));
2677 $result->value = array_slice($result->value, 0, $length);
2679 for ($i = 0; $i < $length; ++$i) {
2680 $result->value[$i] = $result->value[$i] & $x->value[$i];
2683 return $this->_normalize($result);
2689 * @param Math_BigInteger $x
2691 * @internal Implemented per a request by Lluis Pamies i Juarez <lluis _a_ pamies.cat>
2692 * @return Math_BigInteger
2694 function bitwise_or($x)
2696 switch ( MATH_BIGINTEGER_MODE ) {
2697 case MATH_BIGINTEGER_MODE_GMP:
2698 $temp = new Math_BigInteger();
2699 $temp->value = gmp_or($this->value, $x->value);
2701 return $this->_normalize($temp);
2702 case MATH_BIGINTEGER_MODE_BCMATH:
2703 $left = $this->toBytes();
2704 $right = $x->toBytes();
2706 $length = max(strlen($left), strlen($right));
2708 $left = str_pad($left, $length, chr(0), STR_PAD_LEFT);
2709 $right = str_pad($right, $length, chr(0), STR_PAD_LEFT);
2711 return $this->_normalize(new Math_BigInteger($left | $right, 256));
2714 $length = max(count($this->value), count($x->value));
2715 $result = $this->copy();
2716 $result->value = array_pad($result->value, 0, $length);
2717 $x->value = array_pad($x->value, 0, $length);
2719 for ($i = 0; $i < $length; ++$i) {
2720 $result->value[$i] = $this->value[$i] | $x->value[$i];
2723 return $this->_normalize($result);
2727 * Logical Exclusive-Or
2729 * @param Math_BigInteger $x
2731 * @internal Implemented per a request by Lluis Pamies i Juarez <lluis _a_ pamies.cat>
2732 * @return Math_BigInteger
2734 function bitwise_xor($x)
2736 switch ( MATH_BIGINTEGER_MODE ) {
2737 case MATH_BIGINTEGER_MODE_GMP:
2738 $temp = new Math_BigInteger();
2739 $temp->value = gmp_xor($this->value, $x->value);
2741 return $this->_normalize($temp);
2742 case MATH_BIGINTEGER_MODE_BCMATH:
2743 $left = $this->toBytes();
2744 $right = $x->toBytes();
2746 $length = max(strlen($left), strlen($right));
2748 $left = str_pad($left, $length, chr(0), STR_PAD_LEFT);
2749 $right = str_pad($right, $length, chr(0), STR_PAD_LEFT);
2751 return $this->_normalize(new Math_BigInteger($left ^ $right, 256));
2754 $length = max(count($this->value), count($x->value));
2755 $result = $this->copy();
2756 $result->value = array_pad($result->value, 0, $length);
2757 $x->value = array_pad($x->value, 0, $length);
2759 for ($i = 0; $i < $length; ++$i) {
2760 $result->value[$i] = $this->value[$i] ^ $x->value[$i];
2763 return $this->_normalize($result);
2770 * @internal Implemented per a request by Lluis Pamies i Juarez <lluis _a_ pamies.cat>
2771 * @return Math_BigInteger
2773 function bitwise_not()
2775 // calculuate "not" without regard to $this->precision
2776 // (will always result in a smaller number. ie. ~1 isn't 1111 1110 - it's 0)
2777 $temp = $this->toBytes();
2778 $pre_msb = decbin(ord($temp[0]));
2780 $msb = decbin(ord($temp[0]));
2781 if (strlen($msb) == 8) {
2782 $msb = substr($msb, strpos($msb, '0'));
2784 $temp[0] = chr(bindec($msb));
2786 // see if we need to add extra leading 1's
2787 $current_bits = strlen($pre_msb) + 8 * strlen($temp) - 8;
2788 $new_bits = $this->precision - $current_bits;
2789 if ($new_bits <= 0) {
2790 return $this->_normalize(new Math_BigInteger($temp, 256));
2793 // generate as many leading 1's as we need to.
2794 $leading_ones = chr((1 << ($new_bits & 0x7)) - 1) . str_repeat(chr(0xFF), $new_bits >> 3);
2795 $this->_base256_lshift($leading_ones, $current_bits);
2797 $temp = str_pad($temp, ceil($this->bits / 8), chr(0), STR_PAD_LEFT);
2799 return $this->_normalize(new Math_BigInteger($leading_ones | $temp, 256));
2803 * Logical Right Shift
2805 * Shifts BigInteger's by $shift bits, effectively dividing by 2**$shift.
2807 * @param Integer $shift
2808 * @return Math_BigInteger
2810 * @internal The only version that yields any speed increases is the internal version.
2812 function bitwise_rightShift($shift)
2814 $temp = new Math_BigInteger();
2816 switch ( MATH_BIGINTEGER_MODE ) {
2817 case MATH_BIGINTEGER_MODE_GMP:
2821 $two = gmp_init('2');
2824 $temp->value = gmp_div_q($this->value, gmp_pow($two, $shift));
2827 case MATH_BIGINTEGER_MODE_BCMATH:
2828 $temp->value = bcdiv($this->value, bcpow('2', $shift, 0), 0);
2831 default: // could just replace _lshift with this, but then all _lshift() calls would need to be rewritten
2832 // and I don't want to do that...
2833 $temp->value = $this->value;
2834 $temp->_rshift($shift);
2837 return $this->_normalize($temp);
2841 * Logical Left Shift
2843 * Shifts BigInteger's by $shift bits, effectively multiplying by 2**$shift.
2845 * @param Integer $shift
2846 * @return Math_BigInteger
2848 * @internal The only version that yields any speed increases is the internal version.
2850 function bitwise_leftShift($shift)
2852 $temp = new Math_BigInteger();
2854 switch ( MATH_BIGINTEGER_MODE ) {
2855 case MATH_BIGINTEGER_MODE_GMP:
2859 $two = gmp_init('2');
2862 $temp->value = gmp_mul($this->value, gmp_pow($two, $shift));
2865 case MATH_BIGINTEGER_MODE_BCMATH:
2866 $temp->value = bcmul($this->value, bcpow('2', $shift, 0), 0);
2869 default: // could just replace _rshift with this, but then all _lshift() calls would need to be rewritten
2870 // and I don't want to do that...
2871 $temp->value = $this->value;
2872 $temp->_lshift($shift);
2875 return $this->_normalize($temp);
2879 * Logical Left Rotate
2881 * Instead of the top x bits being dropped they're appended to the shifted bit string.
2883 * @param Integer $shift
2884 * @return Math_BigInteger
2887 function bitwise_leftRotate($shift)
2889 $bits = $this->toBytes();
2891 if ($this->precision > 0) {
2892 $precision = $this->precision;
2893 if ( MATH_BIGINTEGER_MODE == MATH_BIGINTEGER_MODE_BCMATH ) {
2894 $mask = $this->bitmask->subtract(new Math_BigInteger(1));
2895 $mask = $mask->toBytes();
2897 $mask = $this->bitmask->toBytes();
2900 $temp = ord($bits[0]);
2901 for ($i = 0; $temp >> $i; ++$i);
2902 $precision = 8 * strlen($bits) - 8 + $i;
2903 $mask = chr((1 << ($precision & 0x7)) - 1) . str_repeat(chr(0xFF), $precision >> 3);
2907 $shift+= $precision;
2909 $shift%= $precision;
2912 return $this->copy();
2915 $left = $this->bitwise_leftShift($shift);
2916 $left = $left->bitwise_and(new Math_BigInteger($mask, 256));
2917 $right = $this->bitwise_rightShift($precision - $shift);
2918 $result = MATH_BIGINTEGER_MODE != MATH_BIGINTEGER_MODE_BCMATH ? $left->bitwise_or($right) : $left->add($right);
2919 return $this->_normalize($result);
2923 * Logical Right Rotate
2925 * Instead of the bottom x bits being dropped they're prepended to the shifted bit string.
2927 * @param Integer $shift
2928 * @return Math_BigInteger
2931 function bitwise_rightRotate($shift)
2933 return $this->bitwise_leftRotate(-$shift);
2937 * Set random number generator function
2939 * $generator should be the name of a random generating function whose first parameter is the minimum
2940 * value and whose second parameter is the maximum value. If this function needs to be seeded, it should
2941 * be seeded prior to calling Math_BigInteger::random() or Math_BigInteger::randomPrime()
2943 * If the random generating function is not explicitly set, it'll be assumed to be mt_rand().
2946 * @see randomPrime()
2947 * @param optional String $generator
2950 function setRandomGenerator($generator)
2952 $this->generator = $generator;
2956 * Generate a random number
2958 * @param optional Integer $min
2959 * @param optional Integer $max
2960 * @return Math_BigInteger
2963 function random($min = false, $max = false)
2965 if ($min === false) {
2966 $min = new Math_BigInteger(0);
2969 if ($max === false) {
2970 $max = new Math_BigInteger(0x7FFFFFFF);
2973 $compare = $max->compare($min);
2976 return $this->_normalize($min);
2977 } else if ($compare < 0) {
2978 // if $min is bigger then $max, swap $min and $max
2984 $generator = $this->generator;
2986 $max = $max->subtract($min);
2987 $max = ltrim($max->toBytes(), chr(0));
2988 $size = strlen($max) - 1;
2992 for ($i = 0; $i < $bytes; ++$i) {
2993 $random.= chr($generator(0, 255));
2996 $blocks = $size >> 1;
2997 for ($i = 0; $i < $blocks; ++$i) {
2998 // mt_rand(-2147483648, 0x7FFFFFFF) always produces -2147483648 on some systems
2999 $random.= pack('n', $generator(0, 0xFFFF));
3002 $temp = new Math_BigInteger($random, 256);
3003 if ($temp->compare(new Math_BigInteger(substr($max, 1), 256)) > 0) {
3004 $random = chr($generator(0, ord($max[0]) - 1)) . $random;
3006 $random = chr($generator(0, ord($max[0]) )) . $random;
3009 $random = new Math_BigInteger($random, 256);
3011 return $this->_normalize($random->add($min));
3015 * Generate a random prime number.
3017 * If there's not a prime within the given range, false will be returned. If more than $timeout seconds have elapsed,
3018 * give up and return false.
3020 * @param optional Integer $min
3021 * @param optional Integer $max
3022 * @param optional Integer $timeout
3023 * @return Math_BigInteger
3025 * @internal See {@link http://www.cacr.math.uwaterloo.ca/hac/about/chap4.pdf#page=15 HAC 4.44}.
3027 function randomPrime($min = false, $max = false, $timeout = false)
3029 $compare = $max->compare($min);
3033 } else if ($compare < 0) {
3034 // if $min is bigger then $max, swap $min and $max
3040 // gmp_nextprime() requires PHP 5 >= 5.2.0 per <http://php.net/gmp-nextprime>.
3041 if ( MATH_BIGINTEGER_MODE == MATH_BIGINTEGER_MODE_GMP && function_exists('gmp_nextprime') ) {
3042 // we don't rely on Math_BigInteger::random()'s min / max when gmp_nextprime() is being used since this function
3043 // does its own checks on $max / $min when gmp_nextprime() is used. When gmp_nextprime() is not used, however,
3044 // the same $max / $min checks are not performed.
3045 if ($min === false) {
3046 $min = new Math_BigInteger(0);
3049 if ($max === false) {
3050 $max = new Math_BigInteger(0x7FFFFFFF);
3053 $x = $this->random($min, $max);
3055 $x->value = gmp_nextprime($x->value);
3057 if ($x->compare($max) <= 0) {
3061 $x->value = gmp_nextprime($min->value);
3063 if ($x->compare($max) <= 0) {
3072 $one = new Math_BigInteger(1);
3073 $two = new Math_BigInteger(2);
3078 $x = $this->random($min, $max);
3079 if ($x->equals($two)) {
3084 if ($x->compare($max) > 0) {
3085 // if $x > $max then $max is even and if $min == $max then no prime number exists between the specified range
3086 if ($min->equals($max)) {
3093 $initial_x = $x->copy();
3096 if ($timeout !== false && time() - $start > $timeout) {
3100 if ($x->isPrime()) {
3106 if ($x->compare($max) > 0) {
3108 if ($x->equals($two)) {
3114 if ($x->equals($initial_x)) {
3121 * Make the current number odd
3123 * If the current number is odd it'll be unchanged. If it's even, one will be added to it.
3125 * @see randomPrime()
3128 function _make_odd()
3130 switch ( MATH_BIGINTEGER_MODE ) {
3131 case MATH_BIGINTEGER_MODE_GMP:
3132 gmp_setbit($this->value, 0);
3134 case MATH_BIGINTEGER_MODE_BCMATH:
3135 if ($this->value[strlen($this->value) - 1] % 2 == 0) {
3136 $this->value = bcadd($this->value, '1');
3140 $this->value[0] |= 1;
3145 * Checks a numer to see if it's prime
3147 * Assuming the $t parameter is not set, this function has an error rate of 2**-80. The main motivation for the
3148 * $t parameter is distributability. Math_BigInteger::randomPrime() can be distributed accross multiple pageloads
3149 * on a website instead of just one.
3151 * @param optional Integer $t
3154 * @internal Uses the
3155 * {@link http://en.wikipedia.org/wiki/Miller%E2%80%93Rabin_primality_test Miller-Rabin primality test}. See
3156 * {@link http://www.cacr.math.uwaterloo.ca/hac/about/chap4.pdf#page=8 HAC 4.24}.
3158 function isPrime($t = false)
3160 $length = strlen($this->toBytes());
3163 // see HAC 4.49 "Note (controlling the error probability)"
3164 if ($length >= 163) { $t = 2; } // floor(1300 / 8)
3165 else if ($length >= 106) { $t = 3; } // floor( 850 / 8)
3166 else if ($length >= 81 ) { $t = 4; } // floor( 650 / 8)
3167 else if ($length >= 68 ) { $t = 5; } // floor( 550 / 8)
3168 else if ($length >= 56 ) { $t = 6; } // floor( 450 / 8)
3169 else if ($length >= 50 ) { $t = 7; } // floor( 400 / 8)
3170 else if ($length >= 43 ) { $t = 8; } // floor( 350 / 8)
3171 else if ($length >= 37 ) { $t = 9; } // floor( 300 / 8)
3172 else if ($length >= 31 ) { $t = 12; } // floor( 250 / 8)
3173 else if ($length >= 25 ) { $t = 15; } // floor( 200 / 8)
3174 else if ($length >= 18 ) { $t = 18; } // floor( 150 / 8)
3178 // ie. gmp_testbit($this, 0)
3179 // ie. isEven() or !isOdd()
3180 switch ( MATH_BIGINTEGER_MODE ) {
3181 case MATH_BIGINTEGER_MODE_GMP:
3182 return gmp_prob_prime($this->value, $t) != 0;
3183 case MATH_BIGINTEGER_MODE_BCMATH:
3184 if ($this->value === '2') {
3187 if ($this->value[strlen($this->value) - 1] % 2 == 0) {
3192 if ($this->value == array(2)) {
3195 if (~$this->value[0] & 1) {
3200 static $primes, $zero, $one, $two;
3202 if (!isset($primes)) {
3204 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59,
3205 61, 67, 71, 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 127, 131, 137,
3206 139, 149, 151, 157, 163, 167, 173, 179, 181, 191, 193, 197, 199, 211, 223, 227,
3207 229, 233, 239, 241, 251, 257, 263, 269, 271, 277, 281, 283, 293, 307, 311, 313,
3208 317, 331, 337, 347, 349, 353, 359, 367, 373, 379, 383, 389, 397, 401, 409, 419,
3209 421, 431, 433, 439, 443, 449, 457, 461, 463, 467, 479, 487, 491, 499, 503, 509,
3210 521, 523, 541, 547, 557, 563, 569, 571, 577, 587, 593, 599, 601, 607, 613, 617,
3211 619, 631, 641, 643, 647, 653, 659, 661, 673, 677, 683, 691, 701, 709, 719, 727,
3212 733, 739, 743, 751, 757, 761, 769, 773, 787, 797, 809, 811, 821, 823, 827, 829,
3213 839, 853, 857, 859, 863, 877, 881, 883, 887, 907, 911, 919, 929, 937, 941, 947,
3214 953, 967, 971, 977, 983, 991, 997
3217 if ( MATH_BIGINTEGER_MODE != MATH_BIGINTEGER_MODE_INTERNAL ) {
3218 for ($i = 0; $i < count($primes); ++$i) {
3219 $primes[$i] = new Math_BigInteger($primes[$i]);
3223 $zero = new Math_BigInteger();
3224 $one = new Math_BigInteger(1);
3225 $two = new Math_BigInteger(2);
3228 if ($this->equals($one)) {
3232 // see HAC 4.4.1 "Random search for probable primes"
3233 if ( MATH_BIGINTEGER_MODE != MATH_BIGINTEGER_MODE_INTERNAL ) {
3234 foreach ($primes as $prime) {
3235 list(, $r) = $this->divide($prime);
3236 if ($r->equals($zero)) {
3237 return $this->equals($prime);
3241 $value = $this->value;
3242 foreach ($primes as $prime) {
3243 list(, $r) = $this->_divide_digit($value, $prime);
3245 return count($value) == 1 && $value[0] == $prime;
3251 $n_1 = $n->subtract($one);
3252 $n_2 = $n->subtract($two);
3255 $r_value = $r->value;
3256 // ie. $s = gmp_scan1($n, 0) and $r = gmp_div_q($n, gmp_pow(gmp_init('2'), $s));
3257 if ( MATH_BIGINTEGER_MODE == MATH_BIGINTEGER_MODE_BCMATH ) {
3259 // if $n was 1, $r would be 0 and this would be an infinite loop, hence our $this->equals($one) check earlier
3260 while ($r->value[strlen($r->value) - 1] % 2 == 0) {
3261 $r->value = bcdiv($r->value, '2', 0);
3265 for ($i = 0, $r_length = count($r_value); $i < $r_length; ++$i) {
3266 $temp = ~$r_value[$i] & 0xFFFFFF;
3267 for ($j = 1; ($temp >> $j) & 1; ++$j);
3272 $s = 26 * $i + $j - 1;
3276 for ($i = 0; $i < $t; ++$i) {
3277 $a = $this->random($two, $n_2);
3278 $y = $a->modPow($r, $n);
3280 if (!$y->equals($one) && !$y->equals($n_1)) {
3281 for ($j = 1; $j < $s && !$y->equals($n_1); ++$j) {
3282 $y = $y->modPow($two, $n);
3283 if ($y->equals($one)) {
3288 if (!$y->equals($n_1)) {
3297 * Logical Left Shift
3299 * Shifts BigInteger's by $shift bits.
3301 * @param Integer $shift
3304 function _lshift($shift)
3306 if ( $shift == 0 ) {
3310 $num_digits = (int) ($shift / 26);
3312 $shift = 1 << $shift;
3316 for ($i = 0; $i < count($this->value); ++$i) {
3317 $temp = $this->value[$i] * $shift + $carry;
3318 $carry = (int) ($temp / 0x4000000);
3319 $this->value[$i] = (int) ($temp - $carry * 0x4000000);
3323 $this->value[] = $carry;
3326 while ($num_digits--) {
3327 array_unshift($this->value, 0);
3332 * Logical Right Shift
3334 * Shifts BigInteger's by $shift bits.
3336 * @param Integer $shift
3339 function _rshift($shift)
3345 $num_digits = (int) ($shift / 26);
3347 $carry_shift = 26 - $shift;
3348 $carry_mask = (1 << $shift) - 1;
3350 if ( $num_digits ) {
3351 $this->value = array_slice($this->value, $num_digits);
3356 for ($i = count($this->value) - 1; $i >= 0; --$i) {
3357 $temp = $this->value[$i] >> $shift | $carry;
3358 $carry = ($this->value[$i] & $carry_mask) << $carry_shift;
3359 $this->value[$i] = $temp;
3362 $this->value = $this->_trim($this->value);
3368 * Removes leading zeros and truncates (if necessary) to maintain the appropriate precision
3370 * @param Math_BigInteger
3371 * @return Math_BigInteger
3375 function _normalize($result)
3377 $result->precision = $this->precision;
3378 $result->bitmask = $this->bitmask;
3380 switch ( MATH_BIGINTEGER_MODE ) {
3381 case MATH_BIGINTEGER_MODE_GMP:
3382 if (!empty($result->bitmask->value)) {
3383 $result->value = gmp_and($result->value, $result->bitmask->value);
3387 case MATH_BIGINTEGER_MODE_BCMATH:
3388 if (!empty($result->bitmask->value)) {
3389 $result->value = bcmod($result->value, $result->bitmask->value);
3395 $value = &$result->value;
3397 if ( !count($value) ) {
3401 $value = $this->_trim($value);
3403 if (!empty($result->bitmask->value)) {
3404 $length = min(count($value), count($this->bitmask->value));
3405 $value = array_slice($value, 0, $length);
3407 for ($i = 0; $i < $length; ++$i) {
3408 $value[$i] = $value[$i] & $this->bitmask->value[$i];
3418 * Removes leading zeros
3420 * @return Math_BigInteger
3423 function _trim($value)
3425 for ($i = count($value) - 1; $i >= 0; --$i) {
3438 * @param $input Array
3439 * @param $multiplier mixed
3443 function _array_repeat($input, $multiplier)
3445 return ($multiplier) ? array_fill(0, $multiplier, $input) : array();
3449 * Logical Left Shift
3451 * Shifts binary strings $shift bits, essentially multiplying by 2**$shift.
3454 * @param $shift Integer
3458 function _base256_lshift(&$x, $shift)
3464 $num_bytes = $shift >> 3; // eg. floor($shift/8)
3465 $shift &= 7; // eg. $shift % 8
3468 for ($i = strlen($x) - 1; $i >= 0; --$i) {
3469 $temp = ord($x[$i]) << $shift | $carry;
3470 $x[$i] = chr($temp);
3471 $carry = $temp >> 8;
3473 $carry = ($carry != 0) ? chr($carry) : '';
3474 $x = $carry . $x . str_repeat(chr(0), $num_bytes);
3478 * Logical Right Shift
3480 * Shifts binary strings $shift bits, essentially dividing by 2**$shift and returning the remainder.
3483 * @param $shift Integer
3487 function _base256_rshift(&$x, $shift)
3490 $x = ltrim($x, chr(0));
3494 $num_bytes = $shift >> 3; // eg. floor($shift/8)
3495 $shift &= 7; // eg. $shift % 8
3499 $start = $num_bytes > strlen($x) ? -strlen($x) : -$num_bytes;
3500 $remainder = substr($x, $start);
3501 $x = substr($x, 0, -$num_bytes);
3505 $carry_shift = 8 - $shift;
3506 for ($i = 0; $i < strlen($x); ++$i) {
3507 $temp = (ord($x[$i]) >> $shift) | $carry;
3508 $carry = (ord($x[$i]) << $carry_shift) & 0xFF;
3509 $x[$i] = chr($temp);
3511 $x = ltrim($x, chr(0));
3513 $remainder = chr($carry >> $carry_shift) . $remainder;
3515 return ltrim($remainder, chr(0));
3518 // one quirk about how the following functions are implemented is that PHP defines N to be an unsigned long
3519 // at 32-bits, while java's longs are 64-bits.
3522 * Converts 32-bit integers to bytes.
3528 function _int2bytes($x)
3530 return ltrim(pack('N', $x), chr(0));
3534 * Converts bytes to 32-bit integers
3540 function _bytes2int($x)
3542 $temp = unpack('Nint', str_pad($x, 4, chr(0), STR_PAD_LEFT));
3543 return $temp['int'];