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

Jim White · 2007-08-30, 23:43

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

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

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:

Størmer's Problem (1897)

Write gpd(n) for the greatest prime divisor of n

Given some prime p, for which n do we have gpd(n(n-1)) <= p?

We can refer to the set of all such values of n for given p as S(p), and S'(p) for the maximum value n for which gpd(n(n-1)) = p (note that these sets are demonstrably finite, and so S'(p) always exists).

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

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

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

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", Illinois J. Math, 8, 1964, pp 57-79), for p 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 S'(23), the greatest value of S for which gpd (S(S-1)) = 23?

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

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

Code:

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

Jim White
Math. Sciences Institute
Australian National University

Jim White · 2007-08-30, 23:58

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

124225935845233319439173
13^6, 17, 37, 53, 67, 79, 83, 101, 127, 137

124225935845233319439174
2, 3^4, 7, 19^2, 23^2, 29, 31, 47, 97^2, 113^3

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

Jim White · 2007-08-31, 00:11

Much of what follows will discuss "Lehmer's Method", which finds all members of S(p) by solving 2^n Pell equations (where n is the number of primes involved), and testing the solutions thereof for p-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 S'(p)

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 S'(p) 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 proof-of-correctness issues

Jim White · 2007-08-31, 09:52

let Q be the set of primes {2, ... p}, and let n = $\pi$ (p) =|Q| be the number of primes in Q

let Q* denote the set of p-smooth numbers

let Q' be the set of all odd, square-free members of Q*, so |Q'| = $2^{n-1}$

let S(p) be the set of all integers $a$ for which $a({a-1}) \in$ Q*

Algorithm to Enumerate S(p) ("Lehmer's Method")

for every $D \in$ Q', find $(x_1, y_1)$ , the fundamental solution to the Pell equation $x^2 - 2Dy^2 = 1$ . If ... $y_1$ ... is p-smooth, then $x_1= 2a - 1$ , and both a and a - 1 are p-smooth, ie $a \in$ S(p)

repeat the step above for 4D, ie solving $x^2 - 4Dy^2 = 1$

in both cases, should the Pell equation produce a p-smooth result, then also inspect the multiple-solution sequence $(x_k, y_k)$ , for $k = {2, 3 ... max(3, (p+1)/2)}$ . Check each $y_k$ ... for p-smoothness, and the corresponding $x_k$ will also produce a member of S(p)

Remarks

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 ... $y_1$ divides all other ... $y_k$ ... 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:

$x_{n+1} = {x_1}{x_n} + Dy_1{y_n}$
$y_{n+1} = {x_1}{y_n} + x_n{y_1}$

$x_{n+1} = 2x_1{x_n} - x_{n-1}$
$y_{n+1} = 2x_1{y_n} - y_{n-1}$

Jim White · 2007-08-31, 10:40

Carl Størmer 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 S(p) was in fact always finite.

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 $D$ than the one we used above is shown to produce all the required values as fundamental solutions.

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

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

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

Mind you, even with Lehmer's reduction of the problem to $O(2^n)$ , 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 $O(10^n)$

R. Gerbicz · 2007-09-01, 08:26

Quote:

Originally Posted by Jim White

[*]in both cases, should the Pell equation produce a p-smooth result, then also inspect the multiple-solution sequence $(x_k, y_k)$ , for $k = {2, 3 ... max(3, (p+1)/2)}$ .

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: http://www.research.att.com/~njas/sequences/Submit.html

Citrix · 2007-09-01, 17:02

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

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

Thanks.

Jim White · 2007-09-01, 22:13

The only alternative method for complete enumeration appears to be a brute-force search of all $S \in$ Q*, a search for which we will obviously need a bound, $S \leq X$

Any results are only provably correct if the bound $X$ 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 Q* recursively, using the bound to limit the depth. For each $S \in$ Q* we simply test $S-1$ for p-smoothness.

Crude as it sounds, and also bearing in mind that our target subset S(p) is but a tiny fraction of the bounded subset of Q*, 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 25 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 $O(2^n)$ , 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 $2^{64}$ to $2^{75}$ . 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 p actually divides the solution

Jim White · 2007-09-01, 22:23

Quote:

Originally Posted by Citrix

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

Yes, there are obvious similarities ...

... and a few potentially significant differences perhaps

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

Cheers

Jim White · 2007-09-01, 22:55

Quote:

Originally Posted by R. Gerbicz

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

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:

Originally Posted by R. Gerbicz

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

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

