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

[R. Gerbicz]

Originally Posted by Jim White

It runs fastest of all so far, but I think that's mainly because it has a faster Pell solver - it uses various 'tricks' to speed up the Pell equation solving process (eg going directly from mid-period to fundamental solution).

Thanks.
Yesterday I've figured out some tricks. I'll try to rewrite my program.

Do you know that period-length can be odd? All of the posted period lengths in the maximal solutions are even but there are examples when the period length of the continued fraction expansion of sqrt(n) is odd ( where n is even ).

[Jim White]

Yes. Ignoring the case period = 1, odd periods are particularly rare when you have even D's (as we do currently, under the Lehmer scheme).

Here I've reported over various runs, the maximum odd period length that turned up, and the maximum D which had an odd period:

p= 83  <none>
 
p= 89
    Odd periods     = 37
    Max length   pl = 31, D = 22952210
    Max D        pl = 27, D = 416392730
 
p= 97
    Odd periods     = 60
    Max length   pl = 31, D = 1715930
    MAx D        pl = 17, D = 3774000501602
 
p=101
    Odd periods     = 72
    Max length   pl = 31, D = 20402
    Max D        pl = 29, D = 23010564730570
 
p=103  <none>
 
p=107  <none>

An odd period can generally be treated as just the mid-point of an even period, I think.

[Jim White]

Regarding Pell equation (CF development, etc) efficiency, the restriction of the CF development to the mid-point is definitely worth implementing, even in our relatively short periods - it's still effective for period lengths of say, 8 or more.

This is generally known as "the mid-point trick", it seems

Another trick is one that seems peculiar to this problem - if we are going to reject period lengths that exceed some max L, and we know that the great majority of cases will be thus rejected, then it makes very good sense to separate the CF expansion into two stages:

1. Periodicity + CF Denominators - compute and store just the CF denominators, aborting any cases where periodicity is not detected in less than L steps

2. Convergents - for acceptable cases only, generate the required convergents (via the recurrence relation multipliers given by the stored CF denominators)

Combine with mid-point shortcut in stage 2, and you have program in a nutshell. That's why it's quite fast, even though some aspects are not fully optimised.

We can also use size bounding, of course, during stage 2. However, it looks like what happened when I applying period length bounds to your code, they didn't make any real difference, as the size bound usually triggered first.

With the period-bounding method described here, it's just the opposite - by rejecting periods > L,we have already filtered out most of those that would run over size (and we have not wasted much computational effort on them at all - this is the main efficiency gain)

So we only develop the convergents for those guys we know that have short periods, and we can still keep an eye on the size, although we expect to find such oversize cases less often.

[Jim White]

Originally Posted by Jim White

We can also use size bounding, of course, during stage 2.

However, it looks like what happened when I applying period length bounds to your code, they didn't make any real difference, as the size bound usually triggered first.

With the period-bounding method described here, it's just the opposite - ... ...

Egad! Must have fallen asleep at the keyboard (again!)

I'm sorry! That should read:

Originally Posted by Jim White

We can also use size bounding, of course, during stage 2.

If I add period length bound checks to Robert's original program, however, they don't make any real difference - his existing size bound check is almost always more effective (catches exceptions earlier) than the priod length bound.

With the period-bounding method described here, it's just the opposite -

[Jim White]

The set of D values that we have to generate can be done in various combinatorial ways. Two examples:

1. Recursive Generation - eg: 3 -> 3,5 -> 3,5,7 -> 3,5,7,11 -> 3,5,11 -> 3,7 - > 3,7,11 -> 3,11 ==> 5 -> ... (a pre-order traversal of the combination tree)

2. Generating Same-size Selections - this involves generating all combinations of size k for k = 1, 2, ... . So pass 1 = {3 -> 5 -> 7 -> 11}, pass 2 = {3,5 -> 3,7 -> 3,11 -> 5,7 -> 5, 11 -> 7, 11} and so on

The first method is very much simpler to code, and neither one makes any significant difference to running time in this problem.

But the second method has at least one interesting feature - it can locate the majority of the targets in a flash!

For example, with n = 30 (p = 113), we are combining the 28 primes {3, 5, ... 109} to form D' (which we then combine with 2 and 113 in the various appropriate ways), of which there are possible combinations (well, give or take 1).

Here's a log from a sample run showing the results for various k :

 k             C(n,k)  Solns  Total    run time
 1                28,    19,    23,      
 2               378,    72,    95,      
 3 **           3276,   191,   286,      
 4 *****       20475,   415,   701,      .1s 
 5 ******      98280,   676,  1377,      .4s
 6 *********  376740,   878,  2255,     1.7s
 7 ********* 1184040,   934,  3189,     5.6s
 8 ********* 3108105,   819,  4008,    15.9s
 9 *****     6906900,   575,  4583,    39.7s
10 ****     13123110,   357,  4940,    92.7s
11 **       21474180,   193,  5133,   181.8s
12 *        30421755,    91,  5224,   311.4s
13          37442160,    36,  5260,   471.8s
14          40116600,    18,  5278,   643.8s
15          37442160,     3,  5281,   804.5s
16          30421755,     2,  5283,   934.9s
17          21474180,     1,  5284,  1026.8s
18          13123110,     0,  5284,  1083.0s
19           6906900,     0,  5284,  1112.7s
..... 

So it took just 16 seconds to find 4000 of the 5300 odd solutions, then nearly 20 minutes to find the the rest ...

