reduce(
  lambda x, y: x and y,
  map(
    lambda n:
        Matrix(n, lambda i,j: lcm(i+1,j+1)).det()
          ==
        factorial(n)*reduce(
          lambda x,y: x*y,
          map(
            lambda p: (1-p)^floor(n/p),
            filter(is_prime, range(n+1))
          ),
          1
        ),
    range(1,101)
  ),
  True
)