Jim White · 2007-09-01, 23:35

The incremental version of the problem is:

Given some prime p, for which n do we have gpd(n(n-1)) = p?

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 p).

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

We look at $D$ values corresponding to each of the $2^{n-1}$ combinations of primes not including the new $p$ , and for each such $D$ we check the Pell equations for both $pD$ and $p^2D$ , 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 p' <= p , so it is useful to have.

2007-08-30, 23:43	#1
Jim White Aug 2007 Canberra, Australia 3×11 Posts	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 http://www.mersenneforum.org/showthread.php?t=5647 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: Størmer's Problem (1897) Write gpd(n) for the greatest prime divisor of n Given some prime p, for which n do we have gpd(n(n-1)) <= p? We can refer to the set of all such values of n for given p as S(p), and S'(p) for the maximum value n for which *gpd(n(n-1)) = p (note that these sets are demonstrably finite,* and so S'(p) always exists). The pairs (n, n-1) are obviously representable as (twice the value of) triangular numbers, via their product, and are also known in musical interval theory as Superparticular Ratios, when represented as rational quotients n / n-1 (these ratios are highly composite, and are as close to possible to 1) ======================================================= OEIS sequences A002072 and A117581 both list the sequence S'(p), only the first one lists the lower value of each pair, ie S'(p) - 1 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", Illinois J. Math, 8, 1964, pp 57-79), for p 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 S'(23), the greatest value of S for which gpd (S(S-1)) = 23? The answer is 5142501, but this is less than S'(19) = 11859211, so it does not appear in the sequence, the latter value is just repeated. Here is a list which applies the definition of S'(p) strictly, and which extends the sequence to the 35th prime (149): Code: N p S'(p) log2(S') ============================================= 1. 2 2 1 2. 3 9 3.1699 3. 5 81 4. 7 4375 5. 11 9801 6. 13 123201 7. 17 336141 8. 19 11859211 23.4995 9. 23 5142501 22.2940 10. 29 177182721 27.4077 11. 31 1611308700 30.5856 12. 37 3463200000 31.6895 13. 41 63927525376 35.8957 14. 43 421138799640 38.6155 15. 47 1109496723126 40.0103 16. 53 1453579866025 40.4027 17. 59 20628591204481 44.2297 18. 61 31887350832897 44.8580 19. 67 12820120234376 43.5435 20. 71 119089041053697 46.7090 21. 73 2286831727304145 51.0223 22. 79 9591468737351909376 63.0565 23. 83 17451620110781857 53.9542 24. 89 166055401586083681 57.2044 25. 97 49956990469100001 55.4715 26. 101 4108258965739505500 61.8332 27. 103 19316158377073923834001 74.0322 28. 107 386539843111191225 58.4234 29. 109 90550606380841216611 66.2954 30. 113 205142063213188103640 67.4752 31. 127 53234795127882729825 65.5290 32. 131 4114304445616636016032 71.8011 33. 137 124225935845233319439174 76.7173 34. 139 3482435534325338530940 71.5606 35. 149 6418869735252139369210 72.4428 Jim White Math. Sciences Institute Australian National University Last fiddled with by Jim White on 2007-08-30 at 23:53

2007-08-31, 00:11	#3
Jim White Aug 2007 Canberra, Australia 41₈ Posts	Lehmer's Method Much of what follows will discuss "Lehmer's Method", which finds all members of S(p) by solving 2^n Pell equations (where n is the number of primes involved), and testing the solutions thereof for p-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 S'(p) 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 S'(p) 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 proof-of-correctness issues Last fiddled with by Jim White on 2007-08-31 at 00:14

2007-08-31, 09:52	#4
Jim White Aug 2007 Canberra, Australia 3·11 Posts	Lehmer method in a Nutshell (updated) let Q be the set of primes {2, ... p}, and let n = $\pi$ (p) =\|Q\| be the number of primes in Q let Q* denote the set of p-smooth numbers let Q' be the set of all odd, square-free members of **Q,*** so \|Q'\| = $2^{n-1}$ let S(p) be the set of all integers $a$ for which $a({a-1}) \in$ Q* Algorithm to Enumerate S(p) ("Lehmer's Method") for every $D \in$ Q', find $(x_1, y_1)$ , the fundamental solution to the Pell equation $x^2 - 2Dy^2 = 1$ . If ... $y_1$ ... is p-smooth, then $x_1= 2a - 1$ , and both a and a - 1 are p-smooth, ie $a \in$ S(p) repeat the step above for 4D, ie solving $x^2 - 4Dy^2 = 1$ in both cases, should the Pell equation produce a p-smooth result, then also inspect the multiple-solution sequence $(x_k, y_k)$ , for $k = {2, 3 ... max(3, (p+1)/2)}$ . Check each $y_k$ ... for p-smoothness, and the corresponding $x_k$ will also produce a member of S(p) Remarks 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 ... $y_1$ divides all other ... $y_k$ ... 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: $x_{n+1} = {x_1}{x_n} + Dy_1{y_n}$ $y_{n+1} = {x_1}{y_n} + x_n{y_1}$ $x_{n+1} = 2x_1{x_n} - x_{n-1}$ $y_{n+1} = 2x_1{y_n} - y_{n-1}$

