https projecteuler netproblem234 for an integer n 4 we define the lower prime square root of n denoted by lpsn as the largest prime n and the upper prime square root of n upsn as the smallest prime n so for example lps4 2 ups4 lps1000 31 ups1000 37 let us call an integer n 4 semidivisible if one of lpsn and upsn divides n but not both the sum of the semidivisible numbers not exceeding 15 is 30 the numbers are 8 10 and 12 15 is not semidivisible because it is a multiple of both lps15 3 and ups15 5 as a further example the sum of the 92 semidivisible numbers up to 1000 is 34825 what is the sum of all semidivisible numbers not exceeding 999966663333 sieve of erotosthenes function to return all the prime numbers up to a certain number https en wikipedia orgwikisieveoferatosthenes primesieve3 2 primesieve50 2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 computes the solution to the problem up to the specified limit solution1000 34825 solution10000 1134942 solution100000 36393008 get numbers divisible by lpscurrent reset the upperbound add the numbers divisible by upscurrent remove the numbers divisible by both ups and lps increment the current number remove twice since it was added by both ups and lps increment the current number setup for next pair sieve of erotosthenes function to return all the prime numbers up to a certain number https en wikipedia org wiki sieve_of_eratosthenes prime_sieve 3 2 prime_sieve 50 2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 computes the solution to the problem up to the specified limit solution 1000 34825 solution 10_000 1134942 solution 100_000 36393008 get numbers divisible by lps current reset the upper_bound add the numbers divisible by ups current remove the numbers divisible by both ups and lps increment the current number remove twice since it was added by both ups and lps increment the current number setup for next pair
import math


def prime_sieve(n: int) -> list:
    
    is_prime = [True] * n
    is_prime[0] = False
    is_prime[1] = False
    is_prime[2] = True

    for i in range(3, int(n**0.5 + 1), 2):
        index = i * 2
        while index < n:
            is_prime[index] = False
            index = index + i

    primes = [2]

    for i in range(3, n, 2):
        if is_prime[i]:
            primes.append(i)

    return primes


def solution(limit: int = 999_966_663_333) -> int:
    
    primes_upper_bound = math.floor(math.sqrt(limit)) + 100
    primes = prime_sieve(primes_upper_bound)

    matches_sum = 0
    prime_index = 0
    last_prime = primes[prime_index]

    while (last_prime**2) <= limit:
        next_prime = primes[prime_index + 1]

        lower_bound = last_prime**2
        upper_bound = next_prime**2

        
        current = lower_bound + last_prime
        while upper_bound > current <= limit:
            matches_sum += current
            current += last_prime

        
        while (upper_bound - next_prime) > limit:
            upper_bound -= next_prime

        
        current = upper_bound - next_prime
        while current > lower_bound:
            matches_sum += current
            current -= next_prime

        
        current = 0
        while upper_bound > current <= limit:
            if current <= lower_bound:
                
                current += last_prime * next_prime
                continue

            if current > limit:
                break

            
            matches_sum -= current * 2

            
            current += last_prime * next_prime

        
        last_prime = next_prime
        prime_index += 1

    return matches_sum


if __name__ == "__main__":
    print(solution())