Skip to content

Instantly share code, notes, and snippets.

@carlosgmartin
Created July 13, 2019 00:08
Show Gist options
  • Save carlosgmartin/d9b9cb7fcd4263c437104f67a9522767 to your computer and use it in GitHub Desktop.
Save carlosgmartin/d9b9cb7fcd4263c437104f67a9522767 to your computer and use it in GitHub Desktop.
Computes the Möbius transform of arithmetic sequences
import numpy as np
import sympy.ntheory
import time
# https://mathoverflow.net/a/227408/74578
def mobius_transform(sequence):
sequence = sequence.copy()
for i in range(1, len(sequence)//2+1):
sequence[i+i-1::i] -= sequence[i-1]
return sequence
n = 10**4
unit = np.zeros(n, dtype=int)
unit[0] = 1
start = time.time()
mobius_transform(unit)
print('{:.3f} seconds to compute Möbius sequence up to {}'.format(time.time() - start, n))
zeros = np.zeros(n, dtype=int)
ones = np.ones(n, dtype=int)
identity = np.arange(n) + 1
totient = np.array([sympy.ntheory.totient(i+1) for i in range(n)])
mobius = np.array([sympy.ntheory.mobius(i+1) for i in range(n)])
divisor_count = np.array(np.array([sympy.ntheory.divisor_count(i+1) for i in range(n)]))
is_square = divisor_count % 2
primeomega = np.array([sympy.ntheory.primeomega(i+1) for i in range(n)])
liouville = (-1) ** primeomega
assert (mobius_transform(zeros) == zeros).all()
assert (mobius_transform(unit) == mobius).all()
assert (mobius_transform(ones) == unit).all()
assert (mobius_transform(identity) == totient).all()
assert (mobius_transform(divisor_count) == ones).all()
assert (mobius_transform(is_square) == liouville).all()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment