[Python] GetPrimes code snippet

GetPrimes code snippet: (naive)

def GetPrimes(n):
  """
  Get all primes up to n.

  Use sieve of Eratosthenes to generate all primes less or equal to n. For
  example, given n = 10, then returns a list of primes [2, 3, 5, 7].

  Args:
    n: An integer representing the upper bound.

  Returns:
    A list containing all primes <= n in the increasing order.
  """
  sqrt_n = int(math.sqrt(n))
  visited = [False] * (n + 1)
  for i in range(2, sqrt_n + 1):
    if not visited[i]:
      for j in range(i * i, n + 1, i):
        visited[j] = True
  primes = []
  for i in range(2, n + 1):
    if not visited[i]:
      primes.append(i)
  return primes

[Python] ModInverse code snippet

ModInverse code snippet: (depends on ExtendedGcd())

def ModInverse(a, m):
  """
  Modular multiplicative inverse algorithm.

  Computes modular multiplicative inverse of an integer a is an integer x such
  that the product ax is congruent to 1 with respect to the modulus m.

  Args:
    a: An integer.
    m: An integer representing modulus.

  Returns:
    An integer x satisfying ax = 1 (mod m). Otherwise throw an exception if
    there is no such integer x.
  """
  g, x, y = ExtendedGcd(a, m)
  if g != 1:
    raise Exception('modular inverse does not exist')
  return x % m

[Python] ExtendedGcd code snippet

ExtendedGcd code snippet:

def ExtendedGcd(a, b):
  """
  Extended Euclidean algorithm.

  Computes, in addition to the greatest common divisor of integers a and b,
  also the coefficients of Bezout's identity, which are integers x and y such
  that a x + b y = gcd(a, b).

  Args:
    a: An integer.
    b: An integer.

  Returns:
    A tuple (g, x, y) satisfying a x + b y = g.
  """
  if a == 0:
    return (b, 0, 1)
  else:
    g, y, x = ExtendedGcd(b % a, a)
    return (g, x - (b // a) * y, y)

[Python] Factorize code snippet

Factorize code snippet: (depends on GetPrimes())

def Factorize(n, primes):
  """
  Factorize n by primes.

  For every prime p <= sqrt(n), we examine whether n is divided by p. If so, we
  compute the maximal integer value e satisfying that n is divided by p^e and n
  is not divide by p^{e + 1}.

  The factorization is represented by a dict(). It means that all prime factors
  of n are not in order. Consider use collections.OrderedDict() instead of
  dict() to keep the order of prime factors.

  Args:
    n: An integer representing a number to be factorized.
    primes: An integer list representing all possible prime factors <= sqrt(n).

  Returns:
    A dictionary where the key representing a prime factor p and the
    corresponding value representing the maximal integer value e satisfying that
    n is divided by p^e and n is not divide by p^{e + 1}.
  """
  factorization = {}
  d = n
  for prime in primes:
    if prime * prime > d:
      break
    e = 0
    while d % prime == 0:
      d //= prime
      e += 1
    if e > 0:
      factorization[prime] = e
  if d > 1:
    factorization[d] = 1

  if primes[-1] * primes[-1] < d:
    raise ValueError
  return factorization

[Python] GetPrimes code snippet

GetPrimes code snippet:

def GetPrimes(n):
  """
  Get all primes up to n.

  Use sieve of Eratosthenes to generate all primes less or equal to n. For
  example, given n = 10, then returns a list of primes [2, 3, 5, 7].

  Args:
    n: An integer representing the upper bound.

  Returns:
    A list containing all primes <= n in the increasing order.
  """
  sqrt_n = int(math.sqrt(n))
  visited = [False] * (n + 1)
  for i in range(2, sqrt_n + 1):
    if not visited[i]:
      for j in range(i * i, n + 1, i):
        visited[j] = True
  primes = []
  for i in range(2, n + 1):
    if not visited[i]:
      primes.append(i)
  return primes

[Python] ModExponentiation code snippet

ModExponentiation code snippet:

def ModExponentiation(b, e, m):
  """
  Modular exponentiation algorithm of time complexity = O(log(e)).

  Args:
    b: An integer representing the base.
    e: An integer representing the exponent.
    m: An integer representing modulus.

  Returns:
    An integer c = b^e (mod m).
  """
  if e == 0:
    return 1
  y = b
  output = 1
  while e > 0:
    if e & 1 == 1:
      output = (output * y) % m
    y = (y * y) % m
    e >>= 1
  return output

[Note] Totient summatory function in sublinear time

Original post: https://math.stackexchange.com/questions/316376/how-to-calculate-these-totient-summation-sums-efficiently

Background:

• Totient summatory function: Wolfram

Python code snippet:

class TotientSummatoryFunction():
  """Totient summatory function class."""

  def __init__(self, bound, mod):
    """Initialize.

    Args:
      bound: An integer.
      mod: An integer.
    """
    self.mod = mod
    self.small_values = None
    self.small_value_bound = None
    self.big_values = {}
    self.InitSmallValues(bound)

  def InitSmallValues(self, bound):
    """Initialize small values for totient summatory function.

    Args:
      bound: An integer.
    """
    n = 3 * int(math.pow(bound, 2 / 3))
    assert(n**3 >= bound**2)
    print('InitSmallValues({n})'.format(n=n))

    phi = [i for i in range(n + 1)]
    for i in range(2, n + 1):
      if phi[i] != i:
        continue
      for j in range(i, n + 1, i):
        phi[j] = phi[j] // i * (i - 1)

    self.small_values = [0 for i in range(n + 1)]
    self.small_values[1] = phi[1] % self.mod
    for i in range(2, n + 1):
      self.small_values[i] = (self.small_values[i - 1] + phi[i]) % self.mod
    self.small_value_bound = n

  def Get(self, n):
    """Get totient summatory function value.

    Args:
      n: An integer.

    Returns:
      An integer.
    """
    if n <= self.small_value_bound:
      return self.small_values[n]

    if n not in self.big_values:
      result = n * (n + 1) // 2
      m = self.IntegerSqrt(n)
      for d in range(1, m + 1):
        result -= self.Get(d) * (n // d - n // (d + 1))
      for k in range(2, n // (m + 1) + 1):
        result -= self.Get(n // k)
      self.big_values[n] = result % self.mod
    return self.big_values[n]

  def IntegerSqrt(self, n):
    """
    Get the integer part of Sqrt[n].

    Args:
      n: An integer.

    Returns:
      An integer.
    """
    guess = (n >> n.bit_length() // 2) + 1
    result = (guess + n // guess) // 2
    while abs(result - guess) > 1:
      guess = result
      result = (guess + n // guess) // 2
    while result * result > n:
      result -= 1
    return result