Fast method for Nth Squarefree number (using mathematica)

Question

I am trying to compute Nth Squarefree numbers using Mathematica. What I am trying to utilize is the SquareFreeQ[i] function.

Here is my solution :

NSqFree[N_] := For[n = 1; i = 1, n <= N + 1, If[SquareFreeQ[i], If[n == N, Print[i] ] n++] i++]

But I am supposed to compute NSqFree[1000000000] but seems like my approach is taking for ever. Any faster method ?

Thanks,

ADDED:

Here an exactly identical topcoder question and the corresponding editorial for the same.

Answer 1

You have to us the Inclusion-Exclusion principle: suppose you want to find the number of square free numbers up to $N$, then from $N$ you have to substract the number of integers divisible by the square of a prime, but then you have to add any multiple of the square of the product of two discinct primes and so on, in formulas the number looked for is $$ N - \sum_{p^2 \le N} \left\lfloor\frac{N}{p^2}\right\rfloor + \sum_{p^2q^2 \le N} \left\lfloor\frac{N}{p^2q^2}\right\rfloor - \sum_{p^2q^2r^2 \le N} \left\lfloor\frac{N}{p^2q^2r^2}\right\rfloor + \cdots -\cdots $$ using the moebius $\mu$ function you can write this last formula as $$ \sum_{n \le \sqrt{N}} \mu(n) \left\lfloor\frac{N}{n^2}\right\rfloor $$ I don't know how to write this in mathematica but it should take a negligible fraction of the time it takes your current method.

Answer 2

You won't do a billion (if I counted the zeros right) loops any time soon, so you need a better approach. Wikipedia has an approximation under the section Distribution of square-free numbers. So you could start with $10^9 \pi^2/6$ and calculate how many below that are square-free, then correct from there. You will need to use inclusion/exclusion, but there are only about 40,000 to consider.

Answer 3

I think that your method can be improved by sieving: if $n$ is not squarefree you don't need to check any of its multiples. Start with a reasonable estimate (like the one Ross Millikan gave) and sieve like you would do in a prime sieve.

Is this viable? I don't know! I can't program in Mathematica, but later today I'll try to implement this in Python.

Answer 4

I'm not an expert of number-theoretic algorithms, but it seems to me that you can employing the Chinese Remainder Theorem to obtain a decent first stab at a "sieving" algorithm.

• Use several registers r[p], each storing a residue modulo p² for some prime p. Define such a register for various small primes p ∈ {2, 3, 5, ... , P} for some suitably large prime P. You will use these to represent an integer R, such that R ≡ r[p] (mod p²). You shouldn't need too many such registers to faithfully represent even reasonably large non-negative numbers (the registers can uniquely determine any integer from 1 to 2² × 3² × 5² × ... × P²).

• Whenever you wish to increment R, increment each of the registers r[p] as well. If R is square-free, none of these registers will be zero modulo the square of the appropriate prime. For sufficiently small integers N, you can even characterize N as being square-free if none of these residues are zero. Put another way, the more registers r[p] you maintain, the larger the range of square-free numbers you can automatically detect using these registers.

Suppose that you want to test for square-freeness up to some upper limit N. What do you need to test square-freeness using nothing but registers such as I've described? What you need is for any composite number less than N to have a prime factor from the list of registers that you maintain; that is, you need registers for each of the primes up to √N.

If you're going to iterate through all integers anyway, you can discover the list of primes that you need to characterize square-freeness at the same time by the Sieve of Erasthotenes, and construct the list of residues modulo squares-of-primes as you go. Any time you find a new prime P, so long as P² < N, add P to the list of primes for which you maintain registers and initialize a register for residues modulo P².

As Ross suggests in his answer, you can use results on the distribution of square-free numbers to obtain an upper bound for the N^th square-free number (it should grow more slowly than N ln(N) or so, in any case). This gives you an upper bound on the number of registers that you have to maintain as you search for square-free integers.

From this, you should actually be able to show a polynomial-time asymptotic upper bound on the runtime of your procedure.