plugins/OStatus/extlib/Math/BigInteger.php

   1 <?php
   2 /* vim: set expandtab tabstop=4 shiftwidth=4 softtabstop=4: */
   3
   4 /**
   5  * Pure-PHP arbitrary precision integer arithmetic library.
   6  *
   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.
   9  *
  10  * PHP versions 4 and 5
  11  *
  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)
  14  *
  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.
  22  *
  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 /
  25  * subtraction).
  26  *
  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)
  29  *
  30  * Useful resources are as follows:
  31  *
  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
  35  *
  36  * Here's an example of how to use this library:
  37  * <code>
  38  * <?php
  39  *    include('Math/BigInteger.php');
  40  *
  41  *    $a = new Math_BigInteger(2);
  42  *    $b = new Math_BigInteger(3);
  43  *
  44  *    $c = $a->add($b);
  45  *
  46  *    echo $c->toString(); // outputs 5
  47  * ?>
  48  * </code>
  49  *
  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.
  54  *
  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.
  59  *
  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,
  63  * MA  02111-1307  USA
  64  *
  65  * @category   Math
  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
  72  */
  73
  74 /**#@+
  75  * Reduction constants
  76  *
  77  * @access private
  78  * @see Math_BigInteger::_reduce()
  79  */
  80 /**
  81  * @see Math_BigInteger::_montgomery()
  82  * @see Math_BigInteger::_prepMontgomery()
  83  */
  84 define('MATH_BIGINTEGER_MONTGOMERY', 0);
  85 /**
  86  * @see Math_BigInteger::_barrett()
  87  */
  88 define('MATH_BIGINTEGER_BARRETT', 1);
  89 /**
  90  * @see Math_BigInteger::_mod2()
  91  */
  92 define('MATH_BIGINTEGER_POWEROF2', 2);
  93 /**
  94  * @see Math_BigInteger::_remainder()
  95  */
  96 define('MATH_BIGINTEGER_CLASSIC', 3);
  97 /**
  98  * @see Math_BigInteger::__clone()
  99  */
 100 define('MATH_BIGINTEGER_NONE', 4);
 101 /**#@-*/
 102
 103 /**#@+
 104  * Array constants
 105  *
 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.
 108  *
 109  * @access private
 110  */
 111 /**
 112  * $result[MATH_BIGINTEGER_VALUE] contains the value.
 113  */
 114 define('MATH_BIGINTEGER_VALUE', 0);
 115 /**
 116  * $result[MATH_BIGINTEGER_SIGN] contains the sign.
 117  */
 118 define('MATH_BIGINTEGER_SIGN', 1);
 119 /**#@-*/
 120
 121 /**#@+
 122  * @access private
 123  * @see Math_BigInteger::_montgomery()
 124  * @see Math_BigInteger::_barrett()
 125  */
 126 /**
 127  * Cache constants
 128  *
 129  * $cache[MATH_BIGINTEGER_VARIABLE] tells us whether or not the cached data is still valid.
 130  */
 131 define('MATH_BIGINTEGER_VARIABLE', 0);
 132 /**
 133  * $cache[MATH_BIGINTEGER_DATA] contains the cached data.
 134  */
 135 define('MATH_BIGINTEGER_DATA', 1);
 136 /**#@-*/
 137
 138 /**#@+
 139  * Mode constants.
 140  *
 141  * @access private
 142  * @see Math_BigInteger::Math_BigInteger()
 143  */
 144 /**
 145  * To use the pure-PHP implementation
 146  */
 147 define('MATH_BIGINTEGER_MODE_INTERNAL', 1);
 148 /**
 149  * To use the BCMath library
 150  *
 151  * (if enabled; otherwise, the internal implementation will be used)
 152  */
 153 define('MATH_BIGINTEGER_MODE_BCMATH', 2);
 154 /**
 155  * To use the GMP library
 156  *
 157  * (if present; otherwise, either the BCMath or the internal implementation will be used)
 158  */
 159 define('MATH_BIGINTEGER_MODE_GMP', 3);
 160 /**#@-*/
 161
 162 /**
 163  * The largest digit that may be used in addition / subtraction
 164  *
 165  * (we do pow(2, 52) instead of using 4503599627370496, directly, because some PHP installations
 166  *  will truncate 4503599627370496)
 167  *
 168  * @access private
 169  */
 170 define('MATH_BIGINTEGER_MAX_DIGIT52', pow(2, 52));
 171
 172 /**
 173  * Karatsuba Cutoff
 174  *
 175  * At what point do we switch between Karatsuba multiplication and schoolbook long multiplication?
 176  *
 177  * @access private
 178  */
 179 define('MATH_BIGINTEGER_KARATSUBA_CUTOFF', 25);
 180
 181 /**
 182  * Pure-PHP arbitrary precision integer arithmetic library. Supports base-2, base-10, base-16, and base-256
 183  * numbers.
 184  *
 185  * @author  Jim Wigginton <terrafrost@php.net>
 186  * @version 1.0.0RC4
 187  * @access  public
 188  * @package Math_BigInteger
 189  */
 190 class Math_BigInteger {
 191     /**
 192      * Holds the BigInteger's value.
 193      *
 194      * @var Array
 195      * @access private
 196      */
 197     var $value;
 198
 199     /**
 200      * Holds the BigInteger's magnitude.
 201      *
 202      * @var Boolean
 203      * @access private
 204      */
 205     var $is_negative = false;
 206
 207     /**
 208      * Random number generator function
 209      *
 210      * @see setRandomGenerator()
 211      * @access private
 212      */
 213     var $generator = 'mt_rand';
 214
 215     /**
 216      * Precision
 217      *
 218      * @see setPrecision()
 219      * @access private
 220      */
 221     var $precision = -1;
 222
 223     /**
 224      * Precision Bitmask
 225      *
 226      * @see setPrecision()
 227      * @access private
 228      */
 229     var $bitmask = false;
 230
 231     /**
 232      * Mode independant value used for serialization.
 233      *
 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.
 237      *
 238      * @see __sleep()
 239      * @see __wakeup()
 240      * @var String
 241      * @access private
 242      */
 243     var $hex;
 244
 245     /**
 246      * Converts base-2, base-10, base-16, and binary strings (eg. base-256) to BigIntegers.
 247      *
 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.
 250      *
 251      * Here's an example:
 252      * <code>
 253      * <?php
 254      *    include('Math/BigInteger.php');
 255      *
 256      *    $a = new Math_BigInteger('0x32', 16); // 50 in base-16
 257      *
 258      *    echo $a->toString(); // outputs 50
 259      * ?>
 260      * </code>
 261      *
 262      * @param optional $x base-10 number or base-$base number if $base set.
 263      * @param optional integer $base
 264      * @return Math_BigInteger
 265      * @access public
 266      */
 267     function Math_BigInteger($x = 0, $base = 10)
 268     {
 269         if ( !defined('MATH_BIGINTEGER_MODE') ) {
 270             switch (true) {
 271                 case extension_loaded('gmp'):
 272                     define('MATH_BIGINTEGER_MODE', MATH_BIGINTEGER_MODE_GMP);
 273                     break;
 274                 case extension_loaded('bcmath'):
 275                     define('MATH_BIGINTEGER_MODE', MATH_BIGINTEGER_MODE_BCMATH);
 276                     break;
 277                 default:
 278                     define('MATH_BIGINTEGER_MODE', MATH_BIGINTEGER_MODE_INTERNAL);
 279             }
 280         }
 281
 282         switch ( MATH_BIGINTEGER_MODE ) {
 283             case MATH_BIGINTEGER_MODE_GMP:
 284                 if (is_resource($x) && get_resource_type($x) == 'GMP integer') {
 285                     $this->value = $x;
 286                     return;
 287                 }
 288                 $this->value = gmp_init(0);
 289                 break;
 290             case MATH_BIGINTEGER_MODE_BCMATH:
 291                 $this->value = '0';
 292                 break;
 293             default:
 294                 $this->value = array();
 295         }
 296
 297         if (empty($x)) {
 298             return;
 299         }
 300
 301         switch ($base) {
 302             case -256:
 303                 if (ord($x[0]) & 0x80) {
 304                     $x = ~$x;
 305                     $this->is_negative = true;
 306                 }
 307             case  256:
 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));
 312                         break;
 313                     case MATH_BIGINTEGER_MODE_BCMATH:
 314                         // round $len to the nearest 4 (thanks, DavidMJ!)
 315                         $len = (strlen($x) + 3) & 0xFFFFFFFC;
 316
 317                         $x = str_pad($x, $len, chr(0), STR_PAD_LEFT);
 318
 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);
 322                         }
 323
 324                         if ($this->is_negative) {
 325                             $this->value = '-' . $this->value;
 326                         }
 327
 328                         break;
 329                     // converts a base-2**8 (big endian / msb) number to base-2**26 (little endian / lsb)
 330                     default:
 331                         while (strlen($x)) {
 332                             $this->value[] = $this->_bytes2int($this->_base256_rshift($x, 26));
 333                         }
 334                 }
 335
 336                 if ($this->is_negative) {
 337                     if (MATH_BIGINTEGER_MODE != MATH_BIGINTEGER_MODE_INTERNAL) {
 338                         $this->is_negative = false;
 339                     }
 340                     $temp = $this->add(new Math_BigInteger('-1'));
 341                     $this->value = $temp->value;
 342                 }
 343                 break;
 344             case  16:
 345             case -16:
 346                 if ($base > 0 && $x[0] == '-') {
 347                     $this->is_negative = true;
 348                     $x = substr($x, 1);
 349                 }
 350
 351                 $x = preg_replace('#^(?:0x)?([A-Fa-f0-9]*).*#', '$1', $x);
 352
 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));
 357                 }
 358
 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;
 364                         break;
 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;
 370                         break;
 371                     default:
 372                         $x = ( strlen($x) & 1 ) ? '0' . $x : $x;
 373                         $temp = new Math_BigInteger(pack('H*', $x), 256);
 374                         $this->value = $temp->value;
 375                 }
 376
 377                 if ($is_negative) {
 378                     $temp = $this->add(new Math_BigInteger('-1'));
 379                     $this->value = $temp->value;
 380                 }
 381                 break;
 382             case  10:
 383             case -10:
 384                 $x = preg_replace('#^(-?[0-9]*).*#', '$1', $x);
 385
 386                 switch ( MATH_BIGINTEGER_MODE ) {
 387                     case MATH_BIGINTEGER_MODE_GMP:
 388                         $this->value = gmp_init($x);
 389                         break;
 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;
 394                         break;
 395                     default:
 396                         $temp = new Math_BigInteger();
 397
 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);
 401
 402                         if ($x[0] == '-') {
 403                             $this->is_negative = true;
 404                             $x = substr($x, 1);
 405                         }
 406
 407                         $x = str_pad($x, strlen($x) + (6 * strlen($x)) % 7, 0, STR_PAD_LEFT);
 408
 409                         while (strlen($x)) {
 410                             $temp = $temp->multiply($multiplier);
 411                             $temp = $temp->add(new Math_BigInteger($this->_int2bytes(substr($x, 0, 7)), 256));
 412                             $x = substr($x, 7);
 413                         }
 414
 415                         $this->value = $temp->value;
 416                 }
 417                 break;
 418             case  2: // base-2 support originally implemented by Lluis Pamies - thanks!
 419             case -2:
 420                 if ($base > 0 && $x[0] == '-') {
 421                     $this->is_negative = true;
 422                     $x = substr($x, 1);
 423                 }
 424
 425                 $x = preg_replace('#^([01]*).*#', '$1', $x);
 426                 $x = str_pad($x, strlen($x) + (3 * strlen($x)) % 4, 0, STR_PAD_LEFT);
 427
 428                 $str = '0x';
 429                 while (strlen($x)) {
 430                     $part = substr($x, 0, 4);
 431                     $str.= dechex(bindec($part));
 432                     $x = substr($x, 4);
 433                 }
 434
 435                 if ($this->is_negative) {
 436                     $str = '-' . $str;
 437                 }
 438
 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;
 442
 443                 break;
 444             default:
 445                 // base not supported, so we'll let $this == 0
 446         }
 447     }
 448
 449     /**
 450      * Converts a BigInteger to a byte string (eg. base-256).
 451      *
 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.
 454      *
 455      * Here's an example:
 456      * <code>
 457      * <?php
 458      *    include('Math/BigInteger.php');
 459      *
 460      *    $a = new Math_BigInteger('65');
 461      *
 462      *    echo $a->toBytes(); // outputs chr(65)
 463      * ?>
 464      * </code>
 465      *
 466      * @param Boolean $twos_compliment
 467      * @return String
 468      * @access public
 469      * @internal Converts a base-2**26 number to base-2**8
 470      */
 471     function toBytes($twos_compliment = false)
 472     {
 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) : '';
 477             }
 478
 479             $temp = $comparison < 0 ? $this->add(new Math_BigInteger(1)) : $this->copy();
 480             $bytes = $temp->toBytes();
 481
 482             if (empty($bytes)) { // eg. if the number we're trying to convert is -1
 483                 $bytes = chr(0);
 484             }
 485
 486             if (ord($bytes[0]) & 0x80) {
 487                 $bytes = chr(0) . $bytes;
 488             }
 489
 490             return $comparison < 0 ? ~$bytes : $bytes;
 491         }
 492
 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) : '';
 497                 }
 498
 499                 $temp = gmp_strval(gmp_abs($this->value), 16);
 500                 $temp = ( strlen($temp) & 1 ) ? '0' . $temp : $temp;
 501                 $temp = pack('H*', $temp);
 502
 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) : '';
 509                 }
 510
 511                 $value = '';
 512                 $current = $this->value;
 513
 514                 if ($current[0] == '-') {
 515                     $current = substr($current, 1);
 516                 }
 517
 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);
 522                 }
 523
 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));
 527         }
 528
 529         if (!count($this->value)) {
 530             return $this->precision > 0 ? str_repeat(chr(0), ($this->precision + 1) >> 3) : '';
 531         }
 532         $result = $this->_int2bytes($this->value[count($this->value) - 1]);
 533
 534         $temp = $this->copy();
 535
 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);
 539         }
 540
 541         return $this->precision > 0 ?
 542             str_pad(substr($result, -(($this->precision + 7) >> 3)), ($this->precision + 7) >> 3, chr(0), STR_PAD_LEFT) :
 543             $result;
 544     }
 545
 546     /**
 547      * Converts a BigInteger to a hex string (eg. base-16)).
 548      *
 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.
 551      *
 552      * Here's an example:
 553      * <code>
 554      * <?php
 555      *    include('Math/BigInteger.php');
 556      *
 557      *    $a = new Math_BigInteger('65');
 558      *
 559      *    echo $a->toHex(); // outputs '41'
 560      * ?>
 561      * </code>
 562      *
 563      * @param Boolean $twos_compliment
 564      * @return String
 565      * @access public
 566      * @internal Converts a base-2**26 number to base-2**8
 567      */
 568     function toHex($twos_compliment = false)
 569     {
 570         return bin2hex($this->toBytes($twos_compliment));
 571     }
 572
 573     /**
 574      * Converts a BigInteger to a bit string (eg. base-2).
 575      *
 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.
 578      *
 579      * Here's an example:
 580      * <code>
 581      * <?php
 582      *    include('Math/BigInteger.php');
 583      *
 584      *    $a = new Math_BigInteger('65');
 585      *
 586      *    echo $a->toBits(); // outputs '1000001'
 587      * ?>
 588      * </code>
 589      *
 590      * @param Boolean $twos_compliment
 591      * @return String
 592      * @access public
 593      * @internal Converts a base-2**26 number to base-2**2
 594      */
 595     function toBits($twos_compliment = false)
 596     {
 597         $hex = $this->toHex($twos_compliment);
 598         $bits = '';
 599         for ($i = 0; $i < strlen($hex); $i+=8) {
 600             $bits.= str_pad(decbin(hexdec(substr($hex, $i, 8))), 32, '0', STR_PAD_LEFT);
 601         }
 602         return $this->precision > 0 ? substr($bits, -$this->precision) : ltrim($bits, '0');
 603     }
 604
 605     /**
 606      * Converts a BigInteger to a base-10 number.
 607      *
 608      * Here's an example:
 609      * <code>
 610      * <?php
 611      *    include('Math/BigInteger.php');
 612      *
 613      *    $a = new Math_BigInteger('50');
 614      *
 615      *    echo $a->toString(); // outputs 50
 616      * ?>
 617      * </code>
 618      *
 619      * @return String
 620      * @access public
 621      * @internal Converts a base-2**26 number to base-10**7 (which is pretty much base-10)
 622      */
 623     function toString()
 624     {
 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') {
 630                     return '0';
 631                 }
 632
 633                 return ltrim($this->value, '0');
 634         }
 635
 636         if (!count($this->value)) {
 637             return '0';
 638         }
 639
 640         $temp = $this->copy();
 641         $temp->is_negative = false;
 642
 643         $divisor = new Math_BigInteger();
 644         $divisor->value = array(10000000); // eg. 10**7
 645         $result = '';
 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;
 649         }
 650         $result = ltrim($result, '0');
 651         if (empty($result)) {
 652             $result = '0';
 653         }
 654
 655         if ($this->is_negative) {
 656             $result = '-' . $result;
 657         }
 658
 659         return $result;
 660     }
 661
 662     /**
 663      * Copy an object
 664      *
 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:
 667      *
 668      * {@link http://php.net/language.oop5.basic#51624}
 669      *
 670      * @access public
 671      * @see __clone()
 672      * @return Math_BigInteger
 673      */
 674     function copy()
 675     {
 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;
 682         return $temp;
 683     }
 684
 685     /**
 686      *  __toString() magic method
 687      *
 688      * Will be called, automatically, if you're supporting just PHP5.  If you're supporting PHP4, you'll need to call
 689      * toString().
 690      *
 691      * @access public
 692      * @internal Implemented per a suggestion by Techie-Michael - thanks!
 693      */
 694     function __toString()
 695     {
 696         return $this->toString();
 697     }
 698
 699     /**
 700      * __clone() magic method
 701      *
 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.
 706      *
 707      * @access public
 708      * @see copy()
 709      * @return Math_BigInteger
 710      */
 711     function __clone()
 712     {
 713         return $this->copy();
 714     }
 715
 716     /**
 717      *  __sleep() magic method
 718      *
 719      * Will be called, automatically, when serialize() is called on a Math_BigInteger object.
 720      *
 721      * @see __wakeup()
 722      * @access public
 723      */
 724     function __sleep()
 725     {
 726         $this->hex = $this->toHex(true);
 727         $vars = array('hex');
 728         if ($this->generator != 'mt_rand') {
 729             $vars[] = 'generator';
 730         }
 731         if ($this->precision > 0) {
 732             $vars[] = 'precision';
 733         }
 734         return $vars;
 735
 736     }
 737
 738     /**
 739      *  __wakeup() magic method
 740      *
 741      * Will be called, automatically, when unserialize() is called on a Math_BigInteger object.
 742      *
 743      * @see __sleep()
 744      * @access public
 745      */
 746     function __wakeup()
 747     {
 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);
 755         }
 756     }
 757
 758     /**
 759      * Adds two BigIntegers.
 760      *
 761      * Here's an example:
 762      * <code>
 763      * <?php
 764      *    include('Math/BigInteger.php');
 765      *
 766      *    $a = new Math_BigInteger('10');
 767      *    $b = new Math_BigInteger('20');
 768      *
 769      *    $c = $a->add($b);
 770      *
 771      *    echo $c->toString(); // outputs 30
 772      * ?>
 773      * </code>
 774      *
 775      * @param Math_BigInteger $y
 776      * @return Math_BigInteger
 777      * @access public
 778      * @internal Performs base-2**52 addition
 779      */
 780     function add($y)
 781     {
 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);
 786
 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);
 791
 792                 return $this->_normalize($temp);
 793         }
 794
 795         $temp = $this->_add($this->value, $this->is_negative, $y->value, $y->is_negative);
 796
 797         $result = new Math_BigInteger();
 798         $result->value = $temp[MATH_BIGINTEGER_VALUE];
 799         $result->is_negative = $temp[MATH_BIGINTEGER_SIGN];
 800
 801         return $this->_normalize($result);
 802     }
 803
 804     /**
 805      * Performs addition.
 806      *
 807      * @param Array $x_value
 808      * @param Boolean $x_negative
 809      * @param Array $y_value
 810      * @param Boolean $y_negative
 811      * @return Array
 812      * @access private
 813      */
 814     function _add($x_value, $x_negative, $y_value, $y_negative)
 815     {
 816         $x_size = count($x_value);
 817         $y_size = count($y_value);
 818
 819         if ($x_size == 0) {
 820             return array(
 821                 MATH_BIGINTEGER_VALUE => $y_value,
 822                 MATH_BIGINTEGER_SIGN => $y_negative
 823             );
 824         } else if ($y_size == 0) {
 825             return array(
 826                 MATH_BIGINTEGER_VALUE => $x_value,
 827                 MATH_BIGINTEGER_SIGN => $x_negative
 828             );
 829         }
 830
 831         // subtract, if appropriate
 832         if ( $x_negative != $y_negative ) {
 833             if ( $x_value == $y_value ) {
 834                 return array(
 835                     MATH_BIGINTEGER_VALUE => array(),
 836                     MATH_BIGINTEGER_SIGN => false
 837                 );
 838             }
 839
 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;
 843
 844             return $temp;
 845         }
 846
 847         if ($x_size < $y_size) {
 848             $size = $x_size;
 849             $value = $y_value;
 850         } else {
 851             $size = $y_size;
 852             $value = $x_value;
 853         }
 854
 855         $value[] = 0; // just in case the carry adds an extra digit
 856
 857         $carry = 0;
 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;
 862
 863             $temp = (int) ($sum / 0x4000000);
 864
 865             $value[$i] = (int) ($sum - 0x4000000 * $temp); // eg. a faster alternative to fmod($sum, 0x4000000)
 866             $value[$j] = $temp;
 867         }
 868
 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]
 874         }
 875
 876         if ($carry) {
 877             for (; $value[$i] == 0x3FFFFFF; ++$i) {
 878                 $value[$i] = 0;
 879             }
 880             ++$value[$i];
 881         }
 882
 883         return array(
 884             MATH_BIGINTEGER_VALUE => $this->_trim($value),
 885             MATH_BIGINTEGER_SIGN => $x_negative
 886         );
 887     }
 888
 889     /**
 890      * Subtracts two BigIntegers.
 891      *
 892      * Here's an example:
 893      * <code>
 894      * <?php
 895      *    include('Math/BigInteger.php');
 896      *
 897      *    $a = new Math_BigInteger('10');
 898      *    $b = new Math_BigInteger('20');
 899      *
 900      *    $c = $a->subtract($b);
 901      *
 902      *    echo $c->toString(); // outputs -10
 903      * ?>
 904      * </code>
 905      *
 906      * @param Math_BigInteger $y
 907      * @return Math_BigInteger
 908      * @access public
 909      * @internal Performs base-2**52 subtraction
 910      */
 911     function subtract($y)
 912     {
 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);
 917
 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);
 922
 923                 return $this->_normalize($temp);
 924         }
 925
 926         $temp = $this->_subtract($this->value, $this->is_negative, $y->value, $y->is_negative);
 927
 928         $result = new Math_BigInteger();
 929         $result->value = $temp[MATH_BIGINTEGER_VALUE];
 930         $result->is_negative = $temp[MATH_BIGINTEGER_SIGN];
 931
 932         return $this->_normalize($result);
 933     }
 934
 935     /**
 936      * Performs subtraction.
 937      *
 938      * @param Array $x_value
 939      * @param Boolean $x_negative
 940      * @param Array $y_value
 941      * @param Boolean $y_negative
 942      * @return Array
 943      * @access private
 944      */
 945     function _subtract($x_value, $x_negative, $y_value, $y_negative)
 946     {
 947         $x_size = count($x_value);
 948         $y_size = count($y_value);
 949
 950         if ($x_size == 0) {
 951             return array(
 952                 MATH_BIGINTEGER_VALUE => $y_value,
 953                 MATH_BIGINTEGER_SIGN => !$y_negative
 954             );
 955         } else if ($y_size == 0) {
 956             return array(
 957                 MATH_BIGINTEGER_VALUE => $x_value,
 958                 MATH_BIGINTEGER_SIGN => $x_negative
 959             );
 960         }
 961
 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;
 966
 967             return $temp;
 968         }
 969
 970         $diff = $this->_compare($x_value, $x_negative, $y_value, $y_negative);
 971
 972         if ( !$diff ) {
 973             return array(
 974                 MATH_BIGINTEGER_VALUE => array(),
 975                 MATH_BIGINTEGER_SIGN => false
 976             );
 977         }
 978
 979         // switch $x and $y around, if appropriate.
 980         if ( (!$x_negative && $diff < 0) || ($x_negative && $diff > 0) ) {
 981             $temp = $x_value;
 982             $x_value = $y_value;
 983             $y_value = $temp;
 984
 985             $x_negative = !$x_negative;
 986
 987             $x_size = count($x_value);
 988             $y_size = count($y_value);
 989         }
 990
 991         // at this point, $x_value should be at least as big as - if not bigger than - $y_value
 992
 993         $carry = 0;
 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;
 998
 999             $temp = (int) ($sum / 0x4000000);
