CCalculations relating to somespecial Harmonic numbers
John Blythe Dobson ([email protected])August 18, 2018
Abstract
We report on the results of a computer search for primes p which dividean Harmonic number H b p/N c with small N > Keywords : Harmonic numbers, anharmonic primes, Fermat quotient, Fer-mat’s Last Theorem
Before its complete proof by Andrew Wiles, a major result for the first caseof Fermat’s last theorem (FLT), that is, the assertion of the impossibility of x p + y p = z p in integers x, y, z , none of which is divisible by a prime p > p must satisfy thecongruence q p (2) := 2 p − − p ≡ p ) , where q p (2) is known as the Fermat quotient of p , base 2. This celebrated resultwas to be generalized in many directions. One of these was the extension ofthe congruence to bases other than 2, the first such step being the proof of ananalogous theorem for the base 3 by Mirimanoff [23] in 1910. Another grewfrom the recognition that some of these criteria could be framed in terms ofcertain special Harmonic numbers of the form H b p/N c := b p/N c X j =1 j (1)for N >
1, and with b·c denoting the greatest-integer function. There are nicehistorical overviews of these developments in [28] and [27] (pp. 155–59), whichgo into greater detail than we attempt here.1 a r X i v : . [ m a t h . N T ] A p r Congruences for H b p/N c involving only Fermatquotients and low-order linear recurrent se-quences At the time when Wieferich’s and Mirimanoff’s results appeared, it was alreadyknown that three of these Harmonic numbers had close connections to the Fer-mat quotient, and satisfied the following congruences (all modulo p ): H b p/ c ≡ − · q p (2) (2) H b p/ c ≡ − · q p (3) (3) H b p/ c ≡ − · q p (2) . (4)All these results are due to Glaisher; those for N = 2 and 4 will be found in([12], pp. 21-22, 23), and that for N = 3 in ([11], p. 50). Although it was notthe next Harmonic number criterion published, it will be convenient to dispensenext with the case N = 6. The apparatus needed to evaluate this Harmonicnumber appears in a 1905 paper of Lerch ([21], p. 476, equations 14 and 15),but the implications of Lerch’s result were long overlooked, and only realizedin 1938 by Emma Lehmer ([20], pp. 356ff), who gave the following congruencemod p : H b p/ c ≡ − · q p (2) − · q p (3) . (5)The fact that the vanishing of H b p/ c mod p is a necessary condition for thefailure of the first case of FLT for the exponent p is an immediate consequenceof the theorems of Wieferich and Mirimanoff. In the present study, we useLehmer’s congruence in the equivalent form H b p/ c ≡ − · q p (432) , (6)obtained from (5) by applying in reverse the logarithmetic and factorizationrules for the Fermat quotient given by Eisenstein [8]. This expression revealsthat divisors p of H b p/ c are instances of the vanishing of the Fermat quotientmod p for composite bases, a problem which has notably been studied in theongoing work of Richard Fischer [10].The four congruences above exhaust the cases of H b p/N c that can be evalu-ated solely in terms of Fermat quotients. Historically, the case N = 5 also hasits origins in this era, and was in fact introduced before that of N = 6. In 1914,Vandiver [35] proved that the vanishing of both H b p/ c and q p (5) are necessaryconditions for p to be an exception to the first case of FLT. Unlike the four casesalready considered, here the connection of the Harmonic number with FLT wasdiscovered before any evaluation of it (beyond the definitional one) was known.It was almost eighty years later that the connection between these results would2ecome apparent. The ingredients needed for the evaluation of this Harmonicnumber were presented in a 1991 paper by Williams ([37], p. 440), and almostsimultaneously by Z. H. Sun ([31], pt. 3, Theorems 3.1 and 3.2); and thoughthey do not write out the formula explicitly, it is clearly implied to be H b p/ c ≡ − · q p (5) − · F p − ( p ) /p (mod p ) , (7)where F p − ( p ) /p is the Fibonacci quotient (OEIS A092330), with F a Fibonaccinumber and (cid:16) p (cid:17) a Jacobi symbol. In light of Vandiver’s theorems on H b p/ c and q p (5), this result immediately established that the vanishing of the Fibonacciquotient mod p was yet another criterion for the failure of the first case of FLTfor the exponent p . The fact was announced shortly afterwards in the celebratedpaper by the Sun brothers [32], which gave fresh impetus to an already vast lit-erature on the Fibonacci quotient; and in its honor the primes p which dividetheir Fibonacci quotient were named Wall-Sun-Sun-Primes. These remain hy-pothetical, as not a single instance has been found despite tests to high limits[24].The remaining known formulae for the type of special Harmonic numbersin which we are interested are of much more recent origin. In some cases theywere discovered simultaneously, or nearly so, by more than one researcher; andwe hope we have not done injustice to any of the participants. Apart from afew published formulae which are apparently in error or underdetermined, wehave the following congruences (all modulo p ) which are undoubtably correct: H b p/ c ≡ − · q p (2) − · U p − ( p ) (2 , − /p (8) H b p/ c ≡ − · q p (2) − · q p (5) − · F p − ( p ) /p (9) H b p/ c ≡ − · q p (2) − · q p (3) − · (cid:18) p (cid:19) · U p − ( p ) (4 , /p (10) H b p/ c ≡ − · q p (2) − · U p − ( p ) (2 , − /p − S − /p (11) H b p/ c ≡ − · q p (2) − · q p (3) − · U p − ( p ) (2 , − /p − · (cid:18) p (cid:19) · U p − ( p ) (4 , /p − · (cid:18) p (cid:19) U p − ( p ) (10 , /p (12)where the (cid:16) · p (cid:17) are Jacobi symbols, • F p − ( p ) /p is as before the Fibonacci quotient (OEIS A092330) • U p − ( p ) (2 , − /p is the Pell quotient (OEIS A000129)3 U p − ( p ) (4 , /p is a quotient derived from the Lucas sequence 1, 4, 15, 56,209, . . . (OEIS A001353) • U p − ( p ) (10 , /p is a quotient derived from the Lucas sequence 1, 10, 99,980, 9701, . . . (OEIS A004189) • (cid:16) p (cid:17) = ( − ( p − / • (cid:16) p (cid:17) = ( − b ( p +1) / c • (cid:16) p (cid:17) = ( − b ( p +5) / c .and S = ( − b p/ c + b p/ c × (cid:0) − C n − − C n + C n +2 (cid:1) if p ≡ ± (cid:0) C n − + C n − C n +1 (cid:1) if p ≡ ± (cid:0) C n − − C n +1 (cid:1) if p ≡ ± (cid:0) − C n − (cid:1) if p ≡ ± , (13)with C = C = C = C = 0 , C = 1 , C = 4 , C = − , C = 14; (14a) C n = 8 C n − − C n − + 16 C n − − C n − . (14b)The result for N = 8 is derived from the 1991 paper by Williams ([37], p.440), with an equivalent result also appearing in Sun ([31], pt. 3, Theorem 3.3).To the best of our knowledge, at the time of the discovery of this formula in1991 the vanishing of H b p/ c modulo p was not recognized as a condition for thefailure of the first case of FLT for the exponent p , and this only became evidentwith the appearance of the 1995 paper of Dilcher and Skula [4] discussed below.The result for N = 10 is due to a 1995 paper by Z. H. Sun ([31], pt. 3,Theorem 3.1). Since the vanishing of all the individual components was by thenknown to be a necessary condition for the failure of the first case of FLT for theexponent p , the same thing was immediately seen to be true for H b p/ c .The result for N = 12 is also derived from the 1991 paper of Williams ([37],p. 440), that for N = 16 is simplified from a 1993 paper by Z. H. Sun ([31], pt.2, Theorem 2.1), and that for N = 24 is from a 2011 paper of Kuzumaki &Urbanowicz ([19], p. 139).In their landmark joint paper of 1995, Dilcher and Skula [4] proved amongother things that the vanishing modulo p of H b p/N c was a necessary criterionfor the failure of the first case of FLT for the exponent p , for all N from 2 to 46.This result gives retrospective interest to the evaluations of H b p/ c and H b p/ c ,and furnishes the main motivation for the present study.All of the Fermat quotients and recurrent sequences mentioned above wereevaluated in PARI using modular arithmetic.4 The remaining cases of H b p/N c Although the calculations in the remaining cases, where no special formula for H b p/N c was known, were initially performed in the obvious way by summing ofmodular inverses over the required range, a better method was later found; thepresentation here for the most part follows that in [6]. The underlying idea wasdeveloped by the twin brothers Zhi-Hong Sun and Zhi-Wei Sun, who worked itout in two papers published just over a decade apart. Z.-H. Sun considered, asa special case of an even more general problem, the lacunary sum of binomialcoefficients T ( N, m ) := m X j =0 j ≡ N ) (cid:18) mj (cid:19) . This sum is well known in the literature in its own right, having been reducedto a series in N in an important paper of 1834 by C. Ramus [26], who showedthat T ( N, m ) = 1 N N − X j =0 (cid:18) jπN (cid:19) m cos jmπN ( m > . (15)This can also be written T ( N, m ) = 1 N N − X j =0 (cid:16) ω j (cid:17) m ( m > , (16)where ω = e πiN is a primitive N th root of unity. This formula was used byHoggatt and Alexanderson [14] to derive results equivalent to the congruences(2) through (8), and Howard and Witt [15] extended this to the case N = 10(9); in the cases N = 5 , ,
10 these results anticipate those cited above.In this section, we require some supplementary notations. H b p/N c is the case k = 0 of a sum studied by Lerch [21] and other writers, s ( k, N ) := b ( k +1) pN c X j = b kpN c +1 j = p j , (17)where it is always assumed that p is sufficiently large that s ( k, N ) contains atleast one element; the provision j = p is necessary when k + 1 = N . It is alsoconvenient to define an alternating version of this sum, s ∗ ( k, N ) := b ( k +1) pN c X j = b kpN c +1 j = p ( − j j . (18)5.-H. Sun ([31], pt. 1) showed that for a prime p , we have N (1 − T ( N, p )) p ≡ ( s (0 , N ) (mod p ) if N is even s ∗ (0 , N ) ≡ − s (1 , N ) (mod p ) if N is odd . (19)Later, Z.-W. Sun ([34]; see also [33]) completed the analysis by introducing theauxiliary expression T ∗ ( N, m ) := m X j =0 j ≡ N ) ( − j/N (cid:18) mj (cid:19) (20)= 2 · T (2 N, m ) − T ( N, m ) . (21)Because our interest is predominantly in the harmonic numbers rather than thesums of binomial coefficients, we motivate this device in a slightly different waythan in the original paper. Notice that by definition, s (0 , N ) ≡ s (0 , N ) + s (1 , N ) (mod p ) , (22)and by the complementarity of the elements of the residue system of a prime p , s ( k, N ) ≡ − s ( N − k − , N ) (mod p ) . Let s ( k, N ) denote the terms with odd denominators in s ( k, N ), and s ( k, N )the terms with even denominators. Then s ∗ (0 , N ) ≡ s (0 , N ) − s ( k, N ) (23) ≡ s (0 , N ) + s ( N − , N ) (24) ≡ s (0 , N ) + 12 s ( N − , N ) (25) ≡ − s (1 , N ) (mod p ) , (26)where the last line follows from Corollary 3.2 of Dilcher and Skula ([3], p. 20),which makes explicit a hint supplied in Lerch [21]. Thus from (22) we have foreven N , s ∗ (0 , N ) ≡ s (0 , N ) − s (0 , N ) (mod p ) , (27)while for odd N s (0 , N ) ≡ s (0 , N ) − s ∗ (0 , N ) (mod p ) . (28)Substituting these last two expressions into (19), we obtain its complement, N (1 − T ∗ ( N, p )) p ≡ ( s (0 , N ) (mod p ) if N is odd s ∗ (0 , N ) ≡ − s (1 , N ) (mod p ) if N is even . (29)6 .1 Recurrent sequence representations (I) Now the sequences associated with s (0 , N ), namely T ( N, m ) ( N even) and T ∗ ( N, m ) ( N odd), satisfy the same simple recurrence relations of order N − a n − (cid:18) N (cid:19) a n − + (cid:18) N (cid:19) a n − − . . . ∓ (cid:18) NN − (cid:19) a n − N +1 = 0 (30)according as N is even or odd; in other words, N − X j =0 ( − j (cid:18) Nj (cid:19) a n − j = 0 . (31)Similary, the sequences associated with s ∗ (0 , N ), namely T ( N, m ) ( N odd) and T ∗ ( N, m ) ( N even), satisfy recurrence relations also of order N − a n − (cid:18) N (cid:19) a n − + (cid:18) N (cid:19) a n − − . . . ∓ (cid:18) NN − (cid:19) a n − N +1 ± a n − N = 0 (32)according as N is even or odd; in other words, ± a n − N + N − X j =0 ( − j (cid:18) Nj (cid:19) a n − j = 0 . (33)These results are easily deduced from [18], p. 12. We will not attempt toreplicate here Z.-W. Sun’s proof that T ( N, m ) can be expressed in terms oflinear recurrent sequences of order not exceeding φ ( m ) /
2, where φ ( · ) is theEuler totient function. However, we will note the advantageous decomposition s (0 , N ) ≡ s (0 , N ) + s ∗ (0 , N ) (mod p ) , (34)where by symmetry the right-hand side has the same value regardless of theparity of N , a value in agreement with (19).A few examples of the characteristic polynomials in difference root form aregiven in Tables 1 and 2 below. Table 1: Recurrent sequences for T ( N, m ) for small values of N N coefficients of the a ’s in a n + a n − + a n − + . . . = 02 1 , −
23 1 , − , , −
24 1 , − , , −
45 1 , − , , − , , −
26 1 , − , , − , , −
67 1 , − , , − , , − , , −
28 1 , − , , − , , − , , −
89 1 , − , , − , , − , , − , , −
210 1 , − , , − , , − , , − , , − able 2: Recurrent sequences for T ∗ ( N, m ) for small values of N N coefficients of the a ’s in a n + a n − + a n − + . . . = 02 1 , − ,
23 1 , − ,
34 1 , − , , − ,
25 1 , − , , − ,
56 1 , − , , − , , − ,
27 1 , − , , − , , − ,
78 1 , − , , − , , − , , − ,
29 1 , − , , − , , − , , − ,
910 1 , − , , − , , − , , − , , − , The most efficient representations of the sums (19) and (29) for computationalpurposes is believed to be in the form of matrices. While the direct evaluationof (1) as a sum of reciprocals (or modular inverses) has an algorithmetical com-plexity of p log pN , that for evaluation by matrix exponentiation is only N log p ,and is clearly to be preferred except when N is large relative to p , and in thepresent study it never is.In what follows, whenever we speak of a term of a sequence being representedby a matrix (or the power of a matrix), we mean by the first term of the first rowof the matrix (or of the resulting power); and not being aware of any standardnotation for this relationship, we have used the symbol . = to express it.For the evaluation of recurrences of order N − N −
1. However, wehave preferred to use suitable matrices of dimension N of a very simple typewhose powers generate the terms of T ( N, m ) and T ∗ ( N, m ), feeling that theslight loss of speed with this approach was compensated for by the great ease ofimplementing it. The matrices for T ( · , m ) are circulant, while the matrices for T ∗ ( · , m ) are skew-circulant, and they differ only in the sign of the first term inthe last row. We note some small examples to illustrate the pattern: T (2 , m ) . = ! m T ∗ (2 , m ) . = − ! m T (3 , m ) . = m T ∗ (3 , m ) . = − m (4 , m ) . = m T ∗ (4 , m ) . = − m T (5 , m ) . = m T ∗ (5 , m ) . = − m . It is much faster to calculate the powers of these matrices than to calculate theassociated sequences directly, and the handling of matrix powers under modulararithmetic is very well implemented in PARI. The fact that the cost of evaluatingpowers of a matrix is in general proportional to the cube of its dimension d isanother reason for introducing (34) above, because the reduction of a calculationof order d to two calculations of order ( d ) involves a fourfold saving.This method has enabled the discovery of the following solutions: for N =9, p = 532199813; for N = 13, p = 427794751; for N = 17, p = 590422517. Several writers have observed that the recurrences for T ( N, m ) considered above(corresponding to s (0 , N ) for even N and s ∗ (0 , N ) for odd N ) can be substan-tially simplified by multisecting the sequences into separate congruence classesof m modulo N . Several OEIS entries illustrate the case m ≡ N ), eventhough they do not always explicitly state the corresponding recurrence (OEISA007613 for N = 3, A070775 for N = 4, A070782 for N = 5, A070967 for N = 6, A094211 for N = 7, A070832 for N = 8, A094213 for N = 9, A070833for N = 10). We in fact avoid the case m ≡ N ) because it is not neededfor testing prime values of m and increases the order of the recurrence by 1when N is even. Table 3 gives the values of the coefficients of some of thesesequences in difference root form. 9 able 3: Recurrent sequences for multisection of T ( N, m ) with m ≡ N ) for small values of N N coefficients of the a ’s in a n + a n − + a n − + . . . = 02 1 , −
43 1 , − , −
84 1 , − , −
645 1 , − , − ,
326 1 , − , − , , − , − , , , − , − , , , − , − , , , − , − , − , , , − − , − , − , − , − , . . . are the nega-tives of 2 n − + n , a shifted-index version of OEIS A005126.Similarly, for T ∗ ( N, m ), we avoid the case m ≡ N ) because it in-creases the order of the recurrence by 1 when N is odd. Here, we must multisectthe sequence according to the congruence classes of m modulo 2 N to accountfor the contribution of the term 2 N in the definition of T ∗ ( N, m ) (20). Table 4gives the values of the coefficients of some of these sequences in difference rootform.
Table 4: Recurrent sequences for multisection of T ∗ ( N, m ) with m ≡ N ) for small values of N N coefficients in a n + a n − + a n − + . . . = 02 1 ,
43 1 ,
274 1 , ,
165 1 , , , , ,
647 1 , , , , , , , , , , , , , , , , Here, we use the standard method of calculating terms of a recurrence by takingpowers of its associated Frobenius companion matrix. The dimension d of thematrix is equal to the number of terms in the tables above, minus 1 (since nocolumn representing the first term is needed); thus d = b ( p + 1) / c for T ( N, m ),and d = b p/ c for T ∗ ( N, m ). We have arbitrarily chosen to use the upper formof the companion matrix, the power of which must be multiplied by a column10ector containing the d starting values of each multisectioned sequence in reverseorder. The required power is b ( p − /N c − d + 1 in the case of T ( N, m ) and b ( p − / (2 N ) c − d + 1 in the case of T ∗ ( N, m ), and the formula is invalid whenthis power is less than 0. Thus, for example, for T ∗ (9 , m ), the expression takesthe form − − − − b p − N c− d +1 · v, where the column vector has the value v = [ − − − − − − − − − − − − p ≡ , , , , ,
17 mod 18.
Several writers have observed that a further reduction of the order of the recur-rences for T ( N, m ) considered above can be accomplished by instead considering J ( N, m ) := N · T ( N, m ) − m , and multisecting the resulting sequences accord-ing to the congruence classes of m modulo N . This technique was exhibited forthe case m ≡ N in a number of OEIS entries by Benoit Cloitre in 2004(OEIS A070775 for N = 4, A070782 for N = 5, A070967 for N = 6, A070832for N = 8, A094213 for N = 9). More recently, it was methodically developedin an important 2014 paper by Russell Jay Hendel [13], where it may be notedthat there is a misprint in equation 1.1 and in the text immediately following(see the paper’s abstract for the correct version). Hendel also created the OEISentry for these “jump sums” (A244608), which links to a table giving the co-efficients in the recursions for N ≤
50. It should also be mentioned that thistechnique is foreshadowed in the congruence for s ∗ (0 ,
9) given by Z.-H. Sun in1993 ([31], pt. 2, Theorem 2.7). In the examples that follow (Table 5), we againavoid the case m ≡ N ) because it is not needed for testing prime valuesof m and increases the order of the recurrence by 1 when N is even; and withthis restriction, the associated recurrences are of order b N − c .11 able 5: Recurrent sequences for J ( N, m ) for small values of N ; cf. OEISA244608 N coefficients in a n + a n − + a n − + . . . = 03 1 ,
14 1 ,
45 1 , , −
16 1 , , −
277 1 , , − , −
18 1 , , − , − , , − , − ,
110 1 , , − , − , , , − , − , ,
112 1 , , − , − , , , , − , − , , , −
114 1 , , − , − , , , − , , , , , , , , , , . . . , arethe Eulerian numbers (OEIS A000295).When m is a prime p , the relationship between Hendel’s sums J ( N, p ) andharmonic numbers is as follows: N − J ( N, p ) − p p ≡ N (1 − T ( N, p )) p ≡ ( s (0 , N ) (mod p ) if N is even s ∗ (0 , N ) (mod p ) if N is odd . (35)There is no corresponding refinement of T ∗ ( N, m ), because substituting N · T ( N, m ) − m into (20) has no effect after the cancellation of opposite terms.The congruence for special harmonic numbers in terms of Hendel’s sums canbe written N − J ( N, p ) − p p = N − J ( N, p ) − − (2 · p − − p ≡ N − − J ( N, p )) p − q p (2) (mod p ) . (36)Thus, Hendel’s result is equivalent to, and provides an independent derivationof, those “classical” congruences for H b p/N c , with N even, that contain a term − q p (2); namely N = 2 , , ,
10 (see congruences 2, 5, 8, 9 above).Regrettably, due to memory overruns, it proved impossible to implementHendel’s algorithm on the machinery available. However, it seems worth record-ing here for its possible value to other researchers.12
The Divisibility of Harmonic Numbers
It may be helpful to distinguish the purpose of the present study with that of themore general problem of the divisibility of Harmonic numbers. Eswarathasanand Levine [9] note that all primes greater than 3 divide the Harmonic numbersof indices p − p ( p − p −
1, and Wolstenholme’s theorem states thatfor a prime p > p | H p − . These results may be inverted to give three generalrules for the divisibility of Harmonic numbers:1. H m is divisible (for even m ) by ( m + 1) if m + 1 is a prime > H m is divisible by √ m +12 if the latter expression is a prime > H m is divisible by √ m + 1 if the latter expression is a prime > harmonic primes as primes that divide onlythe three aforementioned Harmonic numbers (OEIS A092101), and anharmonic primes as those that divide additional Harmonic numbers (OEIS A092102).From the fact that we consider only cases (1) where p | H b p/N c with N > p is anharmonic, and that we are seeking asubset of anharmonic primes that divide Harmonic numbers of relatively smallindex; for example, from Table 10 we see that N = 2 is solved by p = 1093,implying that 1093 | H , while N = 46 is solved by p = 11731, implying that11731 | H . In contrast, the work of Boyd [1] and Rogers [29] entails, in part,finding Harmonic numbers of large index divisible by relatively small primes,such as the case 11 | H . The overlap between ourresults and theirs is slight.Indeed, our results entail only a minority of the anharmonic primes, becausethey require that p divide an Harmonic number H m with m < p −
1, and thusbelong to a subset of the anharmonic primes which has sometimes been calledthe Harmonic irregular primes (see OEIS A092194 and the Wikipedia entry for“Regular prime”). Although the standard definition of such a prime p is thatit divide an Harmonic number H m with n < p −
1, it may be noted that anyprime satisfying this criterion must in fact divide some Harmonic number H m of smaller index m < ( p − /
2, since by symmetry H p − − m ≡ H m mod p . Theinequality sign in this condition is strict because if m = ( p − /
2, then p isa Wieferich prime and by (2) and (4) must likewise divide the smaller H b p/ c .But this inequality also gives a sharp upper bound on the index of the leastHarmonic number divisible by an Harmonic irregular prime p , since among thefirst 999 primes we have m = ( p − / p = 29 , , H m with m ≤ ( p − / p ).In turn, our results entail only a minority of the Harmonic irregular primes,because we must have m = b p/N c for some N . Nonetheless, for some small13rimes this relationship is satisfied by more than one value of N ; for example,with p = 137, p | H , and b p/ c = b p/ c = b p/ c = b p/ c = b p/ c = 5.Taking into account the stringency of the conditions on p , it is perhapsunsurprising that few instances have been found where a given p divides distinctHarmonic numbers of the special type under consideration. By (2) and (4) allWieferich primes p divide both H b p/ c and H b p/ c , but otherwise we have foundonly two instances where a non-Wieferich prime p < , ,
000 divides twodistinct Harmonic numbers of index m in such a way that m = b p/N c is satisfiedfor some N ≤ p = 761, p divides both H (with N = 85through 95) and H (with N = 32, 33), and for p = 845921, p divides both H (with N = 836) and H (with N = 555). Perhaps further solutionsexist with larger N , but we expect such numbers to be rare.There is, however, a sense in which our results do relate to the general prob-lem of the divisibility of Harmonic numbers. While the formulae (2) through(12) are framed in terms of p , they may conversely be seen as divisibility con-ditions on H m . For example, with N = 2 (the Wieferich primes), we seek caseswhere the numerator of H m is divisible by a prime 2 m + 1, 4 m + 1, or 4 m + 3.In the still unresolved case N = 5 we ask whether it is possible for H m to bedivisible by a prime of one of the forms 5 m +1, 5 m +2, 5 m +3, or 5 m +4, and inthe still unresolved case N = 12 whether by a prime of one of the forms 12 m +1,12 m + 5, 12 m + 7, or 12 m + 11. Incidentally, we also tested the numerators of H m for divisibility by any of these linear forms without the restriction that thedivisors be prime, and found only one additional solution for m ≤ , n + 1: H is divisible by 121, which is of course 11 . Forthe same ranges of p as covered by Table 10 below, Table 6 shows, for small N ,all linear forms for which there exists no known p = N · m + r ( r < N ) dividingan Harmonic number H m . 14 able 6: Linear forms for which no known prime p = N · m + r with N ≤
24 divides an Harmonic number H m ; starred rows containall possible forms of p ; p tested to at least 760,000,000 N r ? ?
12 1, 5, 7, 1113 1, 2, 3, 4, 5, 6, 8, 10, 11, 1214 1, 3, 5, 915 1, 2, 4, 8, 11, 13, 1416 5, 7, 9, 11, 1317 1, 2, 3, 4, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16 ?
18 1, 5, 7, 11, 13, 1719 1, 2, 3, 4, 5, 6, 7, 9, 10, 11, 12, 13, 14, 15, 16, 17 ?
20 1, 3, 7, 9, 11, 13, 17, 1921 1, 4, 5, 8, 10, 13, 17, 19, 2022 1, 3, 7, 9, 13, 15, 17, 19, 2123 1, 2, 3, 4, 5, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 2224 11, 13, 23The example noted above of 121 dividing H underscores the fact that if N − N − | H N − and ( N − ≡ N . Thus, although thereare no known prime divisors of H m of the form 12 k + 1, it is not because suchdivisors are algebraically impossible. The following account does not aim to be complete, but we hope that it includesall significant contributions to the problem. In referring to earlier work, forthe sake of brevity we focus mainly on the results of successful searches. Wenote the incorrectness of the statement in [4], pp. 389–390, that H b p/N c ≡ p has a solution p < ,
000 for every N between 2 and 46 other than5, an error perhaps resulting from the inclusion of some vacuous sums with p < N . Indeed, as reported below, for N = 18 , , , ,
43 there are nosolutions with p < , , , N = 12 there is no solution with p < , , , , p of H b p/ c and H b p/ c are the Wieferich primes (OEIS A001220).The first of these, 1093, was found in 1913 by Meissner, and the second, 3511,in 1922 by Beeger. The Wieferich primes have inspired some of the most in-tensive numerical searches ever conducted, but no further instances have beendiscovered. The current search record doubtless belongs to the ongoing test byPrimeGrid [25], which had reached p < , , , , ,
000 as of 1 April2016.The divisors p of H b p/ c are the Mirimanoff primes (OEIS A014127). Thefirst of these, 11, was found by Jacobi in 1828 [16], and the second, 1006003,by Kloss in 1965 [17]. No further instances have been discovered. The currentsearch record appears to be p < , , , , H b p/ c apart from the negativeresult in [4] that it does not vanish modulo p for p < , H b p/ c − H b p/ c . In any case, wehave tested H b p/ c to a much greater limit than that attainable by Schwindt’smethod, without finding any instance where it vanishes modulo p .The divisors of H b p/ c (OEIS A238201, but ignoring solutions with p < p < , p =61. We know of no further computations relating to these numbers until theywere obtained in another guise as the zeros of q p (432) by Richard Fischer [10],whose ongoing test, which had reached p < , , , ,
000 as of 18 August2016, has found two further solutions, included in Tables 9 and 10 below. Thevast range studied by Fischer has not yet been completely retested, but we arecontinuing to replicate his calculations, as to the best of our knowledge they havebeen performed only once. (A similar study by Ležák [22], conducted apparentlywithout awareness of Fischer’s work, extends only to p < , , , The formulae involving Fermat quotients and recurrent sequences cited aboveare all based on exponentiation, and are essentially algorithms of the order log p ,while the direct evaluation of a modular inverse likewise based on exponentiationand so is of order log p , but the calculation of H b p/N c requires the procedure tobe extended over a range proportional to p/N . For those values of H b p/N c notexpressible solely in terms of Harmonic numbers – that is, apart from the caseswith N = 2 , , , N log p . So while several other considerationscome into play, in general terms we are (as noted above) comparing an algorithmof order N log p with one of order p log pN . Thus, roughly speaking, as soon aswe reach values of p > N , the approach involving exponentiation of matricesis to be preferred. As the largest N considered in the main part of this study16s 46, the break-even point for the two methods should be reached by about p = 4 , , N while thesecond is inversely proportional to N , the greatest differences are seen in therun-times for the smallest values of N for which H b p/N c cannot be evaluatedsolely in terms of Fermat quotients and require matrix exponentiation, that is,for N = 5. This study does not attempt to extend existing calculations in the cases N =2 , , ,
6, and relies on known values for those cases in the tables below. AlthoughCikánek [2] also gives similar conditions for the first case of FLT for N up to 94(with certain qualifications), it was decided to concentrate computing resourceson the original range considered by Dilcher and Skula, with N ranging from 2to 46. The limits on N in Figure 1 and in Tables 10 and 11 and are higher,as the values of N between 47 and 52 happen to have fairly small p as leastsolutions, and it was felt they may as well be recorded.This project has been run on varying numbers of CPUs since January 2014,and it would be difficult to state with any precision the number of hours de-voted to individual cases of N . Except for some early exploratory work, allcomputations were performed in PARI, using matrix exponentiation with mod-ular arithmetic to evaluate the required recurrent sequences or matrices, as thecase may be. The search for solutions of p | H b p/N c has now been carried to atleast 383 , ,
000 for all values of N in the range, and solutions have been foundfor all but seven cases out of 45 (see Table 10 below). The search-limits so farattained vary considerably. Those cases for which special formulae exist are asfollows: Table 7: Search limits on divisors p of Harmonic numbers H b p/N c for whicha special formula is known; starred values of N have no knownsolutions N limit of p? ?
12 11,545,400,000,00016 93,700,000,00024 6,691,500,000,000For the remaining values of N the calculations were usually suspended shortlyafter a first solution was obtained, except for N ≤
24, where additional data waswanted for Table 6. These calculations are little different from those involvingspecial formulae, except that the matrices are typically larger. And because17he processing time for matrix exponentiation is proportional to the cube of thedimension of the matrix, the test-limits achieved for the various runs fall off as N increases. The limits so far attained are as follows: Table 8: Search limits on divisors p of Harmonic numbers H b p/N c for whichno special formula is known; starred values of N have no knownsolutions N limit of p
7, 14 8,500,000,0009, ?
18 334,000,000,00011, 22 2,185,000,00013, 15, 17, 19, 21, 23 823,100,000 ?
20 82,880,000,00025–28, 30, 32, 33, 35–42, 44–46 383,950,000 ?
29 82,879,000,000 ?
31 68,720,000,00034 493,000,000 ?
43 27,580,000,000The results of the searches follow. First (Table 9), we report separately thoseresults for which the search-limits are given in Table 2 above, since these haveparticular interest and in some cases coincide with OEIS sequences. Then (Table10) we report all results combined, including those in Table 9.
If one considers the sequence of least divisors p of Harmonic numbers H b p/N c for N = 2 through 52 (most of which are shown in Table 10), arranged in ascendingorder of p and displayed on a logarithmetic scale (see Figure 1), it is readily seenthat the points tend to describe a line of positive curvature. Thus, within thelimits studied, it has been found that the sequence of least divisors p grows ata superexponential rate.Apart from the fact that the solutions for N = 2 and N = 4 coincide, thereis no known number-theoretic reason to doubt that divisors p of H b p/N c arerandomly distributed with respect to N , and the data presented in Table 11provide some empirical support for such a view, yielding a Pearson CorrelationCoefficient between N and p of − . N have at least one solution, then from a statistical point of view the scenariois a “coupon collector’s problem” in which 44 slots (counting N = 2 and N = 4as one) have to be filled with at least one solution, so that the expected numberof solutions (not necessarily minimal) that must be found for all the slots tobe filled is 44 H , or about 192. So far only 78 solutions (some not derivingfrom distinct p ) have been found, but as the test limits for various cases of N vary, we cannot infer much from this fact. However, a rough extrapolation18
10 15 20 25 30 3551015
Figure 1: Least divisors p < , , ,
000 of Harmonic numbers H b p/N c for N = 2 through 52, arranged in ascending order of p ; hori-zontal axis = order of discovery; vertical axis = log p . in Mathematica TM from the distribution of the consecutive minimal solutionssuggests that the 44th minimal solution (if it exists) lies beyond 10 . Thus acomplete solution of the problem would surely lie far beyond the present rangeof computability. Nonetheless, we are continuing the calculations in the hopethat a few more solutions may still be found. We should like to thank the University of Winnipeg Library, and in particularJoffrey Abainza, for providing access to computing resources. We should alsolike to thank T. D. Noe for creating OEIS sequence A238201 based on our [5].19 able 9: Divisors p of Harmonic numbers H b p/N c for which a formula isknown N p
OEIS reference2 1093, 3511 A0012203 11, 1006003 A0141274 1093, 3511 A0012205 —6 61, 1680023, 7308036881 A2382018 269, 8573, 1300709, 11740973, 24107856110 227, 17539, 475015912 —16 38723, 38993, 429254324 137, 577, 247421, 307639, 366019,5262591617, 31251349243
Table 10: Divisors p of Harmonic numbers H b p/N c for all N , 2 through 49 N p N p able 11: Least divisors p < , , ,
000 of Harmonic numbers H b p/N c for N = 2 through 52, arranged in ascending order of p p N
11 361 6137 23, 24, 25, 26, 27227 10269 8509 40, 41, 42521 19677 9761 32, 331093 2, 41297 361423 211439 371553 342113 47, 482267 142473 382843 224127 514139 356967 447537 5211731 4620101 2838723 16134227 15246277 11407893 39609221 45652913 73337469 505611043 4927089407 3043214711 13590422517 1721 eferences [1] David W. Boyd, “A p -adic study of the partial sums of the harmonic series,”Experiment. Math. (1994), 287–302.[2] Petr Cikánek, “A special extension of Wieferich’s criterion,” Math. Comp. (1994), 923–930.[3] K. Dilcher & L. Skula, “Linear relations between certain sums of reciprocalsmodulo p ,” Ann. Sci. Math. Québec (2011), 17–29.[4] Karl Dilcher & Ladislav Skula, “A new criterion for the first case of Fermat’sLast Theorem,” Math. Comp. (1995), 363–392.[5] John Blythe Dobson, “Extended calculations of a special Harmonic number,”preprint at http://arxiv.org/abs/1402.5680 .[6] John Blythe Dobson, “A matrix variation on Ramus’s identity for lacu-nary sums of binomial coefficients,” International Journal of Mathemat-ics and Computer Science 12 (2017), 27–42, available online at http://ijmcs.future-in-tech.net/12.1/R-Dobson.pdf .[7] F. G. Dorais and D. Klyve, “A Wieferich prime search up to p < . ∗ ,”J. Integer Seq. 14 (2011), Art. 11.9.2, 1-14. .[8] [G.] Eisenstein, “Eine neue Gattung zahlentheoretischer Funktionen, welchevon zwei Elementen abhängen und durch gewisse lineare Funktional-Gleichungen definirt werden,” Berichte Königl. Preuß. Akad. Wiss. Berlin (1850), 36–42.[9] Arulappah Eswarathasan and Eugene Levine, “ p -integral harmonic sums,”Discrete Mathematics (1991), 249–257.[10] Richard Fischer, “Fermatquotient b ( p − ≡ p ),” available online at .[11] J. W. L. Glaisher, “A general congruence theorem relating to the Bernoul-lian function,” Proc. London Math. Soc. (1900), 27–56.[12] J. W. L. Glaisher, “On the residues of r p − to modulus p , p , etc.,” Q. J.Math. Oxford (1900–1901), 1–27.[13] Russell Jay Hendel, “A Cayley-Hamilton and circulant approach to jumpsums,” Proceedings of the 16th International Conference on Fibonacci Num-bers and Their Applications (= Fibonacci Quarterly, vol. 52, no. 5, December2014), 12–135.[14] V. E. Hoggatt, Jr., and G. L. Alexanderson, “Sums of partition sets ingeneralized Pascal triangles, I” [no more published?], Fib. Quart. (1976),117–125. 2215] F. T. Howard and Richard Witt, “Lacunary sums of binomial coefficients,”Applications of Fibonacci Numbers (1998), 185–95.[16] C. G. J. Jacobi, “Beantwortung der aufgabe S. 212 dieses Bandes: «Kann α µ − , wenn µ eine Primzahl und α eine ganze Zahl und kleiner als µ undgrößer als 1 ist, durch µµ theilbar sein?»,” J. reine u. angew. Math. (1828),301–303.[17] K. E. Kloss, “Some Number-Theoretic Calculations,” Journal of Researchof the National Bureau of Standards — B. Mathematics and MathematicalPhysics (1965), 335–336.[18] J. Konvalina and Y.-H. Liu, “Arithmetic progression sums of binomial co-efficients,” Appl. Math. Lett. (1997), 11–13.[19] Takako Kuzumaki & Jerzy Urbanowicz, “On congruences for certain bi-nomial coefficients of E. Lehmer’s type,” in Number Theory: Arithmetic inShangri-La. Proceedings of the 6th China–Japan Seminar, Shanghai, China,15–17 August 2011 , 132–140.[20] Emma Lehmer. “On Congruences involving Bernoulli numbers and the quo-tients of Fermat and Wilson,” Ann. of Math. (1938), 350–360.[21] M. Lerch, “Zur Theorie des Fermatschen Quotienten. . . ,” Math. Ann. (1905), 471–490.[22] Petr Ležák, “Hledání Wieferichových Prvočísel,” Kvaternion (2013)103–109, and supplementary computer file at http://sourceforge.net/projects/wieferich/files/results/table.txt .[23] D. Mirimanoff, “Sur le dernier théorème de Fermat,” C. R. Acad. Sci. Paris, (1910) 204–206; also “Sur le dernier théorème de Fermat,” J. reine u.angew. Math., (1911), 309–324.[24] PrimeGrid, “Wall-Sun-Sun Prime Search,” available online at http://prpnet.primegrid.com:13001/ .[25] PrimeGrid, “Wieferich Prime Search,” available online at http://prpnet.primegrid.com:13000/ .[26] [C.] Ramus, “Solution générale d’un problème d’analyse combinatoire,” J.reine u. angew. Math. (1834), 353–55.[27] Paulo Ribenboim,
13 Lectures on Fermat’s Last Theorem . New York, etc.:Springer, 1979.[28] Paulo Ribenboim, “Some Criteria for the First Case of Fermat’s Last The-orem,” Tokyo J. Math. (1978), 149–155.[29] Mathew D. Rogers, “An extension of Boyd’s p -adic algorithm for the har-monic series,” preprint at http://arxiv.org/abs/0708.2439 .2330] H. Schwindt, “Three summation criteria for Fermat’s Last Theorem,”Math. Comp. (1983), 715–716.[31] Zhi-Hong Sun, “[The] combinatorial sum P k ≡ r (mod m ) (cid:0) nk (cid:1) and its ap-plications in Number Theory” (in Chinese), J. Nanjing Univ. Math. Bi-quarterly (1992): 227–240, (1993): 205–118, (1995): 90-102.A very full summary in English is available on the author’s website, at .[32] Zhi-Hong Sun & Zhi-Wei Sun, “Fibonacci numbers and Fermat’s Last The-orem,” Acta Arith. (1992), 371–388.[33] Zhi-Wei Sun, “On sums of binomial coefficients and their applications,”Discrete Mathematics (2008), 4231–4245.[34] Zhi-Wei Sun, “On the sum P k ≡ r (mod m ) (cid:0) nk (cid:1) and related congruences,”Israel Journal of Mathematics (2002), 135–156.[35] H. S. Vandiver, “Extension of the criteria of Wieferich and Mirimanoff inconnection with Fermat’s last theorem,” J. reine u. angew. Math. (1914),314–318.[36] A. Wieferich, “Zum letzten Fermatschen Theorem,” Journal für dieReine und Angewandte Mathematik (1909), 293–302. Available on-line at http://gdzdoc.sub.uni-goettingen.de/sub/digbib/loader?ht=VIEW&did=D255696 .[37] H. C. Williams, “Some formulas concerning the fundamental unit of a realquadratic field,” Discrete Math.92