Find all $k$-smooth pairs

Introduction

In number theory, we say a number is $$\k\$$-smooth when its prime factors are all at most $$\k\$$. For example, 2940 is 7-smooth because $$\2940=2^2\cdot3\cdot5\cdot7^2\$$.

Here, we define a $$\k\$$-smooth pair as two consecutive integers which both are $$\k\$$-smooth. An example of 7-smooth pair will be $$\(4374,4375)\$$ because $$\4374=2\cdot3^7\$$ and $$\4375=5^4\cdot7\$$. Fun fact: This is actually the largest 7-smooth pair.

Størmer proved in 1897 that for every $$\k\$$, there are only finitely many $$\k\$$-smooth pairs, and this fact is known as Størmer's Theorem.

Challenge

Your task is to write a program or function that, given a prime number input $$\k\$$, outputs or returns all $$\k\$$-smooth pairs without duplicate (order within the pair does not matter) in any order you want.

Please be noted that for prime numbers $$\p\$$ and $$\q\$$, assuming $$\p, all $$\p\$$-smooth pairs are also $$\q\$$-smooth pairs.

Sample I/O

Input: 2
Output: (1, 2)

Input: 3
Output: (1, 2), (2, 3), (3, 4), (8, 9)

Input: 5
Output: (1, 2), (2, 3), (3, 4), (4, 5), (5, 6), (8, 9), (9, 10), (15, 16), (24, 25), (80, 81)

Input: 7
Output: (1, 2), (2, 3), (3, 4), (4, 5), (5, 6), (6, 7), (7, 8), (8, 9), (9, 10), (14, 15),
(15, 16), (20, 21), (24, 25), (27, 28), (35, 36), (48, 49), (49, 50), (63, 64),
(80, 81), (125, 126), (224, 225), (2400, 2401), (4374, 4375)


Restriction

The program or function should theoretically terminate in finite time for all inputs. Standard loopholes are disallowed by default.

Winning Criteria

As this is a challenge, shortest valid submission for each language wins.

• Could you please add test cases for 2, 3, and 5? – Jonathan Allan May 12 at 18:01
• @JonathanAllan 2-, 3- and 5- smooth pairs are included in the 7-smooth pairs, but I will add the cases for clarity – Shieru Asakoto May 13 at 5:23
• Is having (1, 2) part of the output mandatory?.. – Kevin Cruijssen May 13 at 10:04
• @KevinCruijssen Yes, all outputs should contain the (1, 2) pair. – Shieru Asakoto May 13 at 13:43

JavaScript (ES7),  234  232 bytes

Finds the solutions by solving Pell equations of the form $$\x^2-2qy^2=1\$$, where $$\q\$$ is a $$\P\$$-smooth square free number.

This is an implementation of Derrick Henry Lehmer's procedure, derived from Størmer's original procedure.

Returns an object whose keys and values describe the $$\P\$$-smooth pairs.

P=>[...Array(P**P)].map((_,n)=>(s=(n,i=0,k=2)=>k>P?n<2:n%k?s(n,i,k+1):s(n/k,i,k+i))(n,1)&&(_=>{for(x=1;(y=((++x*x-1)/n)**.5)%1;);(h=(x,y)=>k--&&h(X*x+n*Y*y,X*y+Y*x,x&s(x=~-x/2)&s(x+1)?r[x]=x+1:0))(X=x,Y=y,k=P<5?3:-~P/2)})(),r={})&&r


Try it online!

How?

The helper function $$\s\$$ tests whether a given integer $$\n\$$ is a $$\P\$$-smooth number when it's called with $$\i=0\$$, or a square free1 $$\P\$$-smooth number when it's called with $$\i=1\$$.

s = (
n,
i = 0,
k = 2
) =>
k > P ?
n < 2
:
n % k ?
s(n, i, k + 1)
:
s(n / k, i, k + i)


We look for all square free1 $$\P\$$-smooth numbers in $$\[1..P^P-1]\$$, where $$\P^P\$$ is used as an upper bound for $$\P!\$$.

P=>[...Array(P ** P)].map((_, n) => s(n, 1) && (...))


For each number $$\n\$$ found above, we look for the fundamental solution of the Pell equation $$\x^2-ny^2=1\$$:

(_ => {
for(x = 1; (y = ((++x * x - 1) / n) ** .5) % 1;);
...
})()


(the above code is the non-recursive version of my answer to this other challenge)

Once the fundamental solution $$\(x_1,y_1)\$$ has been found, we compute the solutions $$\(x_k,y_k)\$$ with $$\k\le max(3,(P+1)/2)\$$, using the recurrence relations:

$$x_{k+1}=x_1x_k+ny_1y_k\\ y_{k+1}=x_1y_k+y_1x_k$$

