Abstract:

Introduction.  Square-free numbers:

   The first few square-free numbers give the sequence:
1, 2, 3, 5, 6, 7, 10, 11, 13, 14, 15, 17, 19, 21, 22, 23, 26, 29, 30, 31, 33, 34, 35, 37, 38, 39, 41, 42, 43, 46, 47, 51, 53, 55, 57, 58, 59, 61, 62, 65, 66, 67, 69, 70, 71, 73, 74, 77, 78, 79, 82, 83, 85, 86, 87, 89, 91, 93, 94, 95, 97, etc. (Sloane's A005117)

   More details on square-free numbers can be found at http://mathworld.wolfram.com/Squarefree.html.
(If a square-free number is used as the argument of the Möbius function, a non-zero value (+1 or -1) is obtained.)

Gaps between square-free numbers:

For any L, it can be shown that there exist infinitely many gaps of length greater than L in the sequence of square-free numbers.  Longer lists of square-free gaps and recent results are available here.  These gaps are series of consecutive squareful numbers (Sloane's A013929).
  Note that the term "squarefull" sometimes denotes a positive integer n such that if p is a prime dividing n, then p² divides n. (See the preprint of Michael Filaseta at http://www.math.sc.edu/~filaseta/paperindex.html, The distribution of squarefull numbers in short intervals, with O. Trifonov, Acta Arith., 67 (1994), 323-333.)

First occurrence of a gap of length L:

   Qgap(L) is listed in the following table for L up to 17.  The third column gives the prime factors that are repeated for each number in the gap.  The values of Qgap(L < 10) and Qgap(11) have been confirmed by different sources (see for example, Sloane's On-Line Encyclopedia of Integer Sequences, sequence A045882).  The sequence Qgap(L) is listed under A051681 in the Encyclopedia of Integer sequences.
 

Upper limits for Qgap(16) to Qgap(24)

* The value of Qgap(16) was found to be 3 215 226 335 143 218 (Z. McGregor-Dorsey, L. Marmet, J. Wetherell, et al., July 22^nd, 2000)
⁺ The value of Qgap(17) was confirmed to be 23 742 453 640 900 972 (E. Wong, Z. McGregor-Dorsey, L. Marmet, et al., July 8^th, 2001)

Basic algorithm: the sieve of Eratosthenes

   A faster algorithm is used by "Mathematica" to determine if a number is square-free.  The method is quite interesting (see:  http://documents.wolfram.com/v5/Add-onsLinks/StandardPackages/NumberTheory/NumberTheoryFunctions.html )

   However, to determine which of many consecutive numbers are square-free, an algorithm based on to the sieve of Eratosthenes is much faster.  The sieve of Eratosthenes is explained on the web site at http://mathworld.wolfram.com/SieveofEratosthenes.html.  It uses a list of numbers from which each composite number is removed.  Once the process is finished, only the prime numbers are in the list.

   To find the square-free numbers using a sieve, a similar technique is used but the algorithm eliminates numbers that are not square-free.  Starting with a list of integers, first cross out the multiples of 4:
 

1

2

3

X

5

6

7

X

9

10

11

X

13

14

15

X

17

18

19

X

21

22

23

X

25

26

...

then the multiples of 9, 25, etc., up to the last number in the list:
 

1

2

3

X

5

6

7

X

X

10

11

X

13

14

15

X

17

X

19

X

21

22

23

X

X

26

...

The remaining numbers are square-free numbers; the gaps are indicated by the series of consecutive "X".

Improvements of the algorithm:

   To implement this algorithm on a computer, it is not necessary to keep the entire list of integers in memory.  An improvement of the algorithm uses instead two shorter arrays to calculate the next non-square-free number after N:

• the first array, called p2, gives the squares of the prime numbers up to the largest number we want to test N_max,
• the second array, called nsqf, gives for each p2[i] the next non-square-free number, that is, the smallest number larger than N that is a multiple of p2[i].  This array can easily be calculated using modulo arithmetic.

   To find the square-free gaps, one finds sequences of non-square-free numbers.  The following example shows the arrays used to find gaps starting from N = 20.
(We use the notation used in the programming language C, where the index of an array starts at 0.)
 

   Using this table, it is easy to find the next non-square-free number: it is the smallest number in the array nsqf, that is, 24.  We set N = 24 and recalculate the array.  This is easy since the only needed operation is to add the corresponding prime squared to the multiple: 24+4 = 28.  The following table is obtained:
 

   Again, the next non-square-free number is the smallest nsqf[i].  By repeating this procedure, N takes the values of all the non-square-free numbers.
   It is advantageous to sort nsqf in increasing order at each step.  This way, the smallest is always nsqf[0].  We set N = 25 and recalculate the array to obtain:
 

   The order of the array primesq has also been changed so that each number p2[i] always corresponds to its multiple nsqf[i].  Repeating the procedure will generate the non-square-free numbers N = 27, 28, 32, 36, 40, etc.  Note that special care has to be taken when some numbers in nsqf are equal - each of these has to be increased by the value of its corresponding p2[i].

   The sort is relatively efficient since after nsqf[0] is given its new value, nsqf[1], nsqf[2] and the following elements are still in increasing order.  The new value has to be moved up the array until its proper place is found.

   If gaps of a given length L_min or more are searched, some nsqf[i] can be skipped.  To show this, we first find the minimum number of squared-primes NP2_min required in a gap of length L:
 

For example, if we chose L = 7, we see that NP2_min[L] = 6 prime factors are required for Qgap(7) starting at N = 217070 (in this case, the prime factors are 2, 3, 5, 7, 11 and 113).

   Continuing with the example given above, we now specifically search gaps with length L_min= 7 or longer.  There is no gap of length L_min= 7 in the interval starting at nsqf[0] and ending at nsqf[5], if nsqf[5] > nsqf[0] + 7. (In general, there is no gap of length L_min in the interval starting at nsqf[0] and ending at nsqf[NP2_min[L_min-1]], if nsqf[NP2_min[L_min-1]] > nsqf[0] + L_min.)

   With N = 40, we have:
 

   Since nsqf[5] = 169 > nsqf[0] + 7 = 51, there is no gap of length 7 in the interval starting at 44 and ending at 169.  We can therefore safely set N = nsqf[5] - L_min +1 = 163 and recalculate the following table:
 

This cuts down on the number of non-square-free that have to be tested and the speed of the calculation is increased.  A factor five in speed was obtained when this was implemented in the program which was used to find Qgap(14) and Qgap(15).

   This variation on Improvement II turns out to be more efficient, especially when it is combined with Improvements IV and V.  The trick is to consider the smallest squared-primes separately from the large squared-primes.
   Most of the time in the algorithm is spent on the process of taking the smallest elements off of nsqf and sorting them back into the array.  If one can reduce the number of non-square-free numbers which are tested, the speed of the algorithm will improve.
   This is actually possible since the smallest squared-primes are not needed to calculate the sieve.  If we have a gap of, say, L = L_min = 7 non-square-free numbers, then at least NP2_min[7] = 6 different primes are found in the gap.  This means that the smallest k₁ = NP2_min[L_min] -1 = 5 squared-primes can be left out of the calculation.  If we search for the 6^th squared-prime, every gap of length L = 7 (or more) will be found.
   We therefore separate the small squared-primes from the large ones, creating a base with k₁ elements.  The table for N = 163 would now look like this:
 

with the base shown in gray.  From the table, one sees that one must test for a gap having L_min= 7 around N = 169.  If none is found, then nsqf[k₁] is increased by 169 and sorted back into the array of large squared-primes to get:
 

   The sort is faster since it is only done on the large squared-primes.  The process is then continued at N = 289.

   To test if there is a gap of L_min around N, one must still know about multiples of the k₁ small squared-primes near N.  This can be done by trial division; even if trial division is slow, it is faster than resorting the base array.  One can also optimize the trial divisions, because the trial divisions for, say, N + 1 and N + 2 are related to each other.  For each small prime p, compute N % p² and store it in a list mod (the "%" symbol is the modulo function in the C language).  Now to see if p² divides N + 1, we just test if this stored value is -1 (mod p²).  To see if p² divides N + 2, we test if this stored value is -2 (mod p²).  (We also precompute the value of -1 (mod p²), -2 (mod p²),  etc.,  for the small set of values which we will possibly need to test.)  Note that it is also necessary to test N - 1, N - 2, etc.  For N = 289, we have the following arrays:
 

   We see that p²[0] and p²[1] divide N - 1 = 288, but this is the only other non-square-free number for this gap of length 2.  We can therefore increase nsqf[k₁] by 289, sort it back into the array (reordering p2 accordingly) and continue with N = 338:
 

   If we include p2[4] in the array of large squared-primes, (so there are only k₂ = NP2_min[L_min] -2 = 4 elements in the base), then we know that a number N = nsqf[k₂] we are testing can be part of a gap only if the next number in nsqf is close to N, that is, if nsqf[k₂+1] is not larger than N + L_min-1.  With N = 338, the arrays are now:
 

   Since nsqf[k₂+1] = 361 is larger than N + L_min-1 = 344,  we can skip to N = 361.  Since N does not pass the closeness test, we do not have to do any computations with the base, saving us a lot of time.

   A chained list can be built with a set of numbers that specify an order for the elements of an array.  If we use a chained list represented by the array called next such that nsqf[next[i]] >= nsqf[i], we can sort the array nsqf without moving any data within tha arrays nsqf or p2.  For a reason that will become obvious later, we choose to have p2 sorted in increasing order.  With N = 361, the arrays would be:
 

   We use an additional variable, "head", which points to the smallest item in the array nsqf.  For the table above, we have head = 7 (next[head] is highlighted).  Since we only need to change two values in the array next to perform a sort, this method is faster.  Searching through the array nsqf now consists of going through the data in the following order:

 for (i=head; tempnsqf>=nsqf[next[i]]; i=next[i]) { ... }

   There is another advantage to this method: since the array p2 is sorted in increasing order, we know that if we have to sort item nsqf[i = head], then nsqf[i-1] <= nsqf[i] + p2[i].  This means we can jump immediately to i = head-1 and start searching from there.  In general, this cuts the search in half! Searching through the list now consists of:

  for (i=head-1; tempnsqf>=nsqf[i2=next[i]]; i=i2) { ... }

   We continue the example with the above table.  Since nsqf[next[head]] = 363 is not larger than N + L_min-1 = 367, we calculated the modulos.  Clearly, there is only a gap of length L = 2 starting at 360.  We therefore set N = 363 and sort tempnsqf = nsqf[head] =nsqf[head] + p2[head] = 722.  We start with i = 6 and the sort requires only one comparison!  The following arrays are obtained:
 

with head now equal to 4.

   A factor three in speed was obtained when this algorithm was implemented in the program.

   Instead of having k₂ elements in the base and look for two large squared-primes that are not too far apart, we can use a base with  k₃ =NP2_min[L_min] -3 = 3 squared-primes, but look to see if nsqf[head] and nsqf[next[next[head]]] are not more than L_min apart.
 

In this example, we look at the difference between nsqf[head] = 363 and nsqf[next[next[head]]] = 507.  Since the difference between the two values is larger than L_min-1, there is no gap of length L_min or longer.

   This improvement gives the program a 37% speed increase with L_min = 14.

   The program becomes slower if four or more large squared-primes are considered (I tested this for L>13, N = 10¹⁴  to 10¹⁴ +10⁹).

Evaluation of order of algorithm:

  First, consider the case when p2[0]=4 which occurs with a probability of 1/4.  Since nsqf[0] has been increased by 4, there is a probability of 4/9 that it will have to be moved above nsqf[j] (if p2[j]=9).  There is an additional probability of 4/25 that nsqf[0] will have to be moved above nsqf[k] (if p2[k]=25), etc.  We therefore get the average number of steps required to place the new nsqf[0] to its correct position:

S(1) = 1/4 × (4/9 + 4/25 + 4/49 + 4/121 + ...) = S_i=2,... 1/(p²(i))

where p(i) is the i^th prime number  (p(1)=2, p(2)=3, p(3)=5, etc.) and p²(i) = p(i)×p(i).

  In the case when p2[0]=9 (which occurs with a probability of 1/9), one move is always necessary to bring nsqf[0] above nsqf[j] (if p2[j]=4).  There is a probability of 9/25 that nsqf[0] will have to be moved above nsqf[k] (if p2[k]=25), etc.  We therefore get:

S(2) = 1/9 × (1 + 9/25 + 9/49 + 9/121 + ...) = 1/(p²(2)) + S_i=3,... 1/(p²(i))

  In general, when p2[0]=p²(m), the average number of moves is:

S(m) = (m-1)/(p²(m)) + S_i=m+1,... 1/(p²(i))

  The sum over all the S(m) gives the average number of moves required to place nsqf[0] to its correct position in the list:

S_m=1,... S(m) = 2 × S_i=2,... (i-1)/(p²(i))

This series converges, as determined with the convergence test 0.224 given in Table of Integrals, Series, and Products by Gradshteyn and Ryzhik (Academic Press, Inc., p. 5).  It converges very slowly to approximately 1.30...  Therefore, given the remainders for N-1, the number of operations required to find out if N is square-free is independent of the value of N.

Acknowledgements:

   Contributed to find Qgap(16): Zach McGregor-Dorsey, Louis Marmet, Joe Wetherell, Gunnard Engebreth, D. Bernier, Erick Wong, Alan Simpson and Nicolas Marmet.

   Contributed to find Qgap(17): E. Wong, Z. McGregor-Dorsey, L. Marmet, Jean-Pierre Bernier, D. Bernier, Nancy Robertson, N. Marmet, Charles Ward and G. Engebreth.

   The search for Qgap(18) is an ongoing team project.  Contributions were made by D. Bernier, L. Marmet, E. Wong, J. Wetherell, Z. McGregor-Dorsey, G. Engebreth, A. Simpson, N. Marmet, N. Robertson, J.-P. Bernier, C.R. Ward, Bruno Le Tual and Horand Gassmann.

   This project started from an idea that was initially suggested to me by David Bernier.

Related web pages:

The Prime Puzzles and Problems Connection.  Carlos Rivera. http://www.primepuzzles.net/problems/prob_028.htm (not updated since Nov. 1999)

B. de Weger and C.E. van de Woestijne, http://www.xs4all.nl/~deweger/publikaties.html
   On the powerfree parts of consecutive integers, Acta Arithmetica 90 [1999], 387 - 395.
 

Copyright © 1999  Louis Marmet