Computes the Möbius transform of arithmetic sequences

mobius.py

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()