# An algorithm to find even sublime numbers

A perfect number is a positive integer whose sum of divisors (except itself) is equal to itself. E.g. 6 (1 + 2 + 3 = 6) and 28 (1 + 2 + 4 + 7 + 14 = 28) are perfect.

A sublime number (OEIS A081357) is a positive integer whose count and sum of divisors (including itself) are both perfect. E.g. 12 is a sublime number because:

• the divisors of 12 are 1, 2, 3, 4, 6, 12
• the count of divisors is 6 (perfect)
• the sum of divisors is 28 (perfect)

The next smallest known sublime number is

6086555670238378989670371734243169622657830773351885970528324860512791691264


and these two numbers are the only known sublime numbers as of 2022. The necessary and sufficient conditions for even sublime numbers have been found in this paper (pdf), but it remains unknown whether odd sublime numbers exist.

The paper outlines the "algorithm" to find even sublime numbers:

Suppose $$\p\$$ is a prime with the following properties,

• $$\q = 2^p − 1\$$ is a prime.
• $$\2^q − 1\$$ is a prime.
• $$\q − 1\$$ can be partitioned into distinct primes $$\l_1, \cdots , l_{p−1}\$$ such that $$\m_i = 2^{l_i} − 1\$$ are also prime for all $$\i\$$.

Then $$\n = 2^{q−1} ( \prod{m_i} )\$$ is a sublime number.

This heavily relies on Mersenne primes, which are primes in the form of $$\2^n-1\$$. Let's define Mersenne exponents be the values of $$\n\$$ such that $$\2^n-1\$$ is prime. The list of known Mersenne exponents can be found on this Wikipedia page; it starts with

2, 3, 5, 7, 13, 17, 19, 31, 61, 89, 107, 127, 521, 607, ...


and only around 50 such exponents are known.

We can rewrite the algorithm in terms of a known set $$\L\$$ of Mersenne exponents:

Loop over $$\p \in L\$$:

• Check that $$\q = 2^p - 1 \in L\$$. If so,
• For each subset $$\\{ l_i \}\$$ of $$\L\$$ of size $$\p-1\$$:
• If it sums to $$\q-1\$$,
• output $$\n = 2^{q-1} ( \prod{(2^{l_i}-1)} )\$$ as a sublime number found.

Since the list of actual Mersenne exponents is fairly limited, we wouldn't get any more sublime numbers by using this algorithm on it. So let's run it on an arbitrary set of positive integers instead.

## Challenge

Given a set of positive integers $$\L\$$, evaluate all possible outputs (values of $$\n\$$) the algorithm above would give you. Since $$\p=1\$$ is kind of an edge case, you may assume that $$\L\$$ does not contain 1.

You may assume that $$\L\$$ is nonempty. If you take $$\L\$$ as a list (as opposed to an unordered set), you may assume that it is sorted. You may output the values of $$\n\$$ in any order, and may output the same value multiple times. If multiple subsets satisfy the condition for a single $$\p\$$, all possible resulting $$\n\$$s must be output.

Standard rules apply. The shortest code in bytes wins.

## Examples

If $$\L=\{2,3\}\$$, $$\p=2\$$ satisfies the condition:

• 2 is in $$\L\$$ and $$\2^2-1 = 3\$$ is also in $$\L\$$.
• A subset of size $$\p-1 = 1\$$ with sum $$\q-1 = 2\$$ exists: $$\\{2\}\$$.
• The output is $$\n = 2^{3-1} (2^2-1) = 12\$$.

For $$\L = \{2, 3, 4, 7\}\$$, $$\p=3\$$ also does:

• 3 is in $$\L\$$ and $$\2^3-1 = 7\$$ is also in $$\L\$$.
• A subset of size $$\p-1 = 2\$$ with sum $$\q-1 = 6\$$ exists: $$\\{2,4\}\$$. ($$\\{3,3\}\$$ does not work since it is just $$\\{3\}\$$.)
• The output is $$\n = 2^{7-1} (2^2-1)(2^4-1) = 2880\$$.

