The Størmer Problem (Lehmer's method etc)
 

It makes sense to start a specific thread for this topic, I think ..

The discussion developed originally in this thread
[URL]http://www.mersenneforum.org/showthread.php?t=5647[/URL]

I'll start this thread by reproducing my first posting there, which has some useful data, and also give a simple definition of the problem for newcomers:


[B][U]Størmer's Problem (1897)[/U][/B]

[SIZE=3][COLOR=black]Write [I][B]gpd[/B][/I]([I][B]n[/B][/I]) for the greatest prime divisor of [I][B]n[/B][/I][/COLOR][/SIZE]


[COLOR=black][SIZE=3]Given some prime [B]p[/B], [/SIZE][SIZE=3]for which [I][B]n[/B][/I] do we have [I][B]gpd[/B][/I]([B]n[/B]([B]n[/B]-[B]1[/B])) <= [B]p[/B]? [/SIZE][/COLOR]


We can refer to the set of all such values of [B][I]n [/I][/B]for given [B]p [/B]as [FONT=Courier New][SIZE=3][COLOR=#000080][FONT=Verdana][SIZE=2][COLOR=black][B]S(p)[/B], and [/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][SIZE=3][COLOR=#000080][SIZE=2][COLOR=black][B]S'[/B]([B]p[/B])[/COLOR][/SIZE] [/COLOR][/SIZE]for the [I]maximum[/I] value [B][I]n [/I][/B]for which [B][I]gpd[/I][/B]([B]n[/B]([B]n[/B]-[B]1[/B])) [B]=[/B] [B]p [/B](note that these sets are demonstrably [I]finite, [/I]and so [B]S'[/B]([B]p[/B]) always exists).

The pairs ([B]n[/B],[B] n-1[/B]) are obviously representable as (twice the value of) triangular numbers, via their product, and are also known in musical interval theory as [I]Superparticular Ratios, [/I]when represented as rational quotients [B]n / n-1 [/B](these ratios are highly composite, and are as close to possible to 1)

=======================================================

OEIS sequences [B]A002072[/B] and [B]A117581[/B] both list the sequence [B]S'(p)[/B], only the first one lists the lower value of each pair, ie [B]S'(p) - 1[/B]

Using the higher value is preferable, I think, since the only published tables, the 1964 results of Dick Lehmer ("On a Problem of Størmer", [I]Illinois J. Math[/I], [B]8[/B], 1964, pp 57-79), for [B]p [/B]up to 41,use that convention.

Also his proofs are given in the same terms, and it corresponds to the convention used by the muso's.

Both sequences as recorded in OEIS have an inconvenient feature, though - there are duplicate values in the sequence, in order (presumably) to keep the sequence increasing.

But this omits useful information. If we ask what is [B][I]S'[/I](23)[/B], the [I]greatest[/I] value of S for which [I]gpd [/I](S(S-1)) = 23?

The answer is 5142501, but this is [I]less[/I] than [B][I]S'[/I](19) = [/B]11859211, so it does not appear in the sequence, the latter value is just repeated.

Here is a list which applies the definition of [B][I]S'[/I](p) [/B]strictly, and which extends the sequence to the 35th prime (149):

[code]
[SIZE=3][FONT=Courier New][COLOR=navy]N p S'(p) log2(S')[/COLOR][/FONT][/SIZE]
[SIZE=3][FONT=Courier New][COLOR=navy]=============================================[/COLOR][/FONT][/SIZE]
[SIZE=3][FONT=Courier New][COLOR=navy] 1. 2 2 1[/COLOR][/FONT][/SIZE]
[SIZE=3][FONT=Courier New][COLOR=navy] 2. 3 9 3.1699[/COLOR][/FONT][/SIZE]
[SIZE=3][FONT=Courier New][COLOR=navy] 3. 5 81[/COLOR][/FONT][/SIZE]
[SIZE=3][FONT=Courier New][COLOR=navy] 4. 7 4375[/COLOR][/FONT][/SIZE]
[SIZE=3][FONT=Courier New][COLOR=navy] 5. 11 9801[/COLOR][/FONT][/SIZE]
[SIZE=3][FONT=Courier New][COLOR=navy] 6. 13 123201[/COLOR][/FONT][/SIZE]
[SIZE=3][FONT=Courier New][COLOR=navy] 7. 17 336141[/COLOR][/FONT][/SIZE]
[SIZE=3][FONT=Courier New][COLOR=navy] 8. 19 11859211 23.4995[/COLOR][/FONT][/SIZE]
[SIZE=3][FONT=Courier New][COLOR=navy] 9. 23 5142501 22.2940[/COLOR][/FONT][/SIZE]
[SIZE=3][FONT=Courier New][COLOR=navy]10. 29 177182721 27.4077[/COLOR][/FONT][/SIZE]
[SIZE=3][FONT=Courier New][COLOR=navy]11. 31 1611308700 30.5856 [/COLOR][/FONT][/SIZE]
[SIZE=3][FONT=Courier New][COLOR=navy]12. 37 3463200000 31.6895[/COLOR][/FONT][/SIZE]
[SIZE=3][FONT=Courier New][COLOR=navy]13. 41 63927525376 35.8957[/COLOR][/FONT][/SIZE]
[SIZE=3][FONT=Courier New][COLOR=navy]14. 43 421138799640 38.6155[/COLOR][/FONT][/SIZE]
[SIZE=3][FONT=Courier New][COLOR=navy]15. 47 1109496723126 40.0103[/COLOR][/FONT][/SIZE]
[SIZE=3][FONT=Courier New][COLOR=navy]16. 53 1453579866025 40.4027[/COLOR][/FONT][/SIZE]
[SIZE=3][FONT=Courier New][COLOR=navy]17. 59 20628591204481 44.2297[/COLOR][/FONT][/SIZE]
[SIZE=3][FONT=Courier New][COLOR=navy]18. 61 31887350832897 44.8580[/COLOR][/FONT][/SIZE]
[SIZE=3][FONT=Courier New][COLOR=navy]19. 67 12820120234376 43.5435[/COLOR][/FONT][/SIZE]
[SIZE=3][FONT=Courier New][COLOR=navy]20. 71 119089041053697 46.7090[/COLOR][/FONT][/SIZE]
[SIZE=3][FONT=Courier New][COLOR=navy]21. 73 2286831727304145 51.0223[/COLOR][/FONT][/SIZE]
[SIZE=3][FONT=Courier New][COLOR=navy]22. 79 9591468737351909376 63.0565 [/COLOR][/FONT][/SIZE]
[SIZE=3][FONT=Courier New][COLOR=navy]23. 83 17451620110781857 53.9542 [/COLOR][/FONT][/SIZE]
[SIZE=3][FONT=Courier New][COLOR=navy]24. 89 166055401586083681 57.2044[/COLOR][/FONT][/SIZE]
[SIZE=3][FONT=Courier New][COLOR=navy]25. 97 49956990469100001 55.4715[/COLOR][/FONT][/SIZE]
[SIZE=3][FONT=Courier New][COLOR=navy]26. 101 4108258965739505500 61.8332[/COLOR][/FONT][/SIZE]
[SIZE=3][FONT=Courier New][COLOR=navy]27. 103 19316158377073923834001 74.0322 [/COLOR][/FONT][/SIZE]
[SIZE=3][FONT=Courier New][COLOR=navy]28. 107 386539843111191225 58.4234[/COLOR][/FONT][/SIZE]
[SIZE=3][FONT=Courier New][COLOR=navy]29. 109 90550606380841216611 66.2954 [/COLOR][/FONT][/SIZE]
[SIZE=3][FONT=Courier New][COLOR=navy]30. 113 205142063213188103640 67.4752[/COLOR][/FONT][/SIZE]
[SIZE=3][FONT=Courier New][COLOR=navy]31. 127 53234795127882729825 65.5290[/COLOR][/FONT][/SIZE]
[SIZE=3][FONT=Courier New][COLOR=navy]32. 131 4114304445616636016032 71.8011[/COLOR][/FONT][/SIZE]
[SIZE=3][FONT=Courier New][COLOR=navy]33. 137 124225935845233319439174 76.7173[/COLOR][/FONT][/SIZE]
[SIZE=3][FONT=Courier New][COLOR=navy]34. 139 3482435534325338530940 71.5606[/COLOR][/FONT][/SIZE]
[SIZE=3][FONT=Courier New][COLOR=navy]35. 149 6418869735252139369210 72.4428[/COLOR][/FONT][/SIZE]
[/code]


Jim White
Math. Sciences Institute
Australian National University

Here is an entry from the table above, showing the factors of the corresponding pair:


[B][FONT=Courier New][COLOR=navy][COLOR=black]124225935845233319439173[/COLOR] [/COLOR][/FONT][/B]
[B][FONT=Courier New][COLOR=navy] 13^6, 17, 37, 53, 67, 79, 83, 101, 127, 137[/COLOR][/FONT][/B]
[B][FONT=Courier New][/FONT][/B]
[B][FONT=Courier New][COLOR=navy][COLOR=black]124225935845233319439174[/COLOR][/COLOR][/FONT][/B]
[B][FONT=Courier New][COLOR=navy]2, 3^4, 7, 19^2, 23^2, 29, 31, 47, 97^2, 113^3[/COLOR][/FONT][/B]

There are no greater values possible for 137-smooth pairs of consecutive integers.

Lehmer's Method
 

Much of what follows will discuss "Lehmer's Method", which finds all members of [B]S(p) [/B]by solving 2^n Pell equations (where n is the number of primes involved), and testing the solutions thereof for [B]p[/B]-smoothness.

So we assume from here on that readers have read that paper and/or know exactly what we're talking about!

There is a brute-force method that can also be used that requires no knowledge of Pell equations or continued fractions, but requires that you have a reasonable upper bound on [B]S'(p)[/B]

Both methods have exponential complexity however, and in reality we have great difficulty extending the known results, as identifying each additional value of the sequence [B]S'(p) [/B]can take anywhere from 2 to around 10 times as long as the previous one.

Currently we seem to be limited in what can be done in a reasonable time to around 25-30 primes.


I'll quickly summarise the two methods in a post below, and what's known about [I]proof-of-correctness [/I]issues

Lehmer method in a Nutshell (updated)
 

[LIST][*]let [B][I][SIZE=3]Q[/SIZE][/I][/B] be the set of primes {2, ... [I][SIZE=3]p[/SIZE][/I]}, and let [I][FONT=Times New Roman][SIZE=4]n[/SIZE][/FONT][/I][B] [I]= [/I][/B][tex]\pi[/tex][SIZE=3]([I]p[/I])[/SIZE] =[SIZE=3][I]|[B]Q[/B]| [/I][/SIZE][SIZE=2]be the number of primes in [/SIZE][SIZE=3][B][I]Q[/I][/B][/SIZE][/LIST][LIST][*]let [I][B][SIZE=3]Q[/SIZE]*[/B] [/I]denote the set of [I][SIZE=3]p[/SIZE][/I]-smooth numbers[/LIST][LIST][*]let [I][B][SIZE=3]Q'[/SIZE][/B] [/I]be the set of all [I]odd, square-free[/I] members of [B][I][SIZE=3]Q[/SIZE]*, [/I][/B]so [SIZE=3]|[B][I]Q'[/I][/B]|[/SIZE] = [tex]2^{n-1}[/tex][/LIST][LIST][*]let [B][I][SIZE=3]S[/SIZE][/I]([/B][I][SIZE=3]p[/SIZE][/I][B]) [/B]be the set of all integers [tex]a[/tex] for which [tex]a({a-1}) \in [/tex] [B][I][SIZE=3]Q[/SIZE]*[/I][/B][/LIST][B][I]Algorithm to Enumerate [SIZE=3]S[/SIZE][/I]([/B][SIZE=3][I]p[/I][/SIZE][B]) ("[I]Lehmer's Method"[/I])[/B][LIST][*]for every [tex]D \in[/tex] [I][B][SIZE=3]Q'[/SIZE][/B], [/I]find [tex] (x_1, y_1) [/tex], the fundamental solution to the Pell equation [tex] x^2 - 2Dy^2 = 1 [/tex]. If ... [tex]y_1[/tex] ... is [SIZE=3][I]p[/I][/SIZE][B]-[/B]smooth, then [tex]x_1= 2a - 1[/tex], and both [I][SIZE=3][FONT=Times New Roman][SIZE=4]a[/SIZE][/FONT] [/SIZE][/I]and [I][SIZE=3][FONT=Times New Roman][SIZE=4]a[/SIZE][/FONT] [/SIZE][/I]- [SIZE=3]1[/SIZE] are [I][SIZE=3]p[/SIZE][/I][B]-[/B]smooth, ie [tex]a \in[/tex] [B][I][SIZE=3]S[/SIZE][/I]([/B][I][SIZE=3]p[/SIZE][/I][B])[/B][/LIST][LIST][*]repeat the step above for 4D, ie solving [tex]x^2 - 4Dy^2 = 1[/tex][/LIST][LIST][*]in both cases, should the Pell equation produce a [I][SIZE=3]p[/SIZE][/I][B]-[/B]smooth result, then also inspect the multiple-solution sequence [tex](x_k, y_k)[/tex], for [tex] k = {2, 3 ... max(3, (p+1)/2)}[/tex]. Check each [tex]y_k[/tex] ... for [I][SIZE=3]p[/SIZE][B]-[/B][/I]smoothness, and the corresponding [tex]x_k[/tex]will also produce a member of [B][I][SIZE=3]S[/SIZE][/I]([/B][I][SIZE=3]p[/SIZE][/I][B])[/B][/LIST][I][U]Remarks[/U][/I]

Remember, there is no need to generate and inspect the solution sequence unless you have got a "hit" with the fundamental solution. That's because ... [tex]y_1[/tex] divides all other ... [tex]y_k[/tex] ... so if it isn't smooth none of them can possibly be smooth.

There are both 1st and 2nd order recurrence relations for iterating the multiple solutions:

[tex]x_{n+1} = {x_1}{x_n} + Dy_1{y_n}[/tex]
[tex]y_{n+1} = {x_1}{y_n} + x_n{y_1}[/tex]

[tex]x_{n+1} = 2x_1{x_n} - x_{n-1}[/tex]
[tex]y_{n+1} = 2x_1{y_n} - y_{n-1}[/tex]

Størmer's Original Solution
 

[B]Carl Størmer [/B]must be given full credit for his demonstration in 1897, that Pell equation properties allowed a finite characterisation of the problem, and came up with an algorithm and the proof of its correctness, which constituted the first ever proof that [B][I][SIZE=3]S[/SIZE][/I]([/B][I][SIZE=3]p[/SIZE][/I][B]) [/B]was in fact always finite[B].[/B]

The properties of the sequences of multiple solutions of a Pell equation were not well known then, and Størmer gave the original method and proof in terms of fundamental solutions alone. A larger set of [tex]D[/tex] than the one we used above is shown to produce all the required values as fundamental solutions.

Størmer showed that solving [tex]x^2 - Dy^2 = 1[/tex] for all [tex]D \in[/tex] [B][I][SIZE=3]Q''[/SIZE][/I][/B], where [SIZE=3][B][I]Q'' [/I][/B][/SIZE]is the set of all [I]cube-free [/I]members of [B][I][SIZE=3]Q*[/SIZE][/I][/B] would necessarily produce all [tex]a \in[/tex] [B][I][SIZE=3]S[/SIZE][/I]([/B][I][SIZE=3]p[/SIZE][/I][B]) [/B]in some solution as [tex]x_1 = 2a-1[/tex]

It boils down to showing that every target value will turn up as fundamental solution for any [tex]D[/tex] that correctly matches the target's prime composition, and which also gets the exponent [I]parities[/I] right too.

The number of equations that need to be checked, though, is much greater than in Lehmer's later refinement, as we have [SIZE=3]|[B][I]Q''[/I][/B]|[/SIZE] = [tex]3^n - 2^n[/tex]

Mind you, even with Lehmer's reduction of the problem to [tex]O(2^n)[/tex], this is not a cheap computation by any means. Lehmer was able to do the first 13 primes on a 1962-vintage IBM704, and with todays 2GHz desktop box we can only advance that by 20 or so primes before hitting the wall ...

As the number of equations to be solved doubles with each prime, so also does the individual cost of each equation, as the magnitudes of the integers involved spirals ever upward.

In reality, the expected running time (with the Lehmer method), appears to be more like [tex]O(10^n)[/tex]

[QUOTE=Jim White;113263][*]in both cases, should the Pell equation produce a [I][SIZE=3]p[/SIZE][/I][B]-[/B]smooth result, then also inspect the multiple-solution sequence [tex](x_k, y_k)[/tex], for [tex] k = {2, 3 ... max(3, (p+1)/2)}[/tex].[/QUOTE]

Why is it need to check up to only max(3,(p+1)/2) ? Is it trivial, because I don't see it.

You can send the extensions of the Sloane's sequences to Neil Sloane as a "b-file", see: [URL="http://www.research.att.com/~njas/sequences/Submit.html"]http://www.research.att.com/~njas/sequences/Submit.html[/URL]

If some one is interested in writing the code, may be we can search for more numbers in this series.

[url]http://www.mersenneforum.org/showthread.php?t=5630[/url]
[url]http://www.research.att.com/~njas/sequences/A116486[/url]

Thanks.:smile:

A Brute-force Method
 

The only alternative method for complete enumeration appears to be a brute-force search of all [tex]S \in[/tex] [B][I][SIZE=3]Q[/SIZE]*[/I][/B], a search for which we will obviously need a [I]bound,[/I] [tex]S \leq X[/tex]

Any results are only provably correct if the bound [tex]X[/tex] is. So we can use known existing bounds, otherwise we'd need something like a bound produced by a lattice basis reduction algorithm.

The method then is to enumerate [B][I][SIZE=3]Q[/SIZE]* [/I][/B]recursively, using the bound to limit the depth. For each [tex]S \in[/tex] [B][I][SIZE=3]Q[/SIZE][SIZE=2]* [/SIZE][/I][/B]we simply test [tex]S-1[/tex] for [I][SIZE=3]p-[/SIZE][/I]smoothness.

Crude as it sounds, and also bearing in mind that our target subset [B][I][SIZE=3]S[/SIZE][/I]([/B][I][SIZE=3]p[/SIZE][/I][B]) [/B]is but a tiny fraction of the bounded subset of [I][B][SIZE=3]Q[/SIZE]*[/B], [/I]in reality this method actually runs much faster than Lehmer's method.

On a 2GHz AMD64 (with GMP for MP-arithmetic where necessary), 17 primes (p \leq 59) takes about 4 hours using Lehmer's method, while the generative search method (armed with the known bound), gets there in just 20 seconds.

Getting to the 18th prime, 61, with the Lehmer method would probably take 30 to 40 hours. The generative search can get to [B]25 [/B]primes (p \leq 97) in about an hour. Once again, we assume a tight bound is available, but we can assume, from what (little) I understand about the "LLL-BRA" method, that this is generally available, and at relatively low cost.

The running time for the generative search appears [tex]O(2^n)[/tex], but only up to the point where it can all be done with native 64-bit integer arithmetic. I expect that the log gradient will begin to rise after that.

I have not yet got enough data to estimate the trend beyond that point (around p = 97). There is a large "spike" at prime 27 = 103, where the upper bound has to be raised from [tex]2^{64}[/tex] to [tex]2^{75}[/tex]. That run I have done, and it took just under 4 hours.

These spikes (due to the not-always-increasing sequence of maximum solutions) make it hard to predict any more, without some rather long and tedious test runs.


I have a preliminary chart of these times, about somewhere, I'll attach it here. The times here are for 'incremental' runs - each run was targeting just those values where the maximum [B]p [/B]actually divides the solution

[quote=Citrix;113318]
[URL]http://www.mersenneforum.org/showthread.php?t=5630[/URL]
[URL]http://www.research.att.com/~njas/sequences/A116486[/URL]
[/quote]

Yes, there are obvious similarities ...

... and a few potentially significant differences perhaps :unsure:


Are these "logarithmically smooth" numbers purely a curiosity, or do they have some application?

Cheers

[quote=R. Gerbicz;113313]You can send the extensions of the Sloane's sequences to Neil Sloane as a "b-file", see: [URL]http://www.research.att.com/~njas/sequences/Submit.html[/URL][/quote]

I plan to do that, also have a very interesting new proposed sequence relating to maximal orthogonal bases in the corresponding smooth-rational vector space

[quote=R. Gerbicz;113313]Why is it need to check up to only max(3,(p+1)/2) ? Is it trivial, because I don't see it.[/quote]

Not difficult, but not obvious, the Pell equation multiple-solution sequences [I]xn, yn [/I]are simple 1st- or 2nd- order recurrence sequences, and Lucas's theorems relating include a "Law of (prime) Apparition" which says that the [tex]n[/tex]-th term in any such sequence must have at least one prime divisor [tex]\geq 2n - 1[/tex]

Incremental Enumeration with Lehmer Method
 

The [I]incremental [/I]version of the problem is:


[COLOR=black][SIZE=3]Given some prime [B]p[/B], [/SIZE][SIZE=3]for which [I][B]n[/B][/I] do we have [I][B]gpd[/B][/I]([B]n[/B]([B]n[/B]-[B]1[/B])) = [B]p[/B]? [/SIZE][/COLOR]


That is, we replace the inequality of the original problem definition with equality.


In the brute-force method, adapting the generative search to perform incrementally is straight-forward (we only visit/generate nodes that include [B]p[/B]).

For the Pell equation approach, however, we need to adapt Lehmer's method slightly.

We look at [tex]D[/tex] values corresponding to each of the [tex]2^{n-1}[/tex] combinations of primes [I]not[/I] including the new [tex]p[/tex], and for each such [tex]D[/tex] we check the Pell equations for both [tex]pD[/tex] and [tex]p^2D[/tex], thus covering both solution parity possibilities for the new divisor.

This modification is, I think, provably correct, and follows from Størmer's Theorem (and works in practice!). It also suggests that multiple solution checking is probably going to be largely unnecessary - to what exact extent I haven't quite got around to figuring out just yet! That's because it really has no significant impact on the running time, which is dominated by the task of obtaining fundamental solutions.

The result generally runs, as expected, in around half the time that it would take to generate the complete set for all [B]p'[/B] <= [B]p[/B] , so it is useful to have.