For each $$\x_k\$$, we test whether $$\x_k\$$ is odd and both $$\(x_k-1)/2\$$ and $$\(x_k+1)/2\$$ are $$\P\$$-smooth. If so, we store them in the object $$\r\$$.

( h = (x, y) =>
k-- &&
h(
X * x + n * Y * y,
X * y + Y * x,
x &
s(x = ~-x / 2) &
s(x + 1) ?
r[x] = x + 1
:
0
)
)(X = x, Y = y, k = P < 5 ? 3 : -~P / 2)


1: Because it does not test the primality of the divisors, the function $$\s\$$ will actually be truthy for some non-square free numbers, even when it's called with $$\i=1\$$. The idea is to filter out most of them so that not too many useless Pell equations are solved.

• Hi Arnauld! I just couldn't wrap my head around these two: x = ~-x / 2 and -~P / 2.Are these some kind of rounding... – rv7 May 17 at 17:49
• @rv7 ~x is a bitwise NOT, which computes -(x+1). Therefore, ~-x is -(-x+1) = x-1 and -~x is -(-(x+1)) = x+1. Like all bitwise operations in JS, only the 32-bit integer part is taken into account. So the they can indeed be used for rounding. But both $x$ and $P$ already are integers here. – Arnauld May 17 at 18:08

Jelly, 16 14 bytes

4*ÆfṀ<ɗƇ‘rƝLÐṂ


Try it online!

Checks for pairs up to $$\4^k\$$ which is inefficient for larger $$\k\$$ but should ensure none are missed.

Thanks to @JonathanAllan for saving 1 byte!

Explanation

4*ÆfṀ<ɗƇ‘rƝLÐṂ  | Monadic link, input k

4*              | 4**k, call this n
ɗƇ        | For each number from 1..n filter those where:
Æf            |   - Prime factors
Ṁ           |   - Maximum
<  ‘       |   - Less than k+1
rƝ     | Inclusive range between neighbouring values
LÐṂ  | Keep only those of minimum length (i.e. adjacent values)

• Are you sure that $4^k$ will always be large enough? In my solution, I used $k!^2$ but JonathanAllan was not sure that it would always be large enough. If $4^k$ does always work, I would be curious to hear an explanation. – Comrade SparklePony May 12 at 20:17
• Thanks for the quick response. I was thinking similarly, but more broadly: "factorial gets pretty high quickly, it's probably big enough." (turns out it wasn't unless I squared it). Congrats on the shorter and more efficient golf, you have my upvote. – Comrade SparklePony May 12 at 20:27
• Note (from oeis.org/A002072) "a(n) < 10^n/n except for n=4. (Conjectured, from experimental data.) - M. F. Hasler, Jan 16 2015". I think we have to stick with the weak bound of Lehmer in projecteuclid.org/download/pdf_1/euclid.ijm/1256067456 (theorem 7) unless we can prove otherwise. – Jonathan Allan May 12 at 22:38
• ...there's an open question over on Mathematics SE asking exactly this too! – Jonathan Allan May 12 at 22:45
• @PeterTaylor that’s for the number of pairs, not the maximum number. The problem is knowing a bound on the maximum number of pairs doesn’t let you stop searching – Nick Kennedy May 13 at 7:15

05AB1E, 8 bytes

°Lü‚ʒfà@


Try it online!

Explanation:

°            # 10 ** the input
Lü‚         # list of pairs up to that number
ʒ        # filtered by...
fà      # the greatest prime factor (of either element of the pair)...
@     # is <= the input


Jelly, 123 bytes

¹©Æ½Ø.;µU×_ƭ/;²®_$÷2ị$}ʋ¥⁸;+®Æ½W¤:/µƬṪ€F¹;Ḋ$LḂ$?ṭ@ṫ-ṚZæ.ʋ¥ƒØ.,U¤-ịWµ1ịżU×®W¤Ɗ$æ.0ị$ṭµ³’H»3¤¡
ÆRŒPP€ḟ2ḤÇ€ẎḢ€+€Ø-HÆfṀ€<ẠʋƇ‘


Try it online!

This is a relatively efficient but long Jelly answer that uses the continued fractions method to solve the fundamental solution for the Pell equations for $$\2×\$$ each k-smooth square-free number, finds $$\\max(3, \frac{k+1}{2})\$$ solutions for each and then checks whether $$\\frac{x-1}{2}, \frac{x+1}{2}\$$ are smooth for each solution. This is Lehmer’s method, as described in the question’s Wikipedia link.

A full program that takes a single argument, $$\k\$$ and returns a list of lists of pairs. The code above doesn’t sort the final output, but the TIO link does.

-11 bytes thanks to nimi