1000
1001             $x_value[$i] = (int) ($sum - 0x4000000 * $temp);
1002             $x_value[$j] = $temp;
1003         }
1004
1005         if ($j == $y_size) { // ie. if $y_size is odd
1006             $sum = $x_value[$i] - $y_value[$i] - $carry;
1007             $carry = $sum < 0;
1008             $x_value[$i] = $carry ? $sum + 0x4000000 : $sum;
1009             ++$i;
1010         }
1011
1012         if ($carry) {
1013             for (; !$x_value[$i]; ++$i) {
1014                 $x_value[$i] = 0x3FFFFFF;
1015             }
1016             --$x_value[$i];
1017         }
1018
1019         return array(
1020             MATH_BIGINTEGER_VALUE => $this->_trim($x_value),
1021             MATH_BIGINTEGER_SIGN => $x_negative
1022         );
1023     }
1024
1025     /**
1026      * Multiplies two BigIntegers
1027      *
1028      * Here's an example:
1029      * <code>
1030      * <?php
1031      *    include('Math/BigInteger.php');
1032      *
1033      *    $a = new Math_BigInteger('10');
1034      *    $b = new Math_BigInteger('20');
1035      *
1036      *    $c = $a->multiply($b);
1037      *
1038      *    echo $c->toString(); // outputs 200
1039      * ?>
1040      * </code>
1041      *
1042      * @param Math_BigInteger $x
1043      * @return Math_BigInteger
1044      * @access public
1045      */
1046     function multiply($x)
1047     {
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);
1052
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);
1057
1058                 return $this->_normalize($temp);
1059         }
1060
1061         $temp = $this->_multiply($this->value, $this->is_negative, $x->value, $x->is_negative);
1062
1063         $product = new Math_BigInteger();
1064         $product->value = $temp[MATH_BIGINTEGER_VALUE];
1065         $product->is_negative = $temp[MATH_BIGINTEGER_SIGN];
1066
1067         return $this->_normalize($product);
1068     }
1069
1070     /**
1071      * Performs multiplication.
1072      *
1073      * @param Array $x_value
1074      * @param Boolean $x_negative
1075      * @param Array $y_value
1076      * @param Boolean $y_negative
1077      * @return Array
1078      * @access private
1079      */
1080     function _multiply($x_value, $x_negative, $y_value, $y_negative)
1081     {
1082         //if ( $x_value == $y_value ) {
1083         //    return array(
1084         //        MATH_BIGINTEGER_VALUE => $this->_square($x_value),
1085         //        MATH_BIGINTEGER_SIGN => $x_sign != $y_value
1086         //    );
1087         //}
1088
1089         $x_length = count($x_value);
1090         $y_length = count($y_value);
1091
1092         if ( !$x_length || !$y_length ) { // a 0 is being multiplied
1093             return array(
1094                 MATH_BIGINTEGER_VALUE => array(),
1095                 MATH_BIGINTEGER_SIGN => false
1096             );
1097         }
1098
1099         return array(
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
1104         );
1105     }
1106
1107     /**
1108      * Performs long multiplication on two BigIntegers
1109      *
1110      * Modeled after 'multiply' in MutableBigInteger.java.
1111      *
1112      * @param Array $x_value
1113      * @param Array $y_value
1114      * @return Array
1115      * @access private
1116      */
1117     function _regularMultiply($x_value, $y_value)
1118     {
1119         $x_length = count($x_value);
1120         $y_length = count($y_value);
1121
1122         if ( !$x_length || !$y_length ) { // a 0 is being multiplied
1123             return array();
1124         }
1125
1126         if ( $x_length < $y_length ) {
1127             $temp = $x_value;
1128             $x_value = $y_value;
1129             $y_value = $temp;
1130
1131             $x_length = count($x_value);
1132             $y_length = count($y_value);
1133         }
1134
1135         $product_value = $this->_array_repeat(0, $x_length + $y_length);
1136
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
1141         // to always be 0
1142
1143         $carry = 0;
1144
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);
1149         }
1150
1151         $product_value[$j] = $carry;
1152
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) {
1156             $carry = 0;
1157
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);
1162             }
1163
1164             $product_value[$k] = $carry;
1165         }
1166
1167         return $product_value;
1168     }
1169
1170     /**
1171      * Performs Karatsuba multiplication on two BigIntegers
1172      *
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}.
1175      *
1176      * @param Array $x_value
1177      * @param Array $y_value
1178      * @return Array
1179      * @access private
1180      */
1181     function _karatsuba($x_value, $y_value)
1182     {
1183         $m = min(count($x_value) >> 1, count($y_value) >> 1);
1184
1185         if ($m < MATH_BIGINTEGER_KARATSUBA_CUTOFF) {
1186             return $this->_regularMultiply($x_value, $y_value);
1187         }
1188
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);
1193
1194         $z2 = $this->_karatsuba($x1, $y1);
1195         $z0 = $this->_karatsuba($x0, $y0);
1196
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);
1202
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]);
1205
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);
1208
1209         return $xy[MATH_BIGINTEGER_VALUE];
1210     }
1211
1212     /**
1213      * Performs squaring
1214      *
1215      * @param Array $x
1216      * @return Array
1217      * @access private
1218      */
1219     function _square($x = false)
1220     {
1221         return count($x) < 2 * MATH_BIGINTEGER_KARATSUBA_CUTOFF ?
1222             $this->_trim($this->_baseSquare($x)) :
1223             $this->_trim($this->_karatsubaSquare($x));
1224     }
1225
1226     /**
1227      * Performs traditional squaring on two BigIntegers
1228      *
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.
1232      *
1233      * @param Array $value
1234      * @return Array
1235      * @access private
1236      */
1237     function _baseSquare($value)
1238     {
1239         if ( empty($value) ) {
1240             return array();
1241         }
1242         $square_value = $this->_array_repeat(0, 2 * count($value));
1243
1244         for ($i = 0, $max_index = count($value) - 1; $i <= $max_index; ++$i) {
1245             $i2 = $i << 1;
1246
1247             $temp = $square_value[$i2] + $value[$i] * $value[$i];
1248             $carry = (int) ($temp / 0x4000000);
1249             $square_value[$i2] = (int) ($temp - 0x4000000 * $carry);
1250
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);
1256             }
1257
1258             // the following line can yield values larger 2**15.  at this point, PHP should switch
1259             // over to floats.
1260             $square_value[$i + $max_index + 1] = $carry;
1261         }
1262
1263         return $square_value;
1264     }
1265
1266     /**
1267      * Performs Karatsuba "squaring" on two BigIntegers
1268      *
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}.
1271      *
1272      * @param Array $value
1273      * @return Array
1274      * @access private
1275      */
1276     function _karatsubaSquare($value)
1277     {
1278         $m = count($value) >> 1;
1279
1280         if ($m < MATH_BIGINTEGER_KARATSUBA_CUTOFF) {
1281             return $this->_baseSquare($value);
1282         }
1283
1284         $x1 = array_slice($value, $m);
1285         $x0 = array_slice($value, 0, $m);
1286
1287         $z2 = $this->_karatsubaSquare($x1);
1288         $z0 = $this->_karatsubaSquare($x0);
1289
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);
1294
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]);
1297
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);
1300
1301         return $xx[MATH_BIGINTEGER_VALUE];
1302     }
1303
1304     /**
1305      * Divides two BigIntegers.
1306      *
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).
1311      *
1312      * Here's an example:
1313      * <code>
1314      * <?php
1315      *    include('Math/BigInteger.php');
1316      *
1317      *    $a = new Math_BigInteger('10');
1318      *    $b = new Math_BigInteger('20');
1319      *
1320      *    list($quotient, $remainder) = $a->divide($b);
1321      *
1322      *    echo $quotient->toString(); // outputs 0
1323      *    echo "\r\n";
1324      *    echo $remainder->toString(); // outputs 10
1325      * ?>
1326      * </code>
1327      *
1328      * @param Math_BigInteger $y
1329      * @return Array
1330      * @access public
1331      * @internal This function is based off of {@link http://www.cacr.math.uwaterloo.ca/hac/about/chap14.pdf#page=9 HAC 14.20}.
1332      */
1333     function divide($y)
1334     {
1335         switch ( MATH_BIGINTEGER_MODE ) {
1336             case MATH_BIGINTEGER_MODE_GMP:
1337                 $quotient = new Math_BigInteger();
1338                 $remainder = new Math_BigInteger();
1339
1340                 list($quotient->value, $remainder->value) = gmp_div_qr($this->value, $y->value);
1341
1342                 if (gmp_sign($remainder->value) < 0) {
1343                     $remainder->value = gmp_add($remainder->value, gmp_abs($y->value));
1344                 }
1345
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();
1350
1351                 $quotient->value = bcdiv($this->value, $y->value, 0);
1352                 $remainder->value = bcmod($this->value, $y->value);
1353
1354                 if ($remainder->value[0] == '-') {
1355                     $remainder->value = bcadd($remainder->value, $y->value[0] == '-' ? substr($y->value, 1) : $y->value, 0);
1356                 }
1357
1358                 return array($this->_normalize($quotient), $this->_normalize($remainder));
1359         }
1360
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));
1369         }
1370
1371         static $zero;
1372         if ( !isset($zero) ) {
1373             $zero = new Math_BigInteger();
1374         }
1375
1376         $x = $this->copy();
1377         $y = $y->copy();
1378
1379         $x_sign = $x->is_negative;
1380         $y_sign = $y->is_negative;
1381
1382         $x->is_negative = $y->is_negative = false;
1383
1384         $diff = $x->compare($y);
1385
1386         if ( !$diff ) {
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()));
1391         }
1392
1393         if ( $diff < 0 ) {
1394             // if $x is negative, "add" $y.
1395             if ( $x_sign ) {
1396                 $x = $y->subtract($x);
1397             }
1398             return array($this->_normalize(new Math_BigInteger()), $this->_normalize($x));
1399         }
1400
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) {
1404             $msb <<= 1;
1405         }
1406         $x->_lshift($shift);
1407         $y->_lshift($shift);
1408         $y_value = &$y->value;
1409
1410         $x_max = count($x->value) - 1;
1411         $y_max = count($y->value) - 1;
1412
1413         $quotient = new Math_BigInteger();
1414         $quotient_value = &$quotient->value;
1415         $quotient_value = $this->_array_repeat(0, $x_max - $y_max + 1);
1416
1417         static $temp, $lhs, $rhs;
1418         if (!isset($temp)) {
1419             $temp = new Math_BigInteger();
1420             $lhs =  new Math_BigInteger();
1421             $rhs =  new Math_BigInteger();
1422         }
1423         $temp_value = &$temp->value;
1424         $rhs_value =  &$rhs->value;
1425
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);
1428
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;
1434         }
1435
1436         for ($i = $x_max; $i >= $y_max + 1; --$i) {
1437             $x_value = &$x->value;
1438             $x_window = array(
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
1442             );
1443             $y_window = array(
1444                 $y_value[$y_max],
1445                 ( $y_max > 0 ) ? $y_value[$y_max - 1] : 0
1446             );
1447
1448             $q_index = $i - $y_max - 1;
1449             if ($x_window[0] == $y_window[0]) {
1450                 $quotient_value[$q_index] = 0x3FFFFFF;
1451             } else {
1452                 $quotient_value[$q_index] = (int) (
1453                     ($x_window[0] * 0x4000000 + $x_window[1])
1454                     /
1455                     $y_window[0]
1456                 );
1457             }
1458
1459             $temp_value = array($y_window[1], $y_window[0]);
1460
1461             $lhs->value = array($quotient_value[$q_index]);
1462             $lhs = $lhs->multiply($temp);
1463
1464             $rhs_value = array($x_window[2], $x_window[1], $x_window[0]);
1465
1466             while ( $lhs->compare($rhs) > 0 ) {
1467                 --$quotient_value[$q_index];
1468
1469                 $lhs->value = array($quotient_value[$q_index]);
1470                 $lhs = $lhs->multiply($temp);
1471             }
1472
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);
1478
1479             $x = $x->subtract($temp);
1480
1481             if ($x->compare($zero) < 0) {
1482                 $temp_value = array_merge($adjust, $y_value);
1483                 $x = $x->add($temp);
1484
1485                 --$quotient_value[$q_index];
1486             }
1487
1488             $x_max = count($x_value) - 1;
1489         }
1490
1491         // unnormalize the remainder
1492         $x->_rshift($shift);
1493
1494         $quotient->is_negative = $x_sign != $y_sign;
1495
1496         // calculate the "common residue", if appropriate
1497         if ( $x_sign ) {
1498             $y->_rshift($shift);
1499             $x = $y->subtract($x);
1500         }
1501
1502         return array($this->_normalize($quotient), $this->_normalize($x));
1503     }
1504
1505     /**
1506      * Divides a BigInteger by a regular integer
1507      *
1508      * abc / x = a00 / x + b0 / x + c / x
1509      *
1510      * @param Array $dividend
1511      * @param Array $divisor
1512      * @return Array
1513      * @access private
1514      */
1515     function _divide_digit($dividend, $divisor)
1516     {
1517         $carry = 0;
1518         $result = array();
1519
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]);
1524         }
1525
1526         return array($result, $carry);
1527     }
1528
1529     /**
1530      * Performs modular exponentiation.
1531      *
1532      * Here's an example:
1533      * <code>
1534      * <?php
1535      *    include('Math/BigInteger.php');
1536      *
1537      *    $a = new Math_BigInteger('10');
1538      *    $b = new Math_BigInteger('20');
1539      *    $c = new Math_BigInteger('30');
1540      *
1541      *    $c = $a->modPow($b, $c);
1542      *
1543      *    echo $c->toString(); // outputs 10
1544      * ?>
1545      * </code>
1546      *
1547      * @param Math_BigInteger $e
1548      * @param Math_BigInteger $n
1549      * @return Math_BigInteger
1550      * @access public
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.
1555      *
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.
1558      *
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?
1563      *
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.
1570      */
1571     function modPow($e, $n)
1572     {
1573         $n = $this->bitmask !== false && $this->bitmask->compare($n) < 0 ? $this->bitmask : $n->abs();
1574
1575         if ($e->compare(new Math_BigInteger()) < 0) {
1576             $e = $e->abs();
1577
1578             $temp = $this->modInverse($n);
1579             if ($temp === false) {
1580                 return false;
1581             }
1582
1583             return $this->_normalize($temp->modPow($e, $n));
1584         }
1585
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);
1590
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);
1595
1596                 return $this->_normalize($temp);
1597         }
1598
1599         if ( empty($e->value) ) {
1600             $temp = new Math_BigInteger();
1601             $temp->value = array(1);
1602             return $this->_normalize($temp);
1603         }
1604
1605         if ( $e->value == array(1) ) {
1606             list(, $temp) = $this->divide($n);
1607             return $this->_normalize($temp);
1608         }
1609
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);
1615         }
1616
1617         return $this->_normalize($this->_slidingWindow($e, $n, MATH_BIGINTEGER_BARRETT));
1618
1619         // is the modulo odd?
1620         if ( $n->value[0] & 1 ) {
1621             return $this->_normalize($this->_slidingWindow($e, $n, MATH_BIGINTEGER_MONTGOMERY));
1622         }
1623         // if it's not, it's even
1624
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;
1630                 $j+= 26 * $i;
1631                 break;
1632             }
1633         }
1634         // at this point, 2^$j * $n/(2^$j) == $n
1635
1636         $mod1 = $n->copy();
1637         $mod1->_rshift($j);
1638         $mod2 = new Math_BigInteger();
1639         $mod2->value = array(1);
1640         $mod2->_lshift($j);
1641
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);
1644
1645         $y1 = $mod2->modInverse($mod1);
1646         $y2 = $mod1->modInverse($mod2);
1647
1648         $result = $part1->multiply($mod2);
1649         $result = $result->multiply($y1);
1650
1651         $temp = $part2->multiply($mod1);
1652         $temp = $temp->multiply($y2);
1653
1654         $result = $result->add($temp);
1655         list(, $result) = $result->divide($n);
1656
1657         return $this->_normalize($result);
1658     }
1659
1660     /**
1661      * Performs modular exponentiation.
1662      *
1663      * Alias for Math_BigInteger::modPow()
1664      *
1665      * @param Math_BigInteger $e
1666      * @param Math_BigInteger $n
1667      * @return Math_BigInteger
1668      * @access public
1669      */
1670     function powMod($e, $n)
1671     {
1672         return $this->modPow($e, $n);
1673     }
1674
1675     /**
1676      * Sliding Window k-ary Modular Exponentiation
1677      *
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.
1682      *
1683      * @param Math_BigInteger $e
1684      * @param Math_BigInteger $n
1685      * @param Integer $mode
1686      * @return Math_BigInteger
1687      * @access private
1688      */
1689     function _slidingWindow($e, $n, $mode)
1690     {
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
1693
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);
1699         }
1700
1701         $e_length = strlen($e_bits);
1702
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);
1706
1707         $n_value = $n->value;
1708
1709         // precompute $this^0 through $this^$window_size
1710         $powers = array();
1711         $powers[1] = $this->_prepareReduce($this->value, $n_value, $mode);
1712         $powers[2] = $this->_squareReduce($powers[1], $n_value, $mode);
1713
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) {
1718             $i2 = $i << 1;
1719             $powers[$i2 + 1] = $this->_multiplyReduce($powers[$i2 - 1], $powers[2], $n_value, $mode);
1720         }
1721
1722         $result = array(1);
1723         $result = $this->_prepareReduce($result, $n_value, $mode);
1724
1725         for ($i = 0; $i < $e_length; ) {
1726             if ( !$e_bits[$i] ) {
1727                 $result = $this->_squareReduce($result, $n_value, $mode);
1728                 ++$i;
1729             } else {
1730                 for ($j = $window_size - 1; $j > 0; --$j) {
1731                     if ( !empty($e_bits[$i + $j]) ) {
1732                         break;
1733                     }
1734                 }
1735
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);
1738                 }
1739
1740                 $result = $this->_multiplyReduce($result, $powers[bindec(substr($e_bits, $i, $j + 1))], $n_value, $mode);
1741
1742                 $i+=$j + 1;
1743             }
1744         }
1745
1746         $temp = new Math_BigInteger();
1747         $temp->value = $this->_reduce($result, $n_value, $mode);
1748
1749         return $temp;
1750     }
1751
1752     /**
1753      * Modular reduction
1754      *
1755      * For most $modes this will return the remainder.
1756      *
1757      * @see _slidingWindow()
1758      * @access private
1759      * @param Array $x
1760      * @param Array $n
1761      * @param Integer $mode
1762      * @return Array
1763      */
1764     function _reduce($x, $n, $mode)
1765     {
1766         switch ($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();
1773                 $lhs->value = $x;
1774                 $rhs = new Math_BigInteger();
1775                 $rhs->value = $n;
1776                 return $x->_mod2($n);
1777             case MATH_BIGINTEGER_CLASSIC:
1778                 $lhs = new Math_BigInteger();
1779                 $lhs->value = $x;
1780                 $rhs = new Math_BigInteger();
1781                 $rhs->value = $n;
1782                 list(, $temp) = $lhs->divide($rhs);
1783                 return $temp->value;
1784             case MATH_BIGINTEGER_NONE:
1785                 return $x;
1786             default:
1787                 // an invalid $mode was provided
1788         }
1789     }
1790
1791     /**
1792      * Modular reduction preperation
1793      *
1794      * @see _slidingWindow()
1795      * @access private
1796      * @param Array $x
1797      * @param Array $n
1798      * @param Integer $mode
1799      * @return Array
1800      */
1801     function _prepareReduce($x, $n, $mode)
1802     {
1803         if ($mode == MATH_BIGINTEGER_MONTGOMERY) {
1804             return $this->_prepMontgomery($x, $n);
1805         }
1806         return $this->_reduce($x, $n, $mode);
1807     }
1808
1809     /**
1810      * Modular multiply
1811      *
1812      * @see _slidingWindow()
1813      * @access private
1814      * @param Array $x
1815      * @param Array $y
1816      * @param Array $n
1817      * @param Integer $mode
1818      * @return Array
1819      */
1820     function _multiplyReduce($x, $y, $n, $mode)
1821     {
1822         if ($mode == MATH_BIGINTEGER_MONTGOMERY) {
1823             return $this->_montgomeryMultiply($x, $y, $n);
1824         }
1825         $temp = $this->_multiply($x, false, $y, false);
1826         return $this->_reduce($temp[MATH_BIGINTEGER_VALUE], $n, $mode);
1827     }
1828
1829     /**
1830      * Modular square
1831      *
1832      * @see _slidingWindow()
1833      * @access private
1834      * @param Array $x
1835      * @param Array $n
1836      * @param Integer $mode
1837      * @return Array
1838      */
1839     function _squareReduce($x, $n, $mode)
1840     {
1841         if ($mode == MATH_BIGINTEGER_MONTGOMERY) {
1842             return $this->_montgomeryMultiply($x, $x, $n);
1843         }
1844         return $this->_reduce($this->_square($x), $n, $mode);
1845     }
1846
1847     /**
1848      * Modulos for Powers of Two
1849      *
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.
1852      *
1853      * @see _slidingWindow()
1854      * @access private
1855      * @param Math_BigInteger
1856      * @return Math_BigInteger
1857      */
1858     function _mod2($n)
1859     {
1860         $temp = new Math_BigInteger();
1861         $temp->value = array(1);
1862         return $this->bitwise_and($n->subtract($temp));
1863     }
1864
1865     /**
1866      * Barrett Modular Reduction
1867      *
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).
1871      *
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."
1875      *
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.
1882      *
1883      * @see _slidingWindow()
1884      * @access private
1885      * @param Array $n
1886      * @param Array $m
1887      * @return Array
1888      */
1889     function _barrett($n, $m)
1890     {
1891         static $cache = array(
1892             MATH_BIGINTEGER_VARIABLE => array(),
1893             MATH_BIGINTEGER_DATA => array()
1894         );
1895
1896         $m_length = count($m);
1897
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();
1902             $lhs->value = $n;
1903             $rhs->value = $m;
1904             list(, $temp) = $lhs->divide($rhs);
1905             return $temp->value;
1906         }
1907
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);
1911         }
1912
1913         // n = 2 * m.length
1914
1915         if ( ($key = array_search($m, $cache[MATH_BIGINTEGER_VARIABLE])) === false ) {
1916             $key = count($cache[MATH_BIGINTEGER_VARIABLE]);
1917             $cache[MATH_BIGINTEGER_VARIABLE][] = $m;
1918
1919             $lhs = new Math_BigInteger();
1920             $lhs_value = &$lhs->value;
1921             $lhs_value = $this->_array_repeat(0, $m_length + ($m_length >> 1));
1922             $lhs_value[] = 1;
1923             $rhs = new Math_BigInteger();
1924             $rhs->value = $m;
1925
1926             list($u, $m1) = $lhs->divide($rhs);
1927             $u = $u->value;
1928             $m1 = $m1->value;
1929
1930             $cache[MATH_BIGINTEGER_DATA][] = array(
1931                 'u' => $u, // m.length >> 1 (technically (m.length >> 1) + 1)
1932                 'm1'=> $m1 // m.length
1933             );
1934         } else {
1935             extract($cache[MATH_BIGINTEGER_DATA][$key]);
1936         }
1937
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
1944
1945         if ($m_length & 1) {
1946             return $this->_regularBarrett($n[MATH_BIGINTEGER_VALUE], $m);
1947         }
1948
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);
1960
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).
1964
1965         $result = $this->_subtract($n[MATH_BIGINTEGER_VALUE], false, $temp[MATH_BIGINTEGER_VALUE], false);
1966
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);
1969         }
1970
1971         return $result[MATH_BIGINTEGER_VALUE];
1972     }
1973
1974     /**
1975      * (Regular) Barrett Modular Reduction
1976      *
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.
1979      *
1980      * @see _slidingWindow()
1981      * @access private
1982      * @param Array $x
1983      * @param Array $n
1984      * @return Array
1985      */
1986     function _regularBarrett($x, $n)
1987     {
1988         static $cache = array(
1989             MATH_BIGINTEGER_VARIABLE => array(),
1990             MATH_BIGINTEGER_DATA => array()
1991         );
1992
1993         $n_length = count($n);
1994
1995         if (count($x) > 2 * $n_length) {
1996             $lhs = new Math_BigInteger();
1997             $rhs = new Math_BigInteger();
1998             $lhs->value = $x;
1999             $rhs->value = $n;
2000             list(, $temp) = $lhs->divide($rhs);
2001             return $temp->value;
2002         }
2003
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);
2010             $lhs_value[] = 1;
2011             $rhs = new Math_BigInteger();
2012             $rhs->value = $n;
2013             list($temp, ) = $lhs->divide($rhs); // m.length
2014             $cache[MATH_BIGINTEGER_DATA][] = $temp->value;
2015         }
2016
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);
2023
2024         // m.length + 1
2025         $result = array_slice($x, 0, $n_length + 1);
2026         // m.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)
2029
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];
2035         }
2036
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);
2041         }
2042
2043         return $result[MATH_BIGINTEGER_VALUE];
2044     }
2045
2046     /**
2047      * Performs long multiplication up to $stop digits
2048      *
2049      * If you're going to be doing array_slice($product->value, 0, $stop), some cycles can be saved.
2050      *
2051      * @see _regularBarrett()
2052      * @param Array $x_value
2053      * @param Boolean $x_negative
2054      * @param Array $y_value
2055      * @param Boolean $y_negative
2056      * @return Array
2057      * @access private
2058      */
2059     function _multiplyLower($x_value, $x_negative, $y_value, $y_negative, $stop)
2060     {
2061         $x_length = count($x_value);
2062         $y_length = count($y_value);
2063
2064         if ( !$x_length || !$y_length ) { // a 0 is being multiplied
2065             return array(
2066                 MATH_BIGINTEGER_VALUE => array(),
2067                 MATH_BIGINTEGER_SIGN => false
2068             );
2069         }
2070
2071         if ( $x_length < $y_length ) {
2072             $temp = $x_value;
2073             $x_value = $y_value;
2074             $y_value = $temp;
2075
2076             $x_length = count($x_value);
2077             $y_length = count($y_value);
2078         }
2079
2080         $product_value = $this->_array_repeat(0, $x_length + $y_length);
2081
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
2086         // to always be 0
2087
2088         $carry = 0;
2089
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);
2094         }
2095
2096         if ($j < $stop) {
2097             $product_value[$j] = $carry;
2098         }
2099
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"
2102
2103         for ($i = 1; $i < $y_length; ++$i) {
2104             $carry = 0;
2105
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);
2110             }
2111
2112             if ($k < $stop) {
2113                 $product_value[$k] = $carry;
2114             }
2115         }
2116
2117         return array(
2118             MATH_BIGINTEGER_VALUE => $this->_trim($product_value),
2119             MATH_BIGINTEGER_SIGN => $x_negative != $y_negative
2120         );
2121     }
2122
2123     /**
2124      * Montgomery Modular Reduction
2125      *
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.
2130      *
2131      * @see _prepMontgomery()
2132      * @see _slidingWindow()
2133      * @access private
2134      * @param Array $x
2135      * @param Array $n
2136      * @return Array
2137      */
2138     function _montgomery($x, $n)
2139     {
2140         static $cache = array(
2141             MATH_BIGINTEGER_VARIABLE => array(),
2142             MATH_BIGINTEGER_DATA => array()
2143         );
2144
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);
2149         }
2150
2151         $k = count($n);
2152
2153         $result = array(MATH_BIGINTEGER_VALUE => $x);
2154
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);
2161         }
2162
2163         $result[MATH_BIGINTEGER_VALUE] = array_slice($result[MATH_BIGINTEGER_VALUE], $k);
2164
2165         if ($this->_compare($result, false, $n, false) >= 0) {
2166             $result = $this->_subtract($result[MATH_BIGINTEGER_VALUE], false, $n, false);
2167         }
2168
2169         return $result[MATH_BIGINTEGER_VALUE];
2170     }
2171
2172     /**
2173      * Montgomery Multiply
2174      *
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}
2177      *
2178      * @see _prepMontgomery()
2179      * @see _montgomery()
2180      * @access private
2181      * @param Array $x
2182      * @param Array $y
2183      * @param Array $m
2184      * @return Array
2185      */
2186     function _montgomeryMultiply($x, $y, $m)
2187     {
2188         $temp = $this->_multiply($x, false, $y, false);
2189         return $this->_montgomery($temp[MATH_BIGINTEGER_VALUE], $m);
2190
2191         static $cache = array(
2192             MATH_BIGINTEGER_VARIABLE => array(),
2193             MATH_BIGINTEGER_DATA => array()
2194         );
2195
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);
2200         }
2201
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);
2215         }
2216         if ($this->_compare($a[MATH_BIGINTEGER_VALUE], false, $m, false) >= 0) {
2217             $a = $this->_subtract($a[MATH_BIGINTEGER_VALUE], false, $m, false);
2218         }
2219         return $a[MATH_BIGINTEGER_VALUE];
2220     }
2221
2222     /**
2223      * Prepare a number for use in Montgomery Modular Reductions
2224      *
2225      * @see _montgomery()
2226      * @see _slidingWindow()
2227      * @access private
2228      * @param Array $x
2229      * @param Array $n
2230      * @return Array
2231      */
2232     function _prepMontgomery($x, $n)
2233     {
2234         $lhs = new Math_BigInteger();
2235         $lhs->value = array_merge($this->_array_repeat(0, count($n)), $x);
2236         $rhs = new Math_BigInteger();
2237         $rhs->value = $n;
2238
2239         list(, $temp) = $lhs->divide($rhs);
2240         return $temp->value;
2241     }
2242
2243     /**
2244      * Modular Inverse of a number mod 2**26 (eg. 67108864)
2245      *
2246      * Based off of the bnpInvDigit function implemented and justified in the following URL:
2247      *
2248      * {@link http://www-cs-students.stanford.edu/~tjw/jsbn/jsbn.js}
2249      *
2250      * The following URL provides more info:
2251      *
2252      * {@link http://groups.google.com/group/sci.crypt/msg/7a137205c1be7d85}
2253      *
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.
2261      *
2262      * Thanks to Pedro Gimeno Fortea for input!
2263      *
2264      * @see _montgomery()
2265      * @access private
2266      * @param Array $x
2267      * @return Integer
2268      */
2269     function _modInverse67108864($x) // 2**26 == 67108864
2270     {
2271         $x = -$x[0];
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;
2278     }
2279
2280     /**
2281      * Calculates modular inverses.
2282      *
2283      * Say you have (30 mod 17 * x mod 17) mod 17 == 1.  x can be found using modular inverses.
2284      *
2285      * Here's an example:
2286      * <code>
2287      * <?php
2288      *    include('Math/BigInteger.php');
2289      *
2290      *    $a = new Math_BigInteger(30);
2291      *    $b = new Math_BigInteger(17);
2292      *
2293      *    $c = $a->modInverse($b);
2294      *    echo $c->toString(); // outputs 4
2295      *
2296      *    echo "\r\n";
2297      *
2298      *    $d = $a->multiply($c);
2299      *    list(, $d) = $d->divide($b);
2300      *    echo $d; // outputs 1 (as per the definition of modular inverse)
2301      * ?>
2302      * </code>
2303      *
2304      * @param Math_BigInteger $n
2305      * @return mixed false, if no modular inverse exists, Math_BigInteger, otherwise.
2306      * @access public
2307      * @internal See {@link http://www.cacr.math.uwaterloo.ca/hac/about/chap14.pdf#page=21 HAC 14.64} for more information.
2308      */
2309     function modInverse($n)
2310     {
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);
2315
2316                 return ( $temp->value === false ) ? false : $this->_normalize($temp);
2317         }
2318
2319         static $zero, $one;
2320         if (!isset($zero)) {
2321             $zero = new Math_BigInteger();
2322             $one = new Math_BigInteger(1);
2323         }
2324
2325         // $x mod $n == $x mod -$n.
2326         $n = $n->abs();
2327
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));
2332         }
2333
2334         extract($this->extendedGCD($n));
2335
2336         if (!$gcd->equals($one)) {
2337             return false;
2338         }
2339
2340         $x = $x->compare($zero) < 0 ? $x->add($n) : $x;
2341
2342         return $this->compare($zero) < 0 ? $this->_normalize($n->subtract($x)) : $this->_normalize($x);
2343     }
2344
2345     /**
2346      * Calculates the greatest common divisor and Bézout's identity.
2347      *
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.
2352      *
2353      * Here's an example:
2354      * <code>
2355      * <?php
2356      *    include('Math/BigInteger.php');
2357      *
2358      *    $a = new Math_BigInteger(693);
2359      *    $b = new Math_BigInteger(609);
2360      *
2361      *    extract($a->extendedGCD($b));
2362      *
2363      *    echo $gcd->toString() . "\r\n"; // outputs 21
2364      *    echo $a->toString() * $x->toString() + $b->toString() * $y->toString(); // outputs 21
2365      * ?>
2366      * </code>
2367      *
2368      * @param Math_BigInteger $n
2369      * @return Math_BigInteger
2370      * @access public
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".
2374      */
2375     function extendedGCD($n)
2376     {
2377         switch ( MATH_BIGINTEGER_MODE ) {
2378             case MATH_BIGINTEGER_MODE_GMP:
2379                 extract(gmp_gcdext($this->value, $n->value));
2380
2381                 return array(
2382                     'gcd' => $this->_normalize(new Math_BigInteger($g)),
2383                     'x'   => $this->_normalize(new Math_BigInteger($s)),
2384                     'y'   => $this->_normalize(new Math_BigInteger($t))
2385                 );
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.
2390
2391                 $u = $this->value;
2392                 $v = $n->value;
2393
2394                 $a = '1';
2395                 $b = '0';
2396                 $c = '0';
2397                 $d = '1';
2398
2399                 while (bccomp($v, '0', 0) != 0) {
2400                     $q = bcdiv($u, $v, 0);
2401
2402                     $temp = $u;
2403                     $u = $v;
2404                     $v = bcsub($temp, bcmul($v, $q, 0), 0);
2405
2406                     $temp = $a;
2407                     $a = $c;
2408                     $c = bcsub($temp, bcmul($a, $q, 0), 0);
2409
2410                     $temp = $b;
2411                     $b = $d;
2412                     $d = bcsub($temp, bcmul($b, $q, 0), 0);
2413                 }
2414
2415                 return array(
2416                     'gcd' => $this->_normalize(new Math_BigInteger($u)),
2417                     'x'   => $this->_normalize(new Math_BigInteger($a)),
2418                     'y'   => $this->_normalize(new Math_BigInteger($b))
2419                 );
2420         }
2421
2422         $y = $n->copy();
2423         $x = $this->copy();
2424         $g = new Math_BigInteger();
2425         $g->value = array(1);
2426
2427         while ( !(($x->value[0] & 1)|| ($y->value[0] & 1)) ) {
2428             $x->_rshift(1);
2429             $y->_rshift(1);
2430             $g->_lshift(1);
2431         }
2432
2433         $u = $x->copy();
2434         $v = $y->copy();
2435
2436         $a = new Math_BigInteger();
2437         $b = new Math_BigInteger();
2438         $c = new Math_BigInteger();
2439         $d = new Math_BigInteger();
2440
2441         $a->value = $d->value = $g->value = array(1);
2442         $b->value = $c->value = array();
2443
2444         while ( !empty($u->value) ) {
2445             while ( !($u->value[0] & 1) ) {
2446                 $u->_rshift(1);
2447                 if ( (!empty($a->value) && ($a->value[0] & 1)) || (!empty($b->value) && ($b->value[0] & 1)) ) {
2448                     $a = $a->add($y);
2449                     $b = $b->subtract($x);
2450                 }
2451                 $a->_rshift(1);
2452                 $b->_rshift(1);
2453             }
2454
2455             while ( !($v->value[0] & 1) ) {
2456                 $v->_rshift(1);
2457                 if ( (!empty($d->value) && ($d->value[0] & 1)) || (!empty($c->value) && ($c->value[0] & 1)) ) {
2458                     $c = $c->add($y);
2459                     $d = $d->subtract($x);
2460                 }
2461                 $c->_rshift(1);
2462                 $d->_rshift(1);
2463             }
2464
2465             if ($u->compare($v) >= 0) {
2466                 $u = $u->subtract($v);
2467                 $a = $a->subtract($c);
2468                 $b = $b->subtract($d);
2469             } else {
2470                 $v = $v->subtract($u);
2471                 $c = $c->subtract($a);
2472                 $d = $d->subtract($b);
2473             }
2474         }
2475
2476         return array(
2477             'gcd' => $this->_normalize($g->multiply($v)),
2478             'x'   => $this->_normalize($c),
2479             'y'   => $this->_normalize($d)
2480         );
2481     }
2482
2483     /**
2484      * Calculates the greatest common divisor
2485      *
2486      * Say you have 693 and 609.  The GCD is 21.
2487      *
2488      * Here's an example:
2489      * <code>
2490      * <?php
2491      *    include('Math/BigInteger.php');
2492      *
2493      *    $a = new Math_BigInteger(693);
2494      *    $b = new Math_BigInteger(609);
2495      *
2496      *    $gcd = a->extendedGCD($b);
2497      *
2498      *    echo $gcd->toString() . "\r\n"; // outputs 21
2499      * ?>
2500      * </code>
2501      *
2502      * @param Math_BigInteger $n
2503      * @return Math_BigInteger
2504      * @access public
2505      */
2506     function gcd($n)
2507     {
2508         extract($this->extendedGCD($n));
2509         return $gcd;
2510     }
2511
2512     /**
2513      * Absolute value.
2514      *
2515      * @return Math_BigInteger
2516      * @access public
2517      */
2518     function abs()
2519     {
2520         $temp = new Math_BigInteger();
2521
2522         switch ( MATH_BIGINTEGER_MODE ) {
2523             case MATH_BIGINTEGER_MODE_GMP:
2524                 $temp->value = gmp_abs($this->value);
2525                 break;
2526             case MATH_BIGINTEGER_MODE_BCMATH:
2527                 $temp->value = (bccomp($this->value, '0', 0) < 0) ? substr($this->value, 1) : $this->value;
2528                 break;
2529             default:
2530                 $temp->value = $this->value;
2531         }
2532
2533         return $temp;
2534     }
2535
2536     /**
2537      * Compares two numbers.
2538      *
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:
2541      *
2542      * $x  > $y: $x->compare($y)  > 0
2543      * $x  < $y: $x->compare($y)  < 0
2544      * $x == $y: $x->compare($y) == 0
2545      *
2546      * Note how the same comparison operator is used.  If you want to test for equality, use $x->equals($y).
2547      *
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.
2550      * @access public
2551      * @see equals()
2552      * @internal Could return $this->subtract($x), but that's not as fast as what we do do.
2553      */
2554     function compare($y)
2555     {
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);
2561         }
2562
2563         return $this->_compare($this->value, $this->is_negative, $y->value, $y->is_negative);
2564     }
2565
2566     /**
2567      * Compares two numbers.
2568      *
2569      * @param Array $x_value
2570      * @param Boolean $x_negative
2571      * @param Array $y_value
2572      * @param Boolean $y_negative
2573      * @return Integer
2574      * @see compare()
2575      * @access private
2576      */
2577     function _compare($x_value, $x_negative, $y_value, $y_negative)
2578     {
2579         if ( $x_negative != $y_negative ) {
2580             return ( !$x_negative && $y_negative ) ? 1 : -1;
2581         }
2582
2583         $result = $x_negative ? -1 : 1;
2584
2585         if ( count($x_value) != count($y_value) ) {
2586             return ( count($x_value) > count($y_value) ) ? $result : -$result;
2587         }
2588         $size = max(count($x_value), count($y_value));
2589
2590         $x_value = array_pad($x_value, $size, 0);
2591         $y_value = array_pad($y_value, $size, 0);
2592
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;
2596             }
2597         }
2598
2599         return 0;
2600     }
2601
2602     /**
2603      * Tests the equality of two numbers.
2604      *
2605      * If you need to see if one number is greater than or less than another number, use Math_BigInteger::compare()
2606      *
2607      * @param Math_BigInteger $x
2608      * @return Boolean
2609      * @access public
2610      * @see compare()
2611      */
2612     function equals($x)
2613     {
2614         switch ( MATH_BIGINTEGER_MODE ) {
2615             case MATH_BIGINTEGER_MODE_GMP:
2616                 return gmp_cmp($this->value, $x->value) == 0;
2617             default:
2618                 return $this->value === $x->value && $this->is_negative == $x->is_negative;
2619         }
2620     }
2621
2622     /**
2623      * Set Precision
2624      *
2625      * Some bitwise operations give different results depending on the precision being used.  Examples include left
2626      * shift, not, and rotates.
2627      *
2628      * @param Math_BigInteger $x
2629      * @access public
2630      * @return Math_BigInteger
2631      */
2632     function setPrecision($bits)
2633     {
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);
2637         } else {
2638             $this->bitmask = new Math_BigInteger(bcpow('2', $bits, 0));
2639         }
2640
2641         $temp = $this->_normalize($this);
2642         $this->value = $temp->value;
2643     }
2644
2645     /**
2646      * Logical And
2647      *
2648      * @param Math_BigInteger $x
2649      * @access public
2650      * @internal Implemented per a request by Lluis Pamies i Juarez <lluis _a_ pamies.cat>
2651      * @return Math_BigInteger
2652      */
2653     function bitwise_and($x)
2654     {
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);
2659
2660                 return $this->_normalize($temp);
2661             case MATH_BIGINTEGER_MODE_BCMATH:
2662                 $left = $this->toBytes();
2663                 $right = $x->toBytes();
2664
2665                 $length = max(strlen($left), strlen($right));
2666
2667                 $left = str_pad($left, $length, chr(0), STR_PAD_LEFT);
2668                 $right = str_pad($right, $length, chr(0), STR_PAD_LEFT);
2669
2670                 return $this->_normalize(new Math_BigInteger($left & $right, 256));
2671         }
2672
2673         $result = $this->copy();
2674
2675         $length = min(count($x->value), count($this->value));
2676
2677         $result->value = array_slice($result->value, 0, $length);
2678
2679         for ($i = 0; $i < $length; ++$i) {
2680             $result->value[$i] = $result->value[$i] & $x->value[$i];
2681         }
2682
2683         return $this->_normalize($result);
2684     }
2685
2686     /**
2687      * Logical Or
2688      *
2689      * @param Math_BigInteger $x
2690      * @access public
2691      * @internal Implemented per a request by Lluis Pamies i Juarez <lluis _a_ pamies.cat>
2692      * @return Math_BigInteger
2693      */
2694     function bitwise_or($x)
2695     {
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);
2700
2701                 return $this->_normalize($temp);
2702             case MATH_BIGINTEGER_MODE_BCMATH:
2703                 $left = $this->toBytes();
2704                 $right = $x->toBytes();
2705
2706                 $length = max(strlen($left), strlen($right));
2707
2708                 $left = str_pad($left, $length, chr(0), STR_PAD_LEFT);
2709                 $right = str_pad($right, $length, chr(0), STR_PAD_LEFT);
2710
2711                 return $this->_normalize(new Math_BigInteger($left | $right, 256));
2712         }
2713
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);
2718
2719         for ($i = 0; $i < $length; ++$i) {
2720             $result->value[$i] = $this->value[$i] | $x->value[$i];
2721         }
2722
2723         return $this->_normalize($result);
2724     }
2725
2726     /**
2727      * Logical Exclusive-Or
2728      *
2729      * @param Math_BigInteger $x
2730      * @access public
2731      * @internal Implemented per a request by Lluis Pamies i Juarez <lluis _a_ pamies.cat>
2732      * @return Math_BigInteger
2733      */
2734     function bitwise_xor($x)
2735     {
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);
2740
2741                 return $this->_normalize($temp);
2742             case MATH_BIGINTEGER_MODE_BCMATH:
2743                 $left = $this->toBytes();
2744                 $right = $x->toBytes();
2745
2746                 $length = max(strlen($left), strlen($right));
2747
2748                 $left = str_pad($left, $length, chr(0), STR_PAD_LEFT);
2749                 $right = str_pad($right, $length, chr(0), STR_PAD_LEFT);
2750
2751                 return $this->_normalize(new Math_BigInteger($left ^ $right, 256));
2752         }
2753
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);
2758
2759         for ($i = 0; $i < $length; ++$i) {
2760             $result->value[$i] = $this->value[$i] ^ $x->value[$i];
2761         }
2762
2763         return $this->_normalize($result);
2764     }
2765
2766     /**
2767      * Logical Not
2768      *
2769      * @access public
2770      * @internal Implemented per a request by Lluis Pamies i Juarez <lluis _a_ pamies.cat>
2771      * @return Math_BigInteger
2772      */
2773     function bitwise_not()
2774     {
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]));
2779         $temp = ~$temp;
2780         $msb = decbin(ord($temp[0]));
2781         if (strlen($msb) == 8) {
2782             $msb = substr($msb, strpos($msb, '0'));
2783         }
2784         $temp[0] = chr(bindec($msb));
2785
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));
2791         }
2792
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);
2796
2797         $temp = str_pad($temp, ceil($this->bits / 8), chr(0), STR_PAD_LEFT);
2798
2799         return $this->_normalize(new Math_BigInteger($leading_ones | $temp, 256));
2800     }
2801
2802     /**
2803      * Logical Right Shift
2804      *
2805      * Shifts BigInteger's by $shift bits, effectively dividing by 2**$shift.
2806      *
2807      * @param Integer $shift
2808      * @return Math_BigInteger
2809      * @access public
2810      * @internal The only version that yields any speed increases is the internal version.
2811      */
2812     function bitwise_rightShift($shift)
2813     {
2814         $temp = new Math_BigInteger();
2815
2816         switch ( MATH_BIGINTEGER_MODE ) {
2817             case MATH_BIGINTEGER_MODE_GMP:
2818                 static $two;
2819
2820                 if (!isset($two)) {
2821                     $two = gmp_init('2');
2822                 }
2823
2824                 $temp->value = gmp_div_q($this->value, gmp_pow($two, $shift));
2825
2826                 break;
2827             case MATH_BIGINTEGER_MODE_BCMATH:
2828                 $temp->value = bcdiv($this->value, bcpow('2', $shift, 0), 0);
2829
2830                 break;
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);
2835         }
2836
2837         return $this->_normalize($temp);
2838     }
2839
2840     /**
2841      * Logical Left Shift
2842      *
2843      * Shifts BigInteger's by $shift bits, effectively multiplying by 2**$shift.
2844      *
2845      * @param Integer $shift
2846      * @return Math_BigInteger
2847      * @access public
2848      * @internal The only version that yields any speed increases is the internal version.
2849      */
2850     function bitwise_leftShift($shift)
2851     {
2852         $temp = new Math_BigInteger();
2853
2854         switch ( MATH_BIGINTEGER_MODE ) {
2855             case MATH_BIGINTEGER_MODE_GMP:
2856                 static $two;
2857
2858                 if (!isset($two)) {
2859                     $two = gmp_init('2');
2860                 }
2861
2862                 $temp->value = gmp_mul($this->value, gmp_pow($two, $shift));
2863
2864                 break;
2865             case MATH_BIGINTEGER_MODE_BCMATH:
2866                 $temp->value = bcmul($this->value, bcpow('2', $shift, 0), 0);
2867
2868                 break;
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);
2873         }
2874
2875         return $this->_normalize($temp);
2876     }
2877
2878     /**
2879      * Logical Left Rotate
2880      *
2881      * Instead of the top x bits being dropped they're appended to the shifted bit string.
2882      *
2883      * @param Integer $shift
2884      * @return Math_BigInteger
2885      * @access public
2886      */
2887     function bitwise_leftRotate($shift)
2888     {
2889         $bits = $this->toBytes();
2890
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();
2896             } else {
2897                 $mask = $this->bitmask->toBytes();
2898             }
2899         } else {
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);
2904         }
2905
2906         if ($shift < 0) {
2907             $shift+= $precision;
2908         }
2909         $shift%= $precision;
2910
2911         if (!$shift) {
2912             return $this->copy();
2913         }
2914
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);
2920     }
2921
2922     /**
2923      * Logical Right Rotate
2924      *
2925      * Instead of the bottom x bits being dropped they're prepended to the shifted bit string.
2926      *
2927      * @param Integer $shift
2928      * @return Math_BigInteger
2929      * @access public
2930      */
2931     function bitwise_rightRotate($shift)
2932     {
2933         return $this->bitwise_leftRotate(-$shift);
2934     }
2935
2936     /**
2937      * Set random number generator function
2938      *
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()
2942      *
2943      * If the random generating function is not explicitly set, it'll be assumed to be mt_rand().
2944      *
2945      * @see random()
2946      * @see randomPrime()
2947      * @param optional String $generator
2948      * @access public
2949      */
2950     function setRandomGenerator($generator)
2951     {
2952         $this->generator = $generator;
2953     }
2954
2955     /**
2956      * Generate a random number
2957      *
2958      * @param optional Integer $min
2959      * @param optional Integer $max
2960      * @return Math_BigInteger
2961      * @access public
2962      */
2963     function random($min = false, $max = false)
2964     {
2965         if ($min === false) {
2966             $min = new Math_BigInteger(0);
2967         }
2968
2969         if ($max === false) {
2970             $max = new Math_BigInteger(0x7FFFFFFF);
2971         }
2972
2973         $compare = $max->compare($min);
2974
2975         if (!$compare) {
2976             return $this->_normalize($min);
2977         } else if ($compare < 0) {
2978             // if $min is bigger then $max, swap $min and $max
2979             $temp = $max;
2980             $max = $min;
2981             $min = $temp;
2982         }
2983
2984         $generator = $this->generator;
2985
2986         $max = $max->subtract($min);
2987         $max = ltrim($max->toBytes(), chr(0));
2988         $size = strlen($max) - 1;
2989         $random = '';
2990
2991         $bytes = $size & 1;
2992         for ($i = 0; $i < $bytes; ++$i) {
2993             $random.= chr($generator(0, 255));
2994         }
2995
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));
3000         }
3001
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;
3005         } else {
3006             $random = chr($generator(0, ord($max[0])    )) . $random;
3007         }
3008
3009         $random = new Math_BigInteger($random, 256);
3010
3011         return $this->_normalize($random->add($min));
3012     }
3013
3014     /**
3015      * Generate a random prime number.
3016      *
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.
3019      *
3020      * @param optional Integer $min
3021      * @param optional Integer $max
3022      * @param optional Integer $timeout
3023      * @return Math_BigInteger
3024      * @access public
3025      * @internal See {@link http://www.cacr.math.uwaterloo.ca/hac/about/chap4.pdf#page=15 HAC 4.44}.
3026      */
3027     function randomPrime($min = false, $max = false, $timeout = false)
3028     {
3029         $compare = $max->compare($min);
3030
3031         if (!$compare) {
3032             return $min;
3033         } else if ($compare < 0) {
3034             // if $min is bigger then $max, swap $min and $max
3035             $temp = $max;
3036             $max = $min;
3037             $min = $temp;
3038         }
3039
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);
3047             }
3048
3049             if ($max === false) {
3050                 $max = new Math_BigInteger(0x7FFFFFFF);
3051             }
3052
3053             $x = $this->random($min, $max);
3054
3055             $x->value = gmp_nextprime($x->value);
3056
3057             if ($x->compare($max) <= 0) {
3058                 return $x;
3059             }
3060
3061             $x->value = gmp_nextprime($min->value);
3062
3063             if ($x->compare($max) <= 0) {
3064                 return $x;
3065             }
3066
3067             return false;
3068         }
3069
3070         static $one, $two;
3071         if (!isset($one)) {
3072             $one = new Math_BigInteger(1);
3073             $two = new Math_BigInteger(2);
3074         }
3075
3076         $start = time();
3077
3078         $x = $this->random($min, $max);
3079         if ($x->equals($two)) {
3080             return $x;
3081         }
3082
3083         $x->_make_odd();
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)) {
3087                 return false;
3088             }
3089             $x = $min->copy();
3090             $x->_make_odd();
3091         }
3092
3093         $initial_x = $x->copy();
3094
3095         while (true) {
3096             if ($timeout !== false && time() - $start > $timeout) {
3097                 return false;
3098             }
3099
3100             if ($x->isPrime()) {
3101                 return $x;
3102             }
3103
3104             $x = $x->add($two);
3105
3106             if ($x->compare($max) > 0) {
3107                 $x = $min->copy();
3108                 if ($x->equals($two)) {
3109                     return $x;
3110                 }
3111                 $x->_make_odd();
3112             }
3113
3114             if ($x->equals($initial_x)) {
3115                 return false;
3116             }
3117         }
3118     }
3119
3120     /**
3121      * Make the current number odd
3122      *
3123      * If the current number is odd it'll be unchanged.  If it's even, one will be added to it.
3124      *
3125      * @see randomPrime()
3126      * @access private
3127      */
3128     function _make_odd()
3129     {
3130         switch ( MATH_BIGINTEGER_MODE ) {
3131             case MATH_BIGINTEGER_MODE_GMP:
3132                 gmp_setbit($this->value, 0);
3133                 break;
3134             case MATH_BIGINTEGER_MODE_BCMATH:
3135                 if ($this->value[strlen($this->value) - 1] % 2 == 0) {
3136                     $this->value = bcadd($this->value, '1');
3137                 }
3138                 break;
3139             default:
3140                 $this->value[0] |= 1;
3141         }
3142     }
3143
3144     /**
3145      * Checks a numer to see if it's prime
3146      *
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.
3150      *
3151      * @param optional Integer $t
3152      * @return Boolean
3153      * @access public
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}.
3157      */
3158     function isPrime($t = false)
3159     {
3160         $length = strlen($this->toBytes());
3161
3162         if (!$t) {
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)
3175             else                     { $t = 27; }
3176         }
3177
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') {
3185                     return true;
3186                 }
3187                 if ($this->value[strlen($this->value) - 1] % 2 == 0) {
3188                     return false;
3189                 }
3190                 break;
3191             default:
3192                 if ($this->value == array(2)) {
3193                     return true;
3194                 }
3195                 if (~$this->value[0] & 1) {
3196                     return false;
3197                 }
3198         }
3199
3200         static $primes, $zero, $one, $two;
3201
3202         if (!isset($primes)) {
3203             $primes = array(
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
3215             );
3216
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]);
3220                 }
3221             }
3222
3223             $zero = new Math_BigInteger();
3224             $one = new Math_BigInteger(1);
3225             $two = new Math_BigInteger(2);
3226         }
3227
3228         if ($this->equals($one)) {
3229             return false;
3230         }
3231
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);
3238                 }
3239             }
3240         } else {
3241             $value = $this->value;
3242             foreach ($primes as $prime) {
3243                 list(, $r) = $this->_divide_digit($value, $prime);
3244                 if (!$r) {
3245                     return count($value) == 1 && $value[0] == $prime;
3246                 }
3247             }
3248         }
3249
3250         $n   = $this->copy();
3251         $n_1 = $n->subtract($one);
3252         $n_2 = $n->subtract($two);
3253
3254         $r = $n_1->copy();
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 ) {
3258             $s = 0;
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);
3262                 ++$s;
3263             }
3264         } else {
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);
3268                 if ($j != 25) {
3269                     break;
3270                 }
3271             }
3272             $s = 26 * $i + $j - 1;
3273             $r->_rshift($s);
3274         }
3275
3276         for ($i = 0; $i < $t; ++$i) {
3277             $a = $this->random($two, $n_2);
3278             $y = $a->modPow($r, $n);
3279
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)) {
3284                         return false;
3285                     }
3286                 }
3287
3288                 if (!$y->equals($n_1)) {
3289                     return false;
3290                 }
3291             }
3292         }
3293         return true;
3294     }
3295
3296     /**
3297      * Logical Left Shift
3298      *
3299      * Shifts BigInteger's by $shift bits.
3300      *
3301      * @param Integer $shift
3302      * @access private
3303      */
3304     function _lshift($shift)
3305     {
3306         if ( $shift == 0 ) {
3307             return;
3308         }
3309
3310         $num_digits = (int) ($shift / 26);
3311         $shift %= 26;
3312         $shift = 1 << $shift;
3313
3314         $carry = 0;
3315
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);
3320         }
3321
3322         if ( $carry ) {
3323             $this->value[] = $carry;
3324         }
3325
3326         while ($num_digits--) {
3327             array_unshift($this->value, 0);
3328         }
3329     }
3330
3331     /**
3332      * Logical Right Shift
3333      *
3334      * Shifts BigInteger's by $shift bits.
3335      *
3336      * @param Integer $shift
3337      * @access private
3338      */
3339     function _rshift($shift)
3340     {
3341         if ($shift == 0) {
3342             return;
3343         }
3344
3345         $num_digits = (int) ($shift / 26);
3346         $shift %= 26;
3347         $carry_shift = 26 - $shift;
3348         $carry_mask = (1 << $shift) - 1;
3349
3350         if ( $num_digits ) {
3351             $this->value = array_slice($this->value, $num_digits);
3352         }
3353
3354         $carry = 0;
3355
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;
3360         }
3361
3362         $this->value = $this->_trim($this->value);
3363     }
3364
3365     /**
3366      * Normalize
3367      *
3368      * Removes leading zeros and truncates (if necessary) to maintain the appropriate precision
3369      *
3370      * @param Math_BigInteger
3371      * @return Math_BigInteger
3372      * @see _trim()
3373      * @access private
3374      */
3375     function _normalize($result)
3376     {
3377         $result->precision = $this->precision;
3378         $result->bitmask = $this->bitmask;
3379
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);
3384                 }
3385
3386                 return $result;
3387             case MATH_BIGINTEGER_MODE_BCMATH:
3388                 if (!empty($result->bitmask->value)) {
3389                     $result->value = bcmod($result->value, $result->bitmask->value);
3390                 }
3391
3392                 return $result;
3393         }
3394
3395         $value = &$result->value;
3396
3397         if ( !count($value) ) {
3398             return $result;
3399         }
3400
3401         $value = $this->_trim($value);
3402
3403         if (!empty($result->bitmask->value)) {
3404             $length = min(count($value), count($this->bitmask->value));
3405             $value = array_slice($value, 0, $length);
3406
3407             for ($i = 0; $i < $length; ++$i) {
3408                 $value[$i] = $value[$i] & $this->bitmask->value[$i];
3409             }
3410         }
3411
3412         return $result;
3413     }
3414
3415     /**
3416      * Trim
3417      *
3418      * Removes leading zeros
3419      *
3420      * @return Math_BigInteger
3421      * @access private
3422      */
3423     function _trim($value)
3424     {
3425         for ($i = count($value) - 1; $i >= 0; --$i) {
3426             if ( $value[$i] ) {
3427                 break;
3428             }
3429             unset($value[$i]);
3430         }
3431
3432         return $value;
3433     }
3434
3435     /**
3436      * Array Repeat
3437      *
3438      * @param $input Array
3439      * @param $multiplier mixed
3440      * @return Array
3441      * @access private
3442      */
3443     function _array_repeat($input, $multiplier)
3444     {
3445         return ($multiplier) ? array_fill(0, $multiplier, $input) : array();
3446     }
3447
3448     /**
3449      * Logical Left Shift
3450      *
3451      * Shifts binary strings $shift bits, essentially multiplying by 2**$shift.
3452      *
3453      * @param $x String
3454      * @param $shift Integer
3455      * @return String
3456      * @access private
3457      */
3458     function _base256_lshift(&$x, $shift)
3459     {
3460         if ($shift == 0) {
3461             return;
3462         }
3463
3464         $num_bytes = $shift >> 3; // eg. floor($shift/8)
3465         $shift &= 7; // eg. $shift % 8
3466
3467         $carry = 0;
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;
3472         }
3473         $carry = ($carry != 0) ? chr($carry) : '';
3474         $x = $carry . $x . str_repeat(chr(0), $num_bytes);
3475     }
3476
3477     /**
3478      * Logical Right Shift
3479      *
3480      * Shifts binary strings $shift bits, essentially dividing by 2**$shift and returning the remainder.
3481      *
3482      * @param $x String
3483      * @param $shift Integer
3484      * @return String
3485      * @access private
3486      */
3487     function _base256_rshift(&$x, $shift)
3488     {
3489         if ($shift == 0) {
3490             $x = ltrim($x, chr(0));
3491             return '';
3492         }
3493
3494         $num_bytes = $shift >> 3; // eg. floor($shift/8)
3495         $shift &= 7; // eg. $shift % 8
3496
3497         $remainder = '';
3498         if ($num_bytes) {
3499             $start = $num_bytes > strlen($x) ? -strlen($x) : -$num_bytes;
3500             $remainder = substr($x, $start);
3501             $x = substr($x, 0, -$num_bytes);
3502         }
3503
3504         $carry = 0;
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);
3510         }
3511         $x = ltrim($x, chr(0));
3512
3513         $remainder = chr($carry >> $carry_shift) . $remainder;
3514
3515         return ltrim($remainder, chr(0));
3516     }
3517
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.
3520
3521     /**
3522      * Converts 32-bit integers to bytes.
3523      *
3524      * @param Integer $x
3525      * @return String
3526      * @access private
3527      */
3528     function _int2bytes($x)
3529     {
3530         return ltrim(pack('N', $x), chr(0));
3531     }
3532
3533     /**
3534      * Converts bytes to 32-bit integers
3535      *
3536      * @param String $x
3537      * @return Integer
3538      * @access private
3539      */
3540     function _bytes2int($x)
3541     {
3542         $temp = unpack('Nint', str_pad($x, 4, chr(0), STR_PAD_LEFT));
3543         return $temp['int'];
3544     }
3545 }