NativeCalculator.php 14 KB

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