The number of possible combinations that have to be considered is listed, and this gives you yet another view of the "fundamental" problem ...

The "peak" solution pass seems to be around n/4 for all cases I've looked at so far.

If you write a generic k from n selection generator, then you can do the selection passes in any order you like, of course ....

[Jim White]

These settings might be useful if you want to save time on a particular run:

n      p    D'2   S'2  v(D') v(S')  L(D)
25     97    64    57   13    17     26
26    101    66    63   14    17     28
27    103    72    76   15    18     28
28    107    74    60   15    19     28
29    109    77    68   16    20     32
30    113    87    69   17    20     36
31    127    81    67   16    20     32
32    131    88    73   16    21     30
33    137    88    78   17    21     30

n = index(p)

D'2 = max reduced producer size (bits). The reduced determinants D' are the odd square-free values. ie. subsets of {3, 5, ... p'}, where p' is the n-1'th prime. This bound can save time by helping to eliminate many of the possible divisor combinations that are supposed to be tested.

S'2 = max solution size (bits). If S' is the max S for which S(S-1) is a hit, then S'2 = . Like knowing max period length, this allows Pell equation solution development to be abandoned early when known to be futile.

v(D') = v(N) is the number of prime divisors of N. v(D') is the highest value for any (reduced) producer D'. In other words, it's pointless combining any more than this when making any D'.

v(S') = the highest number of divisors of any S(S-1) pair. Includes all divisors, so includes 2 and p.

L(D) = the largest period length of any producer.

D'2 and L(D) are very effective bounds (when known!)


As an example, a current program I have takes over 1 hour to run (70 mins) for the case n=32, p=131, using a reasonably good period bound (these seem very predictable, eg n + 10 would usually do the job).

Once I have a sharper value for the period bound, and also know the size of the largest producer I need to generate, re-running the job is quicker, completing in 40 minutes.

[liangbch]

I am also interested in this topic.
I searched all p34-smooth number pair by my c language program and get 63428 solutions. My arithmetic is base on a conjecture, let (n,n+1) is a p34 –smooth number pair, the max of n should at most with 128 bits. The following is the largest 16 solutions (n) that satisfy with both n and n+1 is 139-smooth number

107109384837019296875
118991993751522079210
152295745769656587384
205142063213188103639
507882950033870979071
1043073004436787852800
1338614896790283750000
1492860563840742187104
1805692829069444659200
1811103583756769555714
2187239579855392547199
2888494687736348753790
3482435534325338530939
4114304445616636016031
19316158377073923834000
124225935845233319439173

[liangbch]

I google number 124225935845233319439173 in internet and found a table as bellow, of course, it is spectacular result. The table give the max p_n smooth nubmer pair and n up to 40. I don't understand how to search such result in very very big range.
I believe your state is reasonable. “There is a bigger 157-limit comma”
This table can be found from www.primefan.ru/stuff/math/maxs.xls

Index p k=1 (A002072)
1 2 1
2 3 8
3 5 80
4 7 4374
5 11 9800
6 13 123200
7 17 336140
8 19 11859210
9 23 11859210
10 29 177182720
11 31 1611308699
12 37 3463199999
13 41 63927525375
14 43 421138799639
15 47 1109496723125
16 53 1453579866024
17 59 20628591204480
18 61 31887350832896
19 67 31887350832896
20 71 119089041053696
21 73 2286831727304144
22 79 9591468737351909375
23 83 9591468737351909375
24 89 9591468737351909375
25 97 9591468737351909375
26 101 9591468737351909375
27 103 19316158377073923834000
28 107 19316158377073923834000
29 109 19316158377073923834000
30 113 19316158377073923834000
31 127 19316158377073923834000
32 131 19316158377073923834000
33 137 124225935845233319439173
34 139 124225935845233319439173
35 149 124225935845233319439173
36 151 124225935845233319439173
37 157 6318268857746831540296874
38 163 6318268857746831540296874
39 167 6318268857746831540296874
40 173 6318268857746831540296874

[LaurV]

You (or the guy behind of primefan.ru) should report that here.
My feeling is that the list would not be very difficult to extend, I will give it a try in pari if I find some free time over teh weekend.

[liangbch]

Basically, total 4 OEIS is related with this topic. They are A117581, A002072, A002071 and A145606.
let both m and m+1 is P(n)-smooth number, then A002072= max(m) for each of P(n) smooth number pair, and A117581=max(m+1), so A117581 [i] = A002072[i]+1.
The A002071 is number count for all of m which satisfy m and m+1 is P(n) smooth number.
The A145606 denote that largest number x such that x and x+1 are prime(n)-smooth but not prime(n-1)-smooth.
It is said in both A14506 and A002072 "Jim White, Results to P = 173". If the OEIS owner thinks Jim White's result and computation method is exactly correct, then all of these 4 serial should extend to 40 term. But I think that because Jim White's computation method base on a conjecture, the upperbound should less than a const and this upperbound can not been prove by mathematical theorem, so only a part of serial be given.

Jim White, If I am wrong, please correct me.

[science_man_88]

Originally Posted by Jim White

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):

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

brute force is basically:

gpd(y) = factor(y)[#factor(y)[,1],1]

and a infinite loop of values to check, a less brute force is knowing that for all all primes q<p , all values q*x for x<p are in the set also all powers q^y for q<p are in the set which lead to an infinite set.