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

Jim White · 2007-09-03, 11:07

Here is an updated version of the timing chart posted above, showing Robert Gerbicz's "Bounded Lehmer" performance on the same system.

Forgive the clumsy plotting, the gradient for his plotted times is in reality straight (uncannily so) and at 45 degrees, ie O(2^n). Here are the actual times, and verification of correct incremental counts (NS) ...

Code:

 > N19, p 71, NS =  790, et =   7 sec., Max = 12820120234376,
 > N20, p 73, NS =  943, et =  15 sec., Max = 119089041053697,
 > N21, p 79, NS = 1201, et =  30 sec., Max = 2286831727304145,
 > N22, p 83, NS = 1401, et =  64 sec., Max = 9591468737351909376,
 > N23, p 89, NS = 1738, et = 132 sec., Max = 17451620110781857,
 > N24, p 97, NS = 1955, et = 274 sec., Max = 166055401586083681,
 > N25, p101, NS = 2240, et = 572 sec., Max = 49956990469100001,
 > N26, p103, NS = 2793, et =1182 sec., Max = 4108258965739505500,

The importance of reducing the number of Pell equations that are fully solved is obvious. As the maximum p value rises, the typical Pell equation determinant expands and so do the period lengths. The longer the period length then the greater must be the fundamental solution magnitudes, so having any reasonable bound on the size of that fundamental solution is going to save substantial wasted effort.

Jim White · 2007-09-03, 14:30

The use of a size bound in Lehmer's method as demonstrated by Robert's code above is very effective, in fact more effective than another bound I have been using.

Period-length is seen to be a potentially effective bound on inspection of the period length of the producers, that is, for those the values of D that produced smooth-pair results.

Here's the table given at the top of this thread with an extra column, L(S) indicate the maximum period length of producers in the corresponding run.

Code:

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

I have attached another chart which is just the one above with a 4th "suit" (spades), which shows the times for my implementation of Lehmer with period length bounding

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

If I introduce period bounding into Robert's program it makes no difference at all - actually it runs slightly slower, unless I specify very accurate bounds.

If I allow both bounds to be checked and count the cases, it is obvious that the size bound check is much more effective, kicking in first in well over 90% of cases (even though it is by contrast a fairly relaxed bound).

Armed with this knowledge, I think we might still get a further shift to the right on the chart by blending the two ...

I'll look into that tomorrow!

R. Gerbicz · 2007-09-03, 14:46

In fact there is a much faster way to search the solution if we assume that say S'(p[N])<2^(3*N+1) is really true.
Unfortunately in practice for example for N=35 ( the largest known value ) this isn't give speedup,
, it is still small N value.
The algorithm:

