# 
# Variants of the sieve of Eratosthenes (Python with NumPy)
# by Project Nayuki, 2017. Public domain.
# https://www.nayuki.io/page/the-versatile-sieve-of-eratosthenes
# 

import numpy


# Given an integer limit, this returns a list of Booleans
# where result[k] indicates whether k is a prime number.
def sieve_primeness(limit):
	if limit < 0:
		raise ValueError("Limit must be non-negative")
	result = numpy.ones([limit + 1], dtype=numpy.bool_)
	result[0] = False
	if limit > 0:
		result[1] = False
	for i in range(2, _sqrt(limit) + 1):
		if result[i]:
			result[i * i : : i] = False
	return result


# Given an integer limit, this returns a list of integers
# where result[k] is the smallest prime factor of k.
def sieve_smallest_prime_factor(limit):
	if not (0 <= limit < 2**32 - 1):
		raise ValueError("Limit out of range")
	result = numpy.zeros([limit + 1], dtype=numpy.uint32)
	if limit > 0:
		result[1] = 1
	for i in range(2, _sqrt(limit) + 1):
		if result[i] == 0:
			result[i] = i
			temp = result[i * i : : i]
			temp[temp == 0] = i
	mask = result == 0
	result[mask] = numpy.arange(result.size, dtype=numpy.uint32)[mask]
	return result


# Given an integer limit, this returns a list of integers
# where result[k] is the totient (Euler phi function) of k.
def sieve_totient(limit):
	if not (0 <= limit < 2**32 - 1):
		raise ValueError("Limit out of range")
	result = numpy.arange(limit + 1, dtype=numpy.uint32)
	for i in range(2, len(result)):
		if result[i] == i:
			result[i : : i] -= result[i : : i] // i
	return result


# Given an integer limit, this returns a list of integers where result[k]
# is the number of unique prime factors (omega function) of k.
def sieve_omega(limit):
	if not (0 <= limit < 2**32 - 1):
		raise ValueError("Limit out of range")
	result = numpy.zeros([limit + 1], dtype=numpy.uint32)
	for i in range(2, len(result)):
		if result[i] == 0:
			result[i : : i] += 1
	return result


# Given an integer limit, this returns a list of integers where result[k]
# is the product of the unique prime factors (radical function) of k.
def sieve_radical(limit):
	if not (0 <= limit < 2**32 - 1):
		raise ValueError("Limit out of range")
	result = numpy.ones([limit + 1], dtype=numpy.uint32)
	result[0] = 0
	for i in range(2, len(result)):
		if result[i] == 1:
			result[i : : i] *= i
	return result


# (Private) Returns floor(sqrt(x)) for an integer x >= 0.
def _sqrt(x):
	y = 0
	i = 1 << (x.bit_length() // 2)
	while i > 0:
		y |= i
		if y * y > x:
			y ^= i
		i >>= 1
	return y