## Test cases

L => [n]
 => []
[2, 3] => 
[2, 3, 4] => 
[2, 3, 4, 7] => [12, 2880]
[3, 4, 7] => []
[3, 4, 5, 6, 7, 15] =>
[218480640 = 2^14(2^3-1)(2^4-1)(2^7-1),
223985664 = 2^14(2^3-1)(2^5-1)(2^6-1)]
[2, 3, 5, 7, 15] =>  (suggested by Jonathan Allan)
(actual Mersenne exponents, truncated and 2 removed)
[3, 5, 7, 13, 17, 19, 31, 61, 89, 107, 127, 521, 607] =>
[6086555670238378989670371734243169622657830773351885970528324860512791691264
= 2^126(2^3-1)(2^5-1)(2^7-1)(2^19-1)(2^31-1)(2^61-1)]

• A challenge about primes without any prime computation? That's do X without Y done right! :p Aug 23 at 14:03
• If one considers all subsets of $L$ and checks only whether the subset's sum and $q$ are in $L$ then all of the given test cases pass. I suggest adding [2, 3, 5, 7, 15] (which should produce  rather than [12, 193511424]). (All six answers posted so far pass this test case.) Aug 23 at 16:17

# PARI/GP, 108 bytes

l->forsubset(#l,s,setminus([p=#s+1,q=2^p-1],l)||sum(i=1,#s,l[s[i]])-q--||print(prod(i=1,#s,2^l[s[i]]-1)<<q))

Attempt This Online!

Takes input as sorted lists.

# JavaScript (V8),  119  117 bytes

Given a list of integers, prints the pseudo-sublime numbers.

L=>L.map(g=(p,_,a,s=L.includes(q=2**p-1)*--q,o=2**q)=>--p?a.map((x,i)=>g(p,_,a.slice(i+1),s-x,o*2**x-o)):s||print(o))


Try it online!

### Commented

L =>                    // L[] = input list
L.map(g =               // for each element in L[],
(                     // invoke the recursive callback function g:
p,                  //   p = current element from L[]
_,                  //   (ignored index of this element)
a,                  //   a[] = input list (truncated on subsequent calls)
s = L.includes(     //   test whether L[] includes:
q = 2 ** p - 1    //     q = 2 ** p - 1
) * --q,            //   decrement q and save either 0 or q in s,
//   which is our target sum
o = 2 ** q          //   o = output product, initialized to 2 ** q
) =>                  //
--p ?                 //   decrement p; if it's not 0:
a.map((x, i) =>     //     for each element x at position i in a[]:
g(                //       do a recursive call:
p,              //         new value of p
_,              //         (ignored index)
a.slice(i + 1), //         remove the i + 1 first elements of a[]
s - x,          //         subtract x from s
o * 2 ** x - o  //         multiply o by 2 ** x - 1
)                 //       end of recursive call
)                   //     end of inner map()
:                     //   else (p = 0):
s || print(o)       //     print o if s = 0
)                       // end of outer map()


# Jelly, 27 bytes

L‘2*’ṭƲf’iS{Ṗ
ŒPçƇµ;S2*_ṠP)


A monadic Link that accepts a list of distinct integers greater than $$\1\$$ and yields a list of the implied pseudo-sublime numbers.

Try it online! Or see the test-suite.

### How?

L‘2*’ṭƲf’iS{Ṗ - Helper Link, valid subset?: Subset; L
L             - length of Subset
‘            - increment -> p'
Ʋ       - last four links as a monad - f(p'):
2           -   two
*          -   exponentiate (p')
’         -   decrement -> q'
ṭ        -   tack -> [p', q']
f      - filter keep only those of [p', q'] which are in L
We now have one of: []; [p']; [q']; or [p'=p, q'=q]
’     - decrement -> F (e.g. F=[p-1, q-1])
{  - using left argument, Subset:
S   -   sum
i    - first 1-indexed index of that sum in F or 0 if not found
Ṗ - pop (implicit range [1..that])
This yields:
[] (falsey) when i resulted in 0 or 1; or
 (truthy) when i resulted in 2 (i.e. p'=p, q'=q, and q-1=sum)

ŒPçƇµ;S2*_ṠP) - Link, get pseudo-sublimes: set of integers > 1, L
ŒP            - powerset of L -> all subsets we could make
Ƈ          - filter-keep those for which:
ç           -   call the Helper Link - f(subset, L)
µ       ) - for each:
e.g. [a, b, ..., n]
S       -   sum (Subset)    a + b + ... + n = q-1
;        -   concatenate  -> [    a,     b, ...,     n,    q-1 ]
2      -   two
*     -   exponentiate -> [  2^a,   2^b, ...,   2^n, 2^(q-1)]
Ṡ   -   sign (Subset)-> [    1,     1, ...,     1]
_    -   subtract     -> [2^a-1, 2^b-1, ..., 2^n-1, 2^(q-1)]
P  -   product      -> pseudo-sublime number


# Vyxal, 34 24 bytes

ṗ'L›:E‹"?↔‹n∑ḟɾ;ƛ∑JEn±-Π


Try it Online!

Mess. Well, still a mess, but -10 bytes by porting Jonathan Allan's Jelly answer, upvote that!

# Jelly, 37 36 bytes (thanks @Steffan!)

2*’
§=⁹aÇP€×2*⁹¤ʋ
Ñe³ƊƇµ’³œcç"Ñ’$F¹Ƈ  Try it online! Also a mess. Link 1: computes $$\2^x-1\$$ 2 Two * To the power of the input ’ Decrement  Link 2: given a list of subsets of length $$\p-1\$$ of the input list, and $$\q-1\$$ (as defined in the question), return a list with all valid answers or 0. §=⁹aÇP€×2*⁹¤ʋ § Compute the sum of each sublist (vectorizes at depth 1) = Check if equal to... (vectorizes) ⁹ ...q-1 a Logical AND with... ʋ ...the four previous links as a dyad Ç The previous link as a monad (2^x-1 for each x in each subset) € For each subset: P Compute the product × Multiply by... ¤ ...a nilad followed by some links, as a nilad: 2 Two * To the power of... ⁹ ...q-1  Link 3 (main link): given the input list, outputs all valid answers. Ñe³ƊƇµ’³œcç"Ñ’$F¹Ƈ
Ƈ               Filter the list using:
Ɗ                    The last three links as a monad
Ñ                           The next link as a monad (wraps around to link 1)
e                          Check if it's present in...
³                         ...the input list
µ              Start a new monadic chain
(At this point, the list contains all 'p' in the input list such that 2^p-1 is also in the input list.)
’             Decrement
œc          All combinations of ... of any length in the list (vectorizes over the right argument)
³            The input list
"        Zipping over the left and right arguments:
ç             Call the previous link as a dyad, with the next link as the right argument
\$     The last two links as a monad:
Ñ           The next link (wraps around to link 1) on each valid 'p', thus generating all the candidates for 'q'
’          Decrement (thus generating all the candidates for q-1)
F    Flatten
Ƈ  Filter the list using:
¹   The truth value of each element itself (thus removing all falsy values from the list)

• S€ is just §. Aug 25 at 19:16

# Charcoal, 46 bytes

ＩＥΦＥＸ²ＬθΦθ＆ιＸ²μ∧⁼Σι⊗⊖Ｘ²Ｌι∧№θ⊕Ｌι№θ⊕Σι×Ｘ²ΣιΠ⊖Ｘ²ι


Try it online! Link is to verbose version of code. Explanation:

ＩＥΦＥＸ²ＬθΦθ＆ιＸ²μ


For each subset of the input list, such that...

∧⁼Σι⊗⊖Ｘ²Ｌι


... the sum of each subset is two less than twice two to the power of its size, and...

∧№θ⊕Ｌι


... the size of the subset is one less than one of the integers in the input list, and...

№θ⊕Σι


... the sum of the subset is also one less than one of the integers in the input list, ...

×Ｘ²ΣιΠ⊖Ｘ²ι


... raise two to the power of each of the integers in the subset, decrement the results, take the product, and multiply by two to the power of the sum of the subset.

# Brachylog, 45 bytes

∋P{;2↔^-₁}Q∈?∧P-₁~lċS⊆?∧S++₁Q-₁;2↔^R∧S↰₁ᵐ×;R×


Try it online!

This is a generator that can be used to get all results.

There are surely some shorter ways to express things.

### Explanation

∋P              P is an element of the input
{;2↔^-₁}Q     Q = 2^P - 1
Q∈?   Q is an element of the input
∧
P-₁~lċS         S is a list of length P-1
S⊆?       S is a subset of the input
∧
S++₁Q           The elements of S sum to Q-1
Q-₁;2↔^R    R = 2^(Q-1)
∧
S↰₁ᵐ            Map the first predicate {;2↔^-₁} to each element of S
×           Multiply the result
;R×        Multiply with R


# 05AB1E, 3028 22 bytes

-6 bytes porting @JonathanAllan's Jelly answer, thanks to @Steffan:

æʒg>Do<‚Ã<yOk}εo<PyOo*


Original 30 28 bytes answer:

æεg>©å*®o<åyO®oÍ©Q*iyo<®oªP,


Outputs the results on separated lines to STDOUT.

Explanation:

æ            # Get the powerset of the (implicit) input-list
ʒ           # Filter the inner lists y by:
#  (implicitly push list y)
g          #  Pop and get the length of list y
>         #  Increase it by 1
D        #  Duplicate this length+1
o       #  Take 2 to the power this value
<      #  Decrease it by 1
‚     #  Pair them together: [length+1,2**(length+1)-1]
Ã    #  Only keep those values from the (implicit) input-list
<          #  Decrease each remaining value by 1
yO        #  Push the sum of list y
k       #  Get its 0-based index in the list of 0, 1, or 2 items
#  (only 1 is truthy in 05AB1E)
}ε          # After the filter: map over the remaining inner lists y:
o         #  Take 2 the power each value in the list
<        #  Decrease each by 1
P       #  Take the product of these values
yO     #  Then push the sum of list y
o    #  Take 2 to the power this sum as well
*   #  And multiply it to the earlier product
# (after which the list is output implicitly as result)

æ            # Get the powerset of the (implicit) input-list
ε           # For each over the inner lists y:
#  (implicitly push list y)
g          #  Pop and get the length of list y
>         #  Increase it by 1
©        #  Store it in variable ® (without popping)
å       #  Check that it's in the (implicit) input-list
*      #  Multiply this 0 or 1 to the (implicit) input-list,
#  so we either have a list of 0s, or the unmodified input-list
®          #  Push length+1 from variable ® again
o         #  Pop and push 2 to the power this
<        #  Decrease it by 1
å       #  Check if it's in the earlier list (the input-list or list of 0s)
yO         #  Push the sum of the current list y
®o       #  Push 2 to the power ® again
Í      #  Decrease it by 2
©     #  Store this as new variable ® (without popping)
Q    #  Check if it's equal to the sum
*i         #  If both checks are truthy:
y        #   Push list y again
o       #   Calculate 2 to the power each value in the list
<      #   Then decrease each by 1
®     #   Push variable ®: 2**(length+1)-2
o    #   Get 2 to the power this value
ª   #   Append it to the list
P  #   Get the product of this list
, #   Pop and print it to STDOUT with trailing newline

• 22 bytes by porting Jonathan Allan's answer. Aug 25 at 19:32
• @Steffan Thanks! :) Aug 25 at 19:59