code: def igcd_lehmer(a, b):
      r"""Computes greatest common divisor of two integers.
  
      Explanation
      ===========
  
      Euclid's algorithm for the computation of the greatest
      common divisor ``gcd(a, b)``  of two (positive) integers
      $a$ and $b$ is based on the division identity
      $$ a = q \times b + r$$,
      where the quotient  $q$  and the remainder  $r$  are integers
      and  $0 \le r < b$. Then each common divisor of  $a$  and  $b$
      divides  $r$, and it follows that  ``gcd(a, b) == gcd(b, r)``.
      The algorithm works by constructing the sequence
      r0, r1, r2, ..., where  r0 = a, r1 = b,  and each  rn
      is the remainder from the division of the two preceding
      elements.
  
      In Python, ``q = a // b``  and  ``r = a % b``  are obtained by the
      floor division and the remainder operations, respectively.
      These are the most expensive arithmetic operations, especially
      for large  a  and  b.
  
      Lehmer's algorithm [1]_ is based on the observation that the quotients
      ``qn = r(n-1) // rn``  are in general small integers even
      when  a  and  b  are very large. Hence the quotients can be
      usually determined from a relatively small number of most
      significant bits.
  
      The efficiency of the algorithm is further enhanced by not
      computing each long remainder in Euclid's sequence. The remainders
      are linear combinations of  a  and  b  with integer coefficients
      derived from the quotients. The coefficients can be computed
      as far as the quotients can be determined from the chosen
      most significant parts of  a  and  b. Only then a new pair of
      consecutive remainders is computed and the algorithm starts
      anew with this pair.
  
      References
      ==========
  
      .. [1] https://en.wikipedia.org/wiki/Lehmer%27s_GCD_algorithm
  
      """
      a, b = abs(as_int(a)), abs(as_int(b))
      if a < b:
          a, b = b, a
  
      # The algorithm works by using one or two digit division
      # whenever possible. The outer loop will replace the
      # pair (a, b) with a pair of shorter consecutive elements
      # of the Euclidean gcd sequence until a and b
      # fit into two Python (long) int digits.
      nbits = 2 * sys.int_info.bits_per_digit
  
      while a.bit_length() > nbits and b != 0:
          # Quotients are mostly small integers that can
          # be determined from most significant bits.
          n = a.bit_length() - nbits
          x, y = int(a >> n), int(b >> n)  # most significant bits
  
          # Elements of the Euclidean gcd sequence are linear
          # combinations of a and b with integer coefficients.
          # Compute the coefficients of consecutive pairs
          #     a' = A*a + B*b, b' = C*a + D*b
          # using small integer arithmetic as far as possible.
          A, B, C, D = 1, 0, 0, 1  # initial values
  
          while True:
              # The coefficients alternate in sign while looping.
              # The inner loop combines two steps to keep track
              # of the signs.
  
              # At this point we have
              #   A > 0, B <= 0, C <= 0, D > 0,
              #   x' = x + B <= x < x" = x + A,
              #   y' = y + C <= y < y" = y + D,
              # and
              #   x'*N <= a' < x"*N, y'*N <= b' < y"*N,
              # where N = 2**n.
  
              # Now, if y' > 0, and x"//y' and x'//y" agree,
              # then their common value is equal to  q = a'//b'.
              # In addition,
              #   x'%y" = x' - q*y" < x" - q*y' = x"%y',
              # and
              #   (x'%y")*N < a'%b' < (x"%y')*N.
  
              # On the other hand, we also have  x//y == q,
              # and therefore
              #   x'%y" = x + B - q*(y + D) = x%y + B',
              #   x"%y' = x + A - q*(y + C) = x%y + A',
              # where
              #    B' = B - q*D < 0, A' = A - q*C > 0.
  
              if y + C <= 0:
                  break
              q = (x + A) // (y + C)
  
              # Now  x'//y" <= q, and equality holds if
              #   x' - q*y" = (x - q*y) + (B - q*D) >= 0.
              # This is a minor optimization to avoid division.
              x_qy, B_qD = x - q * y, B - q * D
              if x_qy + B_qD < 0:
                  break
  
              # Next step in the Euclidean sequence.
              x, y = y, x_qy
              A, B, C, D = C, D, A - q * C, B_qD
  
              # At this point the signs of the coefficients
              # change and their roles are interchanged.
              #   A <= 0, B > 0, C > 0, D < 0,
              #   x' = x + A <= x < x" = x + B,
              #   y' = y + D < y < y" = y + C.
  
              if y + D <= 0:
                  break
              q = (x + B) // (y + D)
              x_qy, A_qC = x - q * y, A - q * C
              if x_qy + A_qC < 0:
                  break
  
              x, y = y, x_qy
              A, B, C, D = C, D, A_qC, B - q * D
              # Now the conditions on top of the loop
              # are again satisfied.
              #   A > 0, B < 0, C < 0, D > 0.
  
          if B == 0:
              # This can only happen when y == 0 in the beginning
              # and the inner loop does nothing.
              # Long division is forced.
              a, b = b, a % b
              continue
  
          # Compute new long arguments using the coefficients.
          a, b = A * a + B * b, C * a + D * b
  
      # Small divisors. Finish with the standard algorithm.
      while b:
          a, b = b, a % b
  
      return a
  
file_path: sympy/sympy/core/intfunc.py
imports: import math
  import sys
  from functools import lru_cache
  from .sympify import sympify
  from .singleton import S
  from sympy.external.gmpy import gcd, lcm, sqrt, iroot, bit_scan1, gcdext
  from sympy.utilities.misc import as_int, filldedent
label_desc: Use image node_id 5 for calling a global function with example usage: igcd_lehmer(a, b) and returns: a