Iterating over all squarefree composites $\lt L$

Question

I need to iterate over all integers that are:

1. Composite
2. Squarefree
3. Smaller than a limit $L$

And see for each one if it has a certain property (in addition to these 3). Then I need to sum all such integers with this property. The following code does exactly that:

    public static void main(String[] args) 
    {
        long L=(long) Math.pow(10,12);
        long [] primes = eratosthenes(L/2);
        long sigma=0;
        for (int i=0;i<primes.length-1;i++)
            sigma+=squareFreeCompositesGenerator(L,primes[i],i+1,primes);
        System.out.println(sigma);
    }

    public static long squareFreeCompositesGenerator(long L, long prod, int position, long [] primes)
    {
        long sum=0;
        for (int i=position;i<primes.length&&prod*primes[i]<=L;i++)
        {
            if (property(prod*primes[i]))
                sum+=prod*primes[i];
            sum+=squareFreeCompositesGenerator(limit,prod*primes[i],i+1,primes);
        }
        return sum;
    }

The idea behind this, is that rather than checking all numbers up to $L$ if they're squarefree, I can rather simply create/generate them, which is considerably faster compared to iterating through all numbers up to $L$ and factor and check each one. In this code, a call for eratosthenes(n) returns a sieve of all primes up to $n$. The method squareFreeCompositesGenerator then "creates" all squarefree composites smaller than $L$, and checks if they have the property through the property method. Through recursion, it keeps multiplying the primes without repetition and runs through all possible combinations. The "Achilles heel" of this algorithm though, is in the very beginning, where I require all primes $p \leqslant \frac L2$. This is due to the fact that the largest possible such integer can be of the form $2 \cdot p$ where $p$ is the largest prime up to $ \frac L2$. The program runs swiftly even for $L=10^9$ but for $L$ as big as $10^{12}$, I can't create and store all those primes in reasonable time, and even if I could store them, it would still take a long time.

It is important to note I don't need to sum all squarefree composites $\lt L$, or count how many there are, but rather iterate through all of them and check if they have a certain property.

Is there any way to adjust the code so that it only requires primes up to $\sqrt n$? Or is there perhaps another approach to iterate through all squarefree composites $\lt L$ with a sublinear time complexity?

Answer 1

Try this

 N = 20; mu = [1,-ones(1,N-1)]; primes = [0,ones(1,N-1)]; 
 for n = 2:N, 
   mu(n+n:n:end) = -mu(n) + mu(n+n:n:end); primes(n+n:n:end) = 0; 
 end; 
 squarefree = abs(mu); find(primes==1), find(squarefree==1)

on https://octave-online.net/

Answer 2

I believe you can relatively efficiently do what you're requesting by using that the set of positive integers which are composite and square-free is the same as the set of all integers $\ge 2$ except those which are prime or multiples of prime squares.

To relatively efficiently handle this exclusion, first use a type of Wheel factorization to determine, for some relatively small $n$'th prime $p_n$, the squares of the primes $p_1 = 2$ to $p_n$ multiples up to $(p_n\#)^2$, where $p_n\#$ is the Primorial, i.e., product of all primes up to $p_n$. For example, for $2$ and $3$, going up to $(2 \times 3)^2 = 36$, the multiples of $4$ are $0, 4, 8, 12, 16, 20, 24, 28, 32$ and the multiples of $9$ are $0, 9, 18, 27$. Thus, any integer which has a value modulo $36$ of any of the remaining values of $1, 2, 3, 5, 6, \ldots, 31, 33, 34, 35$ will not be a multiple of $4$ or $9$. You can extend this to handle several more primes. I suggest using something like $(11\#)^2 = (2,310)^2 = 5,336,100$. Also, I recommend you just store the offsets (only one byte per value is needed since they are all quite small), with the last one wrapping around to the first, so you can just simply add the next offset to your current integer being checked.

Next, for the remaining primes up to the square root of $L$, you can store the squares of these primes in $2$ arrays, one for the square value itself (e.g., for updating the second array) and one for the next multiple of the square value to check. Also, keep track of the next prime square larger than the current value you're dealing with. This way, you can limit the iterations in your check loop, plus can just do comparisons, adding and/or updating values, with these operations all being relatively efficient.

With this in place, you can then just use an outer loop to get each next integer which doesn't have a small prime square factor. First, check if this integer is larger than or equal to the next prime (by having some sort of dynamic feed of prime numbers up to your limit $L$ of the next prime available to check against). If it's not a prime, then check if it's a multiple of a prime square (by checking & using those $2$ arrays I mentioned in the previous paragraph). If both of these checks pass, then if your certain property also passes, add the integer to your cumulative sum.

Factoring, which is generally quite computer intensive, is avoided. Apart from possibly the time it takes to check your certain property, one of the most computer intensive operations remaining is finding the primes. Fortunately, there are various techniques to do this quite efficiently. For example, there's Python bindings of a fast C++ primesieve library (with source) at Github. Also, I'm currently doing some primes research where I'm using a modified version of a relatively efficient public domain C++ implementation that primesieve provides of a segmented sieve of Eratosthenes. Note that, as stated there, the algorithm uses $O(n \log \log n)$ operations and $O(\sqrt{n})$ space. It seems your sample code is written in C#, so if you wish, I assume you can either use this in C++ or convert it to C#. On my computer (AMD Ryzen $5$ $1400$ Quad-Core Processor, $3.20$ GHz, $16$ GB RAM, $384$ KB L$1$ cache, $2$ MB L$2$ cache, $8$ MB L$3$ cache & running $64$-bit Windows $10$), it took my somewhat more efficient C++ version $0.672$ hours (i.e., $40.32$ minutes) to generate the primes up to $10^{12}$, and $15.342$ hours to generate the primes up to $2.1 \times 10^{13}$.

Answer 3

There's a discussion here about primes, including python code for popping out squarefree integers. There is a prime generator which serves as the basis for the square factor sieve included.

https://stackoverflow.com/a/33900585/4363639