NativeCalculator.php
446 lines
| 1 | <?php |
| 2 | |
| 3 | declare (strict_types=1); |
| 4 | namespace IAWPSCOPED\Brick\Math\Internal\Calculator; |
| 5 | |
| 6 | use IAWPSCOPED\Brick\Math\Internal\Calculator; |
| 7 | /** |
| 8 | * Calculator implementation using only native PHP code. |
| 9 | * |
| 10 | * @internal |
| 11 | * |
| 12 | * @psalm-immutable |
| 13 | */ |
| 14 | class NativeCalculator extends Calculator |
| 15 | { |
| 16 | /** |
| 17 | * The max number of digits the platform can natively add, subtract, multiply or divide without overflow. |
| 18 | * For multiplication, this represents the max sum of the lengths of both operands. |
| 19 | * |
| 20 | * In addition, it is assumed that an extra digit can hold a carry (1) without overflowing. |
| 21 | * Example: 32-bit: max number 1,999,999,999 (9 digits + carry) |
| 22 | * 64-bit: max number 1,999,999,999,999,999,999 (18 digits + carry) |
| 23 | */ |
| 24 | private int $maxDigits; |
| 25 | /** |
| 26 | * @codeCoverageIgnore |
| 27 | */ |
| 28 | public function __construct() |
| 29 | { |
| 30 | switch (\PHP_INT_SIZE) { |
| 31 | case 4: |
| 32 | $this->maxDigits = 9; |
| 33 | break; |
| 34 | case 8: |
| 35 | $this->maxDigits = 18; |
| 36 | break; |
| 37 | default: |
| 38 | throw new \RuntimeException('The platform is not 32-bit or 64-bit as expected.'); |
| 39 | } |
| 40 | } |
| 41 | public function add(string $a, string $b) : string |
| 42 | { |
| 43 | /** |
| 44 | * @psalm-var numeric-string $a |
| 45 | * @psalm-var numeric-string $b |
| 46 | */ |
| 47 | $result = $a + $b; |
| 48 | if (\is_int($result)) { |
| 49 | return (string) $result; |
| 50 | } |
| 51 | if ($a === '0') { |
| 52 | return $b; |
| 53 | } |
| 54 | if ($b === '0') { |
| 55 | return $a; |
| 56 | } |
| 57 | [$aNeg, $bNeg, $aDig, $bDig] = $this->init($a, $b); |
| 58 | $result = $aNeg === $bNeg ? $this->doAdd($aDig, $bDig) : $this->doSub($aDig, $bDig); |
| 59 | if ($aNeg) { |
| 60 | $result = $this->neg($result); |
| 61 | } |
| 62 | return $result; |
| 63 | } |
| 64 | public function sub(string $a, string $b) : string |
| 65 | { |
| 66 | return $this->add($a, $this->neg($b)); |
| 67 | } |
| 68 | public function mul(string $a, string $b) : string |
| 69 | { |
| 70 | /** |
| 71 | * @psalm-var numeric-string $a |
| 72 | * @psalm-var numeric-string $b |
| 73 | */ |
| 74 | $result = $a * $b; |
| 75 | if (\is_int($result)) { |
| 76 | return (string) $result; |
| 77 | } |
| 78 | if ($a === '0' || $b === '0') { |
| 79 | return '0'; |
| 80 | } |
| 81 | if ($a === '1') { |
| 82 | return $b; |
| 83 | } |
| 84 | if ($b === '1') { |
| 85 | return $a; |
| 86 | } |
| 87 | if ($a === '-1') { |
| 88 | return $this->neg($b); |
| 89 | } |
| 90 | if ($b === '-1') { |
| 91 | return $this->neg($a); |
| 92 | } |
| 93 | [$aNeg, $bNeg, $aDig, $bDig] = $this->init($a, $b); |
| 94 | $result = $this->doMul($aDig, $bDig); |
| 95 | if ($aNeg !== $bNeg) { |
| 96 | $result = $this->neg($result); |
| 97 | } |
| 98 | return $result; |
| 99 | } |
| 100 | public function divQ(string $a, string $b) : string |
| 101 | { |
| 102 | return $this->divQR($a, $b)[0]; |
| 103 | } |
| 104 | public function divR(string $a, string $b) : string |
| 105 | { |
| 106 | return $this->divQR($a, $b)[1]; |
| 107 | } |
| 108 | public function divQR(string $a, string $b) : array |
| 109 | { |
| 110 | if ($a === '0') { |
| 111 | return ['0', '0']; |
| 112 | } |
| 113 | if ($a === $b) { |
| 114 | return ['1', '0']; |
| 115 | } |
| 116 | if ($b === '1') { |
| 117 | return [$a, '0']; |
| 118 | } |
| 119 | if ($b === '-1') { |
| 120 | return [$this->neg($a), '0']; |
| 121 | } |
| 122 | /** @psalm-var numeric-string $a */ |
| 123 | $na = $a * 1; |
| 124 | // cast to number |
| 125 | if (\is_int($na)) { |
| 126 | /** @psalm-var numeric-string $b */ |
| 127 | $nb = $b * 1; |
| 128 | if (\is_int($nb)) { |
| 129 | // the only division that may overflow is PHP_INT_MIN / -1, |
| 130 | // which cannot happen here as we've already handled a divisor of -1 above. |
| 131 | $r = $na % $nb; |
| 132 | $q = ($na - $r) / $nb; |
| 133 | \assert(\is_int($q)); |
| 134 | return [(string) $q, (string) $r]; |
| 135 | } |
| 136 | } |
| 137 | [$aNeg, $bNeg, $aDig, $bDig] = $this->init($a, $b); |
| 138 | [$q, $r] = $this->doDiv($aDig, $bDig); |
| 139 | if ($aNeg !== $bNeg) { |
| 140 | $q = $this->neg($q); |
| 141 | } |
| 142 | if ($aNeg) { |
| 143 | $r = $this->neg($r); |
| 144 | } |
| 145 | return [$q, $r]; |
| 146 | } |
| 147 | public function pow(string $a, int $e) : string |
| 148 | { |
| 149 | if ($e === 0) { |
| 150 | return '1'; |
| 151 | } |
| 152 | if ($e === 1) { |
| 153 | return $a; |
| 154 | } |
| 155 | $odd = $e % 2; |
| 156 | $e -= $odd; |
| 157 | $aa = $this->mul($a, $a); |
| 158 | /** @psalm-suppress PossiblyInvalidArgument We're sure that $e / 2 is an int now */ |
| 159 | $result = $this->pow($aa, $e / 2); |
| 160 | if ($odd === 1) { |
| 161 | $result = $this->mul($result, $a); |
| 162 | } |
| 163 | return $result; |
| 164 | } |
| 165 | /** |
| 166 | * Algorithm from: https://www.geeksforgeeks.org/modular-exponentiation-power-in-modular-arithmetic/ |
| 167 | */ |
| 168 | public function modPow(string $base, string $exp, string $mod) : string |
| 169 | { |
| 170 | // special case: the algorithm below fails with 0 power 0 mod 1 (returns 1 instead of 0) |
| 171 | if ($base === '0' && $exp === '0' && $mod === '1') { |
| 172 | return '0'; |
| 173 | } |
| 174 | // special case: the algorithm below fails with power 0 mod 1 (returns 1 instead of 0) |
| 175 | if ($exp === '0' && $mod === '1') { |
| 176 | return '0'; |
| 177 | } |
| 178 | $x = $base; |
| 179 | $res = '1'; |
| 180 | // numbers are positive, so we can use remainder instead of modulo |
| 181 | $x = $this->divR($x, $mod); |
| 182 | while ($exp !== '0') { |
| 183 | if (\in_array($exp[-1], ['1', '3', '5', '7', '9'])) { |
| 184 | // odd |
| 185 | $res = $this->divR($this->mul($res, $x), $mod); |
| 186 | } |
| 187 | $exp = $this->divQ($exp, '2'); |
| 188 | $x = $this->divR($this->mul($x, $x), $mod); |
| 189 | } |
| 190 | return $res; |
| 191 | } |
| 192 | /** |
| 193 | * Adapted from https://cp-algorithms.com/num_methods/roots_newton.html |
| 194 | */ |
| 195 | public function sqrt(string $n) : string |
| 196 | { |
| 197 | if ($n === '0') { |
| 198 | return '0'; |
| 199 | } |
| 200 | // initial approximation |
| 201 | $x = \str_repeat('9', \intdiv(\strlen($n), 2) ?: 1); |
| 202 | $decreased = \false; |
| 203 | for (;;) { |
| 204 | $nx = $this->divQ($this->add($x, $this->divQ($n, $x)), '2'); |
| 205 | if ($x === $nx || $this->cmp($nx, $x) > 0 && $decreased) { |
| 206 | break; |
| 207 | } |
| 208 | $decreased = $this->cmp($nx, $x) < 0; |
| 209 | $x = $nx; |
| 210 | } |
| 211 | return $x; |
| 212 | } |
| 213 | /** |
| 214 | * Performs the addition of two non-signed large integers. |
| 215 | */ |
| 216 | private function doAdd(string $a, string $b) : string |
| 217 | { |
| 218 | [$a, $b, $length] = $this->pad($a, $b); |
| 219 | $carry = 0; |
| 220 | $result = ''; |
| 221 | for ($i = $length - $this->maxDigits;; $i -= $this->maxDigits) { |
| 222 | $blockLength = $this->maxDigits; |
| 223 | if ($i < 0) { |
| 224 | $blockLength += $i; |
| 225 | /** @psalm-suppress LoopInvalidation */ |
| 226 | $i = 0; |
| 227 | } |
| 228 | /** @psalm-var numeric-string $blockA */ |
| 229 | $blockA = \substr($a, $i, $blockLength); |
| 230 | /** @psalm-var numeric-string $blockB */ |
| 231 | $blockB = \substr($b, $i, $blockLength); |
| 232 | $sum = (string) ($blockA + $blockB + $carry); |
| 233 | $sumLength = \strlen($sum); |
| 234 | if ($sumLength > $blockLength) { |
| 235 | $sum = \substr($sum, 1); |
| 236 | $carry = 1; |
| 237 | } else { |
| 238 | if ($sumLength < $blockLength) { |
| 239 | $sum = \str_repeat('0', $blockLength - $sumLength) . $sum; |
| 240 | } |
| 241 | $carry = 0; |
| 242 | } |
| 243 | $result = $sum . $result; |
| 244 | if ($i === 0) { |
| 245 | break; |
| 246 | } |
| 247 | } |
| 248 | if ($carry === 1) { |
| 249 | $result = '1' . $result; |
| 250 | } |
| 251 | return $result; |
| 252 | } |
| 253 | /** |
| 254 | * Performs the subtraction of two non-signed large integers. |
| 255 | */ |
| 256 | private function doSub(string $a, string $b) : string |
| 257 | { |
| 258 | if ($a === $b) { |
| 259 | return '0'; |
| 260 | } |
| 261 | // Ensure that we always subtract to a positive result: biggest minus smallest. |
| 262 | $cmp = $this->doCmp($a, $b); |
| 263 | $invert = $cmp === -1; |
| 264 | if ($invert) { |
| 265 | $c = $a; |
| 266 | $a = $b; |
| 267 | $b = $c; |
| 268 | } |
| 269 | [$a, $b, $length] = $this->pad($a, $b); |
| 270 | $carry = 0; |
| 271 | $result = ''; |
| 272 | $complement = 10 ** $this->maxDigits; |
| 273 | for ($i = $length - $this->maxDigits;; $i -= $this->maxDigits) { |
| 274 | $blockLength = $this->maxDigits; |
| 275 | if ($i < 0) { |
| 276 | $blockLength += $i; |
| 277 | /** @psalm-suppress LoopInvalidation */ |
| 278 | $i = 0; |
| 279 | } |
| 280 | /** @psalm-var numeric-string $blockA */ |
| 281 | $blockA = \substr($a, $i, $blockLength); |
| 282 | /** @psalm-var numeric-string $blockB */ |
| 283 | $blockB = \substr($b, $i, $blockLength); |
| 284 | $sum = $blockA - $blockB - $carry; |
| 285 | if ($sum < 0) { |
| 286 | $sum += $complement; |
| 287 | $carry = 1; |
| 288 | } else { |
| 289 | $carry = 0; |
| 290 | } |
| 291 | $sum = (string) $sum; |
| 292 | $sumLength = \strlen($sum); |
| 293 | if ($sumLength < $blockLength) { |
| 294 | $sum = \str_repeat('0', $blockLength - $sumLength) . $sum; |
| 295 | } |
| 296 | $result = $sum . $result; |
| 297 | if ($i === 0) { |
| 298 | break; |
| 299 | } |
| 300 | } |
| 301 | // Carry cannot be 1 when the loop ends, as a > b |
| 302 | \assert($carry === 0); |
| 303 | $result = \ltrim($result, '0'); |
| 304 | if ($invert) { |
| 305 | $result = $this->neg($result); |
| 306 | } |
| 307 | return $result; |
| 308 | } |
| 309 | /** |
| 310 | * Performs the multiplication of two non-signed large integers. |
| 311 | */ |
| 312 | private function doMul(string $a, string $b) : string |
| 313 | { |
| 314 | $x = \strlen($a); |
| 315 | $y = \strlen($b); |
| 316 | $maxDigits = \intdiv($this->maxDigits, 2); |
| 317 | $complement = 10 ** $maxDigits; |
| 318 | $result = '0'; |
| 319 | for ($i = $x - $maxDigits;; $i -= $maxDigits) { |
| 320 | $blockALength = $maxDigits; |
| 321 | if ($i < 0) { |
| 322 | $blockALength += $i; |
| 323 | /** @psalm-suppress LoopInvalidation */ |
| 324 | $i = 0; |
| 325 | } |
| 326 | $blockA = (int) \substr($a, $i, $blockALength); |
| 327 | $line = ''; |
| 328 | $carry = 0; |
| 329 | for ($j = $y - $maxDigits;; $j -= $maxDigits) { |
| 330 | $blockBLength = $maxDigits; |
| 331 | if ($j < 0) { |
| 332 | $blockBLength += $j; |
| 333 | /** @psalm-suppress LoopInvalidation */ |
| 334 | $j = 0; |
| 335 | } |
| 336 | $blockB = (int) \substr($b, $j, $blockBLength); |
| 337 | $mul = $blockA * $blockB + $carry; |
| 338 | $value = $mul % $complement; |
| 339 | $carry = ($mul - $value) / $complement; |
| 340 | $value = (string) $value; |
| 341 | $value = \str_pad($value, $maxDigits, '0', \STR_PAD_LEFT); |
| 342 | $line = $value . $line; |
| 343 | if ($j === 0) { |
| 344 | break; |
| 345 | } |
| 346 | } |
| 347 | if ($carry !== 0) { |
| 348 | $line = $carry . $line; |
| 349 | } |
| 350 | $line = \ltrim($line, '0'); |
| 351 | if ($line !== '') { |
| 352 | $line .= \str_repeat('0', $x - $blockALength - $i); |
| 353 | $result = $this->add($result, $line); |
| 354 | } |
| 355 | if ($i === 0) { |
| 356 | break; |
| 357 | } |
| 358 | } |
| 359 | return $result; |
| 360 | } |
| 361 | /** |
| 362 | * Performs the division of two non-signed large integers. |
| 363 | * |
| 364 | * @return string[] The quotient and remainder. |
| 365 | */ |
| 366 | private function doDiv(string $a, string $b) : array |
| 367 | { |
| 368 | $cmp = $this->doCmp($a, $b); |
| 369 | if ($cmp === -1) { |
| 370 | return ['0', $a]; |
| 371 | } |
| 372 | $x = \strlen($a); |
| 373 | $y = \strlen($b); |
| 374 | // we now know that a >= b && x >= y |
| 375 | $q = '0'; |
| 376 | // quotient |
| 377 | $r = $a; |
| 378 | // remainder |
| 379 | $z = $y; |
| 380 | // focus length, always $y or $y+1 |
| 381 | for (;;) { |
| 382 | $focus = \substr($a, 0, $z); |
| 383 | $cmp = $this->doCmp($focus, $b); |
| 384 | if ($cmp === -1) { |
| 385 | if ($z === $x) { |
| 386 | // remainder < dividend |
| 387 | break; |
| 388 | } |
| 389 | $z++; |
| 390 | } |
| 391 | $zeros = \str_repeat('0', $x - $z); |
| 392 | $q = $this->add($q, '1' . $zeros); |
| 393 | $a = $this->sub($a, $b . $zeros); |
| 394 | $r = $a; |
| 395 | if ($r === '0') { |
| 396 | // remainder == 0 |
| 397 | break; |
| 398 | } |
| 399 | $x = \strlen($a); |
| 400 | if ($x < $y) { |
| 401 | // remainder < dividend |
| 402 | break; |
| 403 | } |
| 404 | $z = $y; |
| 405 | } |
| 406 | return [$q, $r]; |
| 407 | } |
| 408 | /** |
| 409 | * Compares two non-signed large numbers. |
| 410 | * |
| 411 | * @return int [-1, 0, 1] |
| 412 | */ |
| 413 | private function doCmp(string $a, string $b) : int |
| 414 | { |
| 415 | $x = \strlen($a); |
| 416 | $y = \strlen($b); |
| 417 | $cmp = $x <=> $y; |
| 418 | if ($cmp !== 0) { |
| 419 | return $cmp; |
| 420 | } |
| 421 | return \strcmp($a, $b) <=> 0; |
| 422 | // enforce [-1, 0, 1] |
| 423 | } |
| 424 | /** |
| 425 | * Pads the left of one of the given numbers with zeros if necessary to make both numbers the same length. |
| 426 | * |
| 427 | * The numbers must only consist of digits, without leading minus sign. |
| 428 | * |
| 429 | * @return array{string, string, int} |
| 430 | */ |
| 431 | private function pad(string $a, string $b) : array |
| 432 | { |
| 433 | $x = \strlen($a); |
| 434 | $y = \strlen($b); |
| 435 | if ($x > $y) { |
| 436 | $b = \str_repeat('0', $x - $y) . $b; |
| 437 | return [$a, $b, $x]; |
| 438 | } |
| 439 | if ($x < $y) { |
| 440 | $a = \str_repeat('0', $y - $x) . $a; |
| 441 | return [$a, $b, $y]; |
| 442 | } |
| 443 | return [$a, $b, $x]; |
| 444 | } |
| 445 | } |
| 446 |