q 1=[1]
q n=(:)<*>q.div n$[x|x<-[2..n],mod n x==0]!!0 f k|let r=all(<=k).q=[(n,n+1)|n<-[1..4^k],r n,r(n+1)]  Try it online! • q n calculates a list of all prime factors of n • f k generates a list of $$\k\$$-smooth pairs for a given k by filtering a list of all pairs • You can loop through [2..n] within p and inline it into q. Try it online! – nimi May 13 at 19:48 Jelly, 24 bytes ³!²R‘Ė ÇÆFḢ€€€’<³FȦ$€Tị¢


Try it online!

This takes a long time for 7, but it computes much faster if you remove the squaring of the factorial: Try it online!

Explanation:

³!²R‘Ė                Generates a list like [[1,2],[2,3],...]
³!²                  Take the square of the factorial of the input
R                 Range 1 through through the above number.
‘Ė               Decrement and enumerate, yielding desired list

ÇÆFḢ€€€’<³FȦ$€Tị¢ Ç Get the list of pairs ÆF Get the prime factors of each number Ḣ€€€ Get the base of each ’<³ Is each base less than or equal to the input? FȦ$€       Check that all bases of a pair fit the above.
T      Get a list of the truthy indexes
ị¢    Index into the original list of pairs
Implicit output


-3 bytes thanks to @JonathanAllen

• I don't read Jelly, can you give an explanation on how this works? – Embodiment of Ignorance May 12 at 17:57
• I don't think this works - isn't (8,9) a 3-smooth pair since $8=2^3$ and $9=3^2$? – Jonathan Allan May 12 at 18:00
• I'm not sure it is though. What makes you think that will hold? – Jonathan Allan May 12 at 18:31
• @JonathanAllan Naive optimism and the fact for all examples I've seen (admittedly not many), the largest pair is less than k! (except for 3, which has a small factorial because it is a small number). – Comrade SparklePony May 12 at 18:38
• The upper bound you're using is on the maximum number used in a pair, not on the number of pairs (you can't implement an upper bound on the number of pairs in this way as you wont know when to stop looking!) See theorem 7 for an upper bound on the product of the largest pair. – Jonathan Allan May 12 at 19:34

Python 3 + sympy, 116 bytes

import sympy
def f(k):x=[i for i in range(2,4**k)if max(sympy.factorint(i))<=k];return[(y,y+1)for y in x if y+1in x]


Try it online!

Python 3 + sympy, 111 bytes

from sympy import*
def f(k):
b=1
for i in range(2,4**k):
x=max(factorint(i))<=k
if x&b:print(i-1,i)
b=x


Try it online!

Two variations on my Jelly answer but in Python 3. They both define a function which accepts an argument k. The first returns a list of tuples of the pairs that meet the criteria. The second prints them to stdout.

Wolfram Language (Mathematica), 241 bytes

uses Pell equations

(s=#;v@a_:=Max[#&@@@#&/@FactorInteger@a]<=s;Select[{#-1,#+1}/2&/@(t={};k=y=1;While[k<=Max[3,(s+1)/2],If[IntegerQ[x=Sqrt[1+2y^2#]],t~AppendTo~x;k++];y++];t),v@#&]&/@Join[{1},Select[Range[3,Times@@Prime@Range@PrimePi@s],SquareFreeQ@#&&v@#&]])&


Try it online!

Pyth, 15 bytes

f!>#QsPMTCtBS^4


Try it online!

Uses Nick Kennedy's observation that no output number will be larger than 4^k.

05AB1E, 16 bytes

°LʒfàI>‹}Xšü‚ʒ¥


Try it online (extremely inefficient, so times out for $$\n\gt3\$$..). Here a slightly faster alternative, although still pretty slow..

Explanation:

°                # Take 10 to the power of the (implicit) input
L               # Create a list in the range [1, 10^input]
ʒ              # Filter this list by:
fà            #  Get the maximum prime factor
I>‹         #  And check if it's smaller than or equal to the input
}Xš      # After the filter: prepend 1 again
ü‚    # Create pairs
ʒ   # And filter these pairs by:
¥$$$$ #  Where the forward difference / delta is 1


Stax, 14 bytes

Θ",²aÇu,á‼⌐çLB


Run and debug it

This is not the shortest possible program, but it begins producing output as soon as matching pairs are found. It does terminate eventually, but output is produced as it's found.

Ruby, 89+8 = 97 bytes

Uses the -rprime flag. For each number $$\i\$$ from 1 to $$\4^n\$$, map it to [i, i+1] if both are $$\n\$$-smooth, otherwise map it to false, then prune all false from the list.

->n{g=->x{x.prime_division.all?{|b,_|b<=n}};(1..4**n).map{|i|g[i]&&g[i+1]&&[i,i+1]}-[!1]}


Try it online!