@@ -232,15 +232,11 @@ abstract public function modPow(string $base, string $exp, string $mod): string;
232232 */
233233 public function gcd (string $ a , string $ b ): string
234234 {
235- if ($ a = == '0 ' ) {
236- return $ this ->abs ( $ b ) ;
235+ while ($ b ! == '0 ' ) {
236+ [ $ a , $ b ] = [ $ b , $ this ->divR ( $ a , $ b )] ;
237237 }
238238
239- if ($ b === '0 ' ) {
240- return $ this ->abs ($ a );
241- }
242-
243- return $ this ->gcd ($ b , $ this ->divR ($ a , $ b ));
239+ return $ this ->abs ($ a );
244240 }
245241
246242 /**
@@ -610,16 +606,23 @@ final protected function init(string $a, string $b): array
610606 */
611607 private function gcdExtended (string $ a , string $ b ): array
612608 {
613- if ($ a === '0 ' ) {
614- return [$ b , '0 ' , '1 ' ];
609+ // Iterative extended Euclidean algorithm (recursion would exhaust memory on large
610+ // inputs, see gcd()), maintaining the invariants:
611+ // $a * $x0 + $b * $y0 = $r0
612+ // $a * $x1 + $b * $y1 = $r1
613+ [$ r0 , $ r1 ] = [$ a , $ b ];
614+ [$ x0 , $ x1 ] = ['1 ' , '0 ' ];
615+ [$ y0 , $ y1 ] = ['0 ' , '1 ' ];
616+
617+ while ($ r1 !== '0 ' ) {
618+ [$ q , $ r ] = $ this ->divQR ($ r0 , $ r1 );
619+
620+ [$ r0 , $ r1 ] = [$ r1 , $ r ];
621+ [$ x0 , $ x1 ] = [$ x1 , $ this ->sub ($ x0 , $ this ->mul ($ q , $ x1 ))];
622+ [$ y0 , $ y1 ] = [$ y1 , $ this ->sub ($ y0 , $ this ->mul ($ q , $ y1 ))];
615623 }
616624
617- [$ gcd , $ x1 , $ y1 ] = $ this ->gcdExtended ($ this ->mod ($ b , $ a ), $ a );
618-
619- $ x = $ this ->sub ($ y1 , $ this ->mul ($ this ->divQ ($ b , $ a ), $ x1 ));
620- $ y = $ x1 ;
621-
622- return [$ gcd , $ x , $ y ];
625+ return [$ r0 , $ x0 , $ y0 ];
623626 }
624627
625628 /**
0 commit comments