suppose that we know: a=S'(p[n])<2^(3*n+1), but x=S'(p[n])=2*a-1 and x^2-2*D*y^2=1 or x^2-4*D*y^2=1
From this 2^(6*n+4)>4*a^2>x^2>2*D*y^2>=2*D so D<2^(6*n+3)
Is it possible, that all prime numbers up to p[n] divides D ? For large n this can't be true:
because from prime number's theorem p[n] is about n*log(n) and the product of all primes up to
n*log(n) is e^(n*log(n))=n^n ( also from prime number's theorem ) and for large n this is
larger than the previous bound: 2^(6*n+3)

I haven't calculated how many squarefree numbers are there up to 2^(6*n+3) whose smoothness
is p[n]~n*log(n). But I think this number is much smaller than 2^n, my guess is about
c^(n*log(log(n))/log(n)) ( where c>1 ).
This is similar to Dickman's problem.

By backtracking it is possible to generate only p[n] smooth and squarefree D values and then
solve for only these values the Pell's equations.

Another approach what I've tried yesterday: by bactracking generate all p[n] smooth numbers
up to 2^(3*n+1) that are divisible by p[n] and check if k+1 or k-1 is also smooth or not.
This program was about 100 times slower than my posted program. The above method I think is
much better ( we use only squarefree p[n] smooth numbers but yes there is a larger
computation to solve the pell's equations).

R. Gerbicz · 2007-09-03, 14:54

Quote:

Originally Posted by Jim White

I have attached another chart which is just the one above with a 4th "suit" (spades), which shows the times for my implementation of Lehmer with period length bounding
...

If I introduce period bounding into Robert's program it makes no difference at all - actually it runs slightly slower, unless I specify very accurate bounds.

Probably the reason for that period length bounding and bounding the p[n] values is almost the same is a deep theorem: http://en.wikipedia.org/wiki/Khinchin%27s_constant

Citrix · 2007-09-03, 15:52

R. Gerbicz,
Can your program be modified to search for log smooth numbers? If yes, could you provide an exe file.

Thank you.

R. Gerbicz · 2007-09-03, 22:27

Quote:

Originally Posted by Citrix

R. Gerbicz,
Can your program be modified to search for log smooth numbers? If yes, could you provide an exe file.

Thank you.

For more than 2 consecutive numbers this approach is (much) slower than my posted method in Smooth power trios thread.
And that algorithm already contain some very fast methods for example identify smooth numbers by sieving, what isn't possible when we are using Pell's equations.

Jim White · 2007-09-04, 04:50

Robert, I made a stupid error in the timing chart above, and your time-line is one position to the right of where it should be - I had your program displaying summary lines like this

> N24, p 97, NS = 1955, et = 274 sec., Max = 166055401586083681,

... I plotted that time against p=97, but from the counter (1955) it's obviously the result for p=89

Corrected plot is attached ..

I've figured out the gist of your code - it just took me a while to understand the comments

Thanks for your latest feedback, I'll throw some ideas I have in the ring later today after I've read them properly ...

Agree Pell Eqn solving is the real work - the earlier that we can abort a solution the better ....

Even 100% optimal code is going to run in exponential time, $c^n$ . We have got c down to around 2, and since the number of equations to be considered is $2^n$ that seems unlikely to change much.

To open a door to a seriously larger set of primes seems only likely with some of combinatorial optimisation feature, like some way to say "you need not bother looking at any node below this one".

A better depth bound ... I have interesting data on properties of solution sets, number of prime divisors in a and a+1, these distributions are very compact and illustrate why 99.99% of the time is spent looking for 0.01% of the targets!

Jim White · 2007-09-04, 08:19

Quote:

Originally Posted by R. Gerbicz

Probably the reason for that period length bounding and bounding the p[n] values is almost the same is a deep theorem: http://en.wikipedia.org/wiki/Khinchin%27s_constant

Yes, it could be something like that ... not that I really understand the subject that well

I do have a copy of a 1966 paper co-authored by Daniel Shanks with the rather long title "Questions Concerning Khintchine's Constant and the Efficient Computation of Regular Continued Fractions" ....

Sadly, like so many others, it is still in my pile marked "Papers I Really Should Read". This pile, by the way, exhibits superlinear growth, dwarfing the other piles such as "Papers I have Read", and the even smaller "Papers I have Read and Understood"

I try not to think about the possibility that the size of that last pile might be converging to a limit!

I see on that Wiki page that the CF for e satisfies the relationship, which is interesting because of the geometric mean aspect - I wonder if that might possibly relate to the fact that the arithmetic-geometric-mean (AGM) method happens to be a good way to compute the function log(x) for arbitrary real x > 0?

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

Quote:

Originally Posted by R. Gerbicz

By backtracking it is possible to generate only p[n] smooth and squarefree D values and then
solve for only these values the Pell's equations.

That's more or less what the "Spades" program does, if I understand you correctly. To do the 'incremental' enumeration for 97, say, I generate all square-free q in {3,5,7,11 ... 89}*, then look at the 4 cases

$D \in \lbrace \sp 2 \times 97q,\sp 2 \times 97^2q, \sp4 \times 97q,\sp 4 \times 97^2q\rbrace$

The appearance of spades above is a result of my trying to find a way to insert spacing into a tex string - I was trying "\sp" and got a surprise, in spades! I never did find a way to insert spaces yet, so they will do instead!

I've got two different generators I can use for the square-free combo's - a recursive descent, and a "breadth-first" (ie all combns of 2 primes, then all combns of 3, etc).

Both are suitable for bounding via the number of primes ... I've been using rough but effective values based on what I've seen.

I've been looking at the distributions of numbers of divisors in the target values, rather than the D values themselves, because the target values are a fixed quantuity, whereas D sets can be formed different ways, and there may turn out to be better "producer superset" than the one we have now.

On the other hand, maybe not!

Jim White · 2007-09-04, 09:14

Quote:

Originally Posted by R. Gerbicz

Another approach what I've tried yesterday: by bactracking generate all p[n] smooth numbers
up to 2^(3*n+1) that are divisible by p[n] and check if k+1 or k-1 is also smooth or not.
This program was about 100 times slower than my posted program. The above method I think is
much better ( we use only squarefree p[n] smooth numbers but yes there is a larger
computation to solve the pell's equations).

I called that the "Brute force" or "generative search" above - the plot line marked $\cl$ (aha!) on the chart is for a program which does just that ...

The spikes are because of the very tight known bounds which have those expensive little jumps ...

Bounding on the solution size (or the D size) is very much more effective, and being out by a factor of 10 is unlikely to be of noticable impact (maybe an extra solution step in some subset of all the Pell equations)

Jim White · 2007-09-04, 10:05

From this $2^{6n+4}>4a^2>x^2>2Dy^2 >=2D$ , so we have $D<2^{6n+3}$

Here is an example, the biggest producer (of a solution) for p = 103

$D = 1670718379371535594099740$

D = 2^2, 3, 5, 13, 17, 19, 23, 29, 31, 37, 41, 53, 59, 67, 97, 101, 103
Y = 142
X = 91771952961741

The maximum number of prime divisors for D is here 17 (of the 27 in the full set 2...103)

And the producer of the biggest solution:

$D = 500884459304740$

Y = 1726163680847580
X = 19316158377073923834001

There is only 1 producer value here with 17 divisors, but there are C(25,15) = 4,903,140 ways to select the composition, and it would be nice indeed to be able to find that combination by some other means than checking 4 x 4,903,140 cases ...

Here's the distribution for number of solutions v numbers of prime divisors of the producer D, for p = 103. The divisor counts are inclusive of the two constant divisors {2, 103}:

Code:

npd =  3, solns =  22
npd =  4, solns =  76
npd =  5, solns = 152
npd =  6, solns = 355
npd =  7, solns = 502
npd =  8, solns = 598
npd =  9, solns = 593
npd = 10, solns = 467
npd = 11, solns = 299
npd = 12, solns = 149
npd = 13, solns =  74
npd = 14, solns =  41
npd = 15, solns =   6
npd = 16, solns =   3
npd = 17, solns =   1

This well-balanced pattern appears to persist across all the cases I've looked at so far .... the number of solutions with few divisors is matched by the number with many, but the former are very much easier to find, of course.

The biggest D here is only about 80 bits, not much bigger than the biggest solution a = 19316158377073923834001 which is about 74 bits.

It would be nice, wouldn't it, if the greatest D required was generally not much bigger than the max solution size

I'll try and get more info from what results I have ....

2007-09-03, 11:07	#23
Jim White Aug 2007 Canberra, Australia 3×11 Posts	Timing chart update (Gerbicz's "Bounded Lehmer") Here is an updated version of the timing chart posted above, showing Robert Gerbicz's "Bounded Lehmer" performance on the same system. Forgive the clumsy plotting, the gradient for his plotted times is in reality straight (uncannily so) and at 45 degrees, ie O(2^n). Here are the actual times, and verification of correct incremental counts (NS) ... Code: > N19, p 71, NS = 790, et = 7 sec., Max = 12820120234376, > N20, p 73, NS = 943, et = 15 sec., Max = 119089041053697, > N21, p 79, NS = 1201, et = 30 sec., Max = 2286831727304145, > N22, p 83, NS = 1401, et = 64 sec., Max = 9591468737351909376, > N23, p 89, NS = 1738, et = 132 sec., Max = 17451620110781857, > N24, p 97, NS = 1955, et = 274 sec., Max = 166055401586083681, > N25, p101, NS = 2240, et = 572 sec., Max = 49956990469100001, > N26, p103, NS = 2793, et =1182 sec., Max = 4108258965739505500, The importance of reducing the number of Pell equations that are fully solved is obvious. As the maximum p value rises, the typical Pell equation determinant expands and so do the period lengths. The longer the period length then the greater must be the fundamental solution magnitudes, so having any reasonable bound on the size of that fundamental solution is going to save substantial wasted effort. Attached Thumbnails Last fiddled with by Jim White on 2007-09-03 at 11:09

2007-09-03, 14:30	#24
Jim White Aug 2007 Canberra, Australia 3·11 Posts	Period length as an effective bound The use of a size bound in Lehmer's method as demonstrated by Robert's code above is very effective, in fact more effective than another bound I have been using. Period-length is seen to be a potentially effective bound on inspection of the period length of the producers, that is, for those the values of D that produced smooth-pair results. Here's the table given at the top of this thread with an extra column, L(S) indicate the maximum period length of producers in the corresponding run. Code: N p S'(p) log2(S') L(S) ================================================= 1. 2 2 1 2. 3 9 3.1699 2 3. 5 81 4 4. 7 4375 6 5. 11 9801 8 6. 13 123201 10 7. 17 336141 10 8. 19 11859211 23.4995 12 9. 23 5142501 22.2940 16 10. 29 177182721 27.4077 14 11. 31 1611308700 30.5856 20 12. 37 3463200000 31.6895 18 13. 41 63927525376 35.8957 14 14. 43 421138799640 38.6155 18 15. 47 1109496723126 40.0103 20 16. 53 1453579866025 40.4027 18 17. 59 20628591204481 44.2297 20 18. 61 31887350832897 44.8580 22 19. 67 12820120234376 43.5435 24 20. 71 119089041053697 46.7090 20 21. 73 2286831727304145 51.0223 22 22. 79 9591468737351909376 63.0565 22 23. 83 17451620110781857 53.9542 28 24. 89 166055401586083681 57.2044 26 25. 97 49956990469100001 55.4715 24 26. 101 4108258965739505500 61.8332 28 27. 103 19316158377073923834001 74.0322 28 28. 107 386539843111191225 58.4234 30 29. 109 90550606380841216611 66.2954 30 30. 113 205142063213188103640 67.4752 36 31. 127 53234795127882729825 65.5290 32 32. 131 4114304445616636016032 71.8011 28 33. 137 124225935845233319439174 76.7173 30 34. 139 3482435534325338530940 71.5606 32 35. 149 6418869735252139369210 72.4428 34 I have attached another chart which is just the one above with a 4th "suit" (spades), which shows the times for my implementation of Lehmer with period length bounding 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). If I introduce period bounding into Robert's program it makes no difference at all - actually it runs slightly slower, unless I specify very accurate bounds. If I allow both bounds to be checked and count the cases, it is obvious that the size bound check is much more effective, kicking in first in well over 90% of cases (even though it is by contrast a fairly relaxed bound). Armed with this knowledge, I think we might still get a further shift to the right on the chart by blending the two ... I'll look into that tomorrow! Attached Thumbnails

2007-09-04, 04:50	#29
Jim White Aug 2007 Canberra, Australia 41₈ Posts	Robert, I made a stupid error in the timing chart above, and your time-line is one position to the right of where it should be - I had your program displaying summary lines like this > N24, p 97, NS = 1955, et = 274 sec., Max = 166055401586083681, ... I plotted that time against p=97, but from the counter (1955) it's obviously the result for p=89 Corrected plot is attached .. I've figured out the gist of your code - it just took me a while to understand the comments Thanks for your latest feedback, I'll throw some ideas I have in the ring later today after I've read them properly ... Agree Pell Eqn solving is the real work - the earlier that we can abort a solution the better .... Even 100% optimal code is going to run in exponential time, $c^n$ . We have got c down to around 2, and since the number of equations to be considered is $2^n$ that seems unlikely to change much. To open a door to a seriously larger set of primes seems only likely with some of combinatorial optimisation feature, like some way to say "you need not bother looking at any node below this one". A better depth bound ... I have interesting data on properties of solution sets, number of prime divisors in a and a+1, these distributions are very compact and illustrate why 99.99% of the time is spent looking for 0.01% of the targets! Attached Thumbnails

2007-09-04, 10:05	#33
Jim White Aug 2007 Canberra, Australia 41₈ Posts	From this $2^{6n+4}>4a^2>x^2>2Dy^2 >=2D$ , so we have $D<2^{6n+3}$ Here is an example, the biggest producer (of a solution) for p = 103 $D = 1670718379371535594099740$ D = 2^2, 3, 5, 13, 17, 19, 23, 29, 31, 37, 41, 53, 59, 67, 97, 101, 103 Y = 142 X = 91771952961741 The maximum number of prime divisors for D is here 17 (of the 27 in the full set 2...103) And the producer of the biggest solution: $D = 500884459304740$ Y = 1726163680847580 X = 19316158377073923834001 There is only 1 producer value here with 17 divisors, but there are C(25,15) = 4,903,140 ways to select the composition, and it would be nice indeed to be able to find that combination by some other means than checking 4 x 4,903,140 cases ... Here's the distribution for number of solutions v numbers of prime divisors of the producer D, for p = 103. The divisor counts are inclusive of the two constant divisors {2, 103}: Code: npd = 3, solns = 22 npd = 4, solns = 76 npd = 5, solns = 152 npd = 6, solns = 355 npd = 7, solns = 502 npd = 8, solns = 598 npd = 9, solns = 593 npd = 10, solns = 467 npd = 11, solns = 299 npd = 12, solns = 149 npd = 13, solns = 74 npd = 14, solns = 41 npd = 15, solns = 6 npd = 16, solns = 3 npd = 17, solns = 1 This well-balanced pattern appears to persist across all the cases I've looked at so far .... the number of solutions with few divisors is matched by the number with many, but the former are very much easier to find, of course. The biggest D here is only about 80 bits, not much bigger than the biggest solution a = 19316158377073923834001 which is about 74 bits. It would be nice, wouldn't it, if the greatest D required was generally not much bigger than the max solution size I'll try and get more info from what results I have .... Last fiddled with by Jim White on 2007-09-04 at 10:15

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

2007-09-03, 14:46	#25
R. Gerbicz "Robert Gerbicz" Oct 2005 Hungary 2²×7×53 Posts	In fact there is a much faster way to search the solution if we assume that say S'(p[N])<2^(3N+1) is really true. Unfortunately in practice for example for N=35 ( the largest known value ) this isn't give speedup, , it is still small N value. The algorithm: suppose that we know: a=S'(p[n])<2^(3n+1), but x=S'(p[n])=2a-1 and x^2-2Dy^2=1 or x^2-4Dy^2=1 From this 2^(6n+4)>4a^2>x^2>2Dy^2>=2D so D<2^(6n+3) Is it possible, that all prime numbers up to p[n] divides D ? For large n this can't be true: because from prime number's theorem p[n] is about nlog(n) and the product of all primes up to nlog(n) is e^(nlog(n))=n^n ( also from prime number's theorem ) and for large n this is larger than the previous bound: 2^(6n+3) I haven't calculated how many squarefree numbers are there up to 2^(6n+3) whose smoothness is p[n]~nlog(n). But I think this number is much smaller than 2^n, my guess is about c^(nlog(log(n))/log(n)) ( where c>1 ). This is similar to Dickman's problem. By backtracking it is possible to generate only p[n] smooth and squarefree D values and then solve for only these values the Pell's equations. Another approach what I've tried yesterday: by bactracking generate all p[n] smooth numbers up to 2^(3*n+1) that are divisible by p[n] and check if k+1 or k-1 is also smooth or not. This program was about 100 times slower than my posted program. The above method I think is much better ( we use only squarefree p[n] smooth numbers but yes there is a larger computation to solve the pell's equations).

2007-09-03, 15:52	#27
Citrix Jun 2003 2·7·113 Posts	R. Gerbicz, Can your program be modified to search for log smooth numbers? If yes, could you provide an exe file. Thank you.