2007-08-31, 10:40	#5
Jim White Aug 2007 Canberra, Australia 3·11 Posts	Størmer's Original Solution Carl Størmer 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 S(p) was in fact always finite. 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 $D$ than the one we used above is shown to produce all the required values as fundamental solutions. Størmer showed that solving $x^2 - Dy^2 = 1$ for all $D \in$ *Q'', where Q''* is the set of all cube-free members of Q* would necessarily produce all $a \in$ S(p) in some solution as $x_1 = 2a-1$ It boils down to showing that every target value will turn up as fundamental solution for any $D$ that correctly matches the target's prime composition, and which also gets the exponent parities right too. The number of equations that need to be checked, though, is much greater than in Lehmer's later refinement, as we have \|*Q''\| = $3^n - 2^n$ Mind you, even with Lehmer's reduction of the problem to $O(2^n)$ , 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 $O(10^n)$ Last fiddled with by Jim White on 2007-08-31 at 10:54*

2007-09-01, 22:13	#8
Jim White Aug 2007 Canberra, Australia 33₁₀ Posts	A Brute-force Method The only alternative method for complete enumeration appears to be a brute-force search of all $S \in$ Q, a search for which we will obviously need a bound,* $S \leq X$ Any results are only provably correct if the bound $X$ 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 Q* recursively, using the bound to limit the depth. For each $S \in$ Q* we simply test $S-1$ for p-smoothness. Crude as it sounds, and also bearing in mind that our target subset S(p) is but a tiny fraction of the bounded subset of Q,* 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 25 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 $O(2^n)$ , 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 $2^{64}$ to $2^{75}$ . 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 p actually divides the solution Attached Thumbnails Last fiddled with by akruppa on 2007-09-05 at 12:46 Reason: fixed some TeX

2007-08-30, 23:58	#2
Jim White Aug 2007 Canberra, Australia 3×11 Posts	Here is an entry from the table above, showing the factors of the corresponding pair: 124225935845233319439173 13^6, 17, 37, 53, 67, 79, 83, 101, 127, 137 124225935845233319439174 2, 3^4, 7, 19^2, 23^2, 29, 31, 47, 97^2, 113^3 There are no greater values possible for 137-smooth pairs of consecutive integers.

2007-09-01, 17:02	#7
Citrix Jun 2003 2×7×113 Posts	If some one is interested in writing the code, may be we can search for more numbers in this series. http://www.mersenneforum.org/showthread.php?t=5630 http://www.research.att.com/~njas/sequences/A116486 Thanks.

2007-09-01, 23:35	#11
Jim White Aug 2007 Canberra, Australia 3×11 Posts	Incremental Enumeration with Lehmer Method The incremental version of the problem is: Given some prime p, for which n do we have gpd(n(n-1)) = p? 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 p). For the Pell equation approach, however, we need to adapt Lehmer's method slightly. We look at $D$ values corresponding to each of the $2^{n-1}$ combinations of primes not including the new $p$ , and for each such $D$ we check the Pell equations for both $pD$ and $p^2D$ , 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 p' <= p , so it is useful to have. Last fiddled with by Jim White on 2007-09-01 at 23:39

Similar Threads
Thread	Thread Starter	Forum	Replies	Last Post
B1 and B2 in P-1 method	Miszka	Math	13	2013-12-27 20:23
New Method	Unregistered	Miscellaneous Math	14	2013-05-24 10:55
Lucas-Lehmer	Dougal	Information & Answers	9	2009-02-06 10:25
Emma Lehmer	R.D. Silverman	Information & Answers	1	2007-08-28 06:20
Størmer's theorem	grandpascorpion	Factoring	0	2006-09-10 04:59