All pairs of consecutive 47-smooth numbers, calculated by my Python code for Lehmer's improvement to Størmer's theorem. That is, this file lists pairs x, x+1 such that all prime factors of both x and x+1 are at most 47. There are 1502 such pairs, the largest of which is formed by 1109496723125 = 54 × 7 × 17 × 192 × 312 × 43 and 1109496723126 = 2 × 3 × 112 × 23 × 292 × 412 × 47.
It took nearly two days on a dual-core Pentium, so the next prime 53 may be in reach but would likely take more computer time than I'm willing to devote to it. Judging by A002072, faster codes are out there.
There was a gmp-d library to use the multiprecision GMP with D, but it seems vanished. But your program look simple enough, so maybe you can translate it to C that calls the GMP (I use GMP with Python + Psyco, but if you have many small numbers it's not a significative improvements on the python built-in longs). ShedSkin doesn't support GMP yet, but maybe it's not difficult to use it.
I did make one change to my code since running the computation that should speed things up significantly: if the first solution to a given Pell equation is already not smooth, the rest of them for the same equation won't be either, so we can stop early. As there are 32768 Pell equations (for the 47-smooth calculation) but only 1502 smooth pairs, this early escape should happen very frequently.
For the record, Lehmer did make this point clear - in Section 6 (Implementation notes) he says "since every y[n] is divisible by y, it is useless to examine multiple solutions if y ... [is not smooth]"
Fixing this will not necessarily translate into dramatic speed increases in your program - the additional imposed cost of iterating the multiple solutions, assuming you check no more than M as specified, is likely to be modest compared wih the cost of computing the Pell solution.
For maximum speed clearly GMP (suitably CPU-targeted) is going to give the best results on most common CPU's.
Here are my times recently recorded for P = 47 and similar cases. The test rig is an AMD64 2.2 GHz.
Columns are no of primes, highest P, number of eqns, NS = number of pairs found, job time (seconds).
13 p = 41 8K NS = 869, .870 s 14 p = 43 16K NS = 1153, 7.490 s 15 p = 47 32K NS = 1502, 78.910 s 16 p = 53 64K NS = 1930, 1163.150 s 17 P = 59 128K NS = 2454, 15153~ 4hrs 18 P = 61 256K NS = 3106, 144838~ 40hrs
The first couple of runs look good, but the growth in cost of 10x per additional prime is the dominant, and crippling feature ... and that value "10" may also tutn out to be growing ...
Each extra prime doubles the number of equations, but also raises the average "weight" of the D values, so numbers and periods get heavier and heavier....
And, as you may have already realised, you can't handle each additional prime "incrementally" - if I have all 43-smooth pairs, and want to get the additional ones for 47, I still have to redo all the previous equations (with D 43-smooth) as some of them will produce 47-smooth pairs
So while an incremental approach (forcing 47 as a divisor of D) will generally identify _most_ new solutions, it will not be complete.
I don't totally understand why A002071 only lists counts for up to P = 47, yet A002072 lists maximum pair values for up to P = 97.
I sent an enquiry to eppstein (as A002071 suggests), but have had no reply ... from my results, one wonders just how long the A002071 results took to produce, and just how it was done
I'll try another enquiry, this time to Fred Schneider ... see if he has some tricks up his sleeve?
Erratum: "one wonders just how long the A002072 results took ...."
I don't totally understand why A002071 only lists counts for up to P = 47, yet A002072 lists maximum pair values for up to P = 97. I sent an enquiry to eppstein (as A002071 suggests), but have had no reply
That would be me. I've been traveling and largely away from my email. Though in the best of times I've been known not to reply to my emails...
Anyway, A002072 links to some code by Don Reble that may have been used to compute it — if so, I don't know why it would be so much faster than the code I used to compute the additional terms in A002071.
I clearly failed to notice the significance of your user id "0xDE", it does not signify Diophantine (or DIfferential) Equations! :-)
I thought perhaps my email enquiry might have been skimmed by a spam-detector (understandable in these times) ...
My recent progress: I have computed the complete solution sets up to n=13 (p <=97), and the maximum values I encountered do tally with those givem in A002702 (which, by the way, I would prefer to see give the explicit maximal value for each prime, rather than repeating values)
The program I've developed took only 2 mins to do this, and involves a growth factor of just 2+e rather than 10+e for each additional prime.
The 40 hours it cost to produce the results for n=18 (p<=61) should now buy me a solution set for n=31 (p<=127), and I have just started a run with that in mind.
My values for A002071 for n=13 to 25 are:
13 p = 41 NS = 869 14 p = 43 NS = 1153 15 p = 47 NS = 1502 16 p = 53 NS = 1930 17 P = 59 NS = 2454 18 P = 61 NS = 3106 19 P = 67 NS = 3896 20 P = 71 NS = 4839 21 P = 73 NS = 6040 22 P = 79 NS = 7441 23 P = 83 NS = 9179 24 P = 89 NS = 11134 25 P = 97 NS = 13374
a) Are you able to confirm any of these?
b) Are you aware of any known values beyond n=25, p=97? One of the problems with my method is that (currently) it is not yet provably correct - I need to sharpen Lehmer's "Theorem" before I can claim the p=127 results will be definitive. Meantime all I can say is that they will almost certainly be correct!
This work is not purely for fun - my current research is motivated by a broader practical application, in which one approach to the problem is dependent on having an "efficient" method of finding these "smooth-pairs".
(MSI, ANU, Canberra)
The first three terms you already know from A002071. I'm not really in a position to do anything compute-intensive while traveling. But given the 10x growth in my code's runtime it seems that it might take several days to compute even the next term after that.
No need to waste any time on it, then ...
I'll send those figures over to Mr Schneider (last indicated correspondent at A002072) for comment ...
Regarding the "express method", this is not a new method, just a nice refinement of Lehmer's method.
I only "discovered" the trick a few days ago, hence the uncertainty as to its verifiability ... I've been preoccupied with its implementation and performance ... I am devoting this weekend, however, to a look at whether it can be shown to be correct, since ...
I'm giving an informal talk (as PhD candidates are obliged to do from time to time) on Monday, and I now wish to abandon my previous topic and talk about this fascinating problem ...
so I'll be knocking up some notes, which I should be able to provide shortly after...
The job for n=31, p <= 127, took just 3h 20m.
Predicted A00271 terms:
P = 97 NS = 13374 P = 101 NS = 16167 P = 103 NS = 19505 P = 107 NS = 23361 P = 109 NS = 27939 P = 113 NS = 33208 P = 127 NS = 39240
When setting up that job, I missed a "section", in which 43 additional pairs were to be found (a sparsely populated section it was, too, requiring 2 further hours to check).
P = 97 NS = 13374, max = 49956990469100001 P = 101 NS = 16167, max = 4108258965739505500 P = 103 NS = 19508, max = 19316158377073923834001 P = 107 NS = 23372, max = 386539843111191225 P = 109 NS = 27956, max = 90550606380841216611 P = 113 NS = 33250, max = 205142063213188103640 P = 127 NS = 39312, max = 53234795127882729825
max for each P is the value s+1 in the highest pair (s, s+1), with P dividing either s or s+1.