text:
Amount of permutations on an NxNxN Rubik's Cube
Introduction:
A 3x3x3 Rubik's Cube has \$43,252,003,274,489,856,000\$ possible permutations, which is approximately 43 quintillion. You may have heard about this number before, but how is it actually calculated?
A 3x3x3 Rubik's Cube has six sides, each with nine stickers. Looking at the (external) pieces instead of stickers however, we have six center pieces; eight corners pieces; and twelve edge pieces. Since the centers can't be moved, we can ignore those in the calculations. As for the corners and edges:

There are \$8!\$ (\$40,320\$) ways to arrange the eight corners. Each corner has three possible orientations, although only seven (of the eight) can be oriented independently; the orientation of the eighth/final corner depends on the preceding seven, given \$3^7\$ (\$2,187\$) possibilities.
There are \$\frac{12!}{2}\$ (\$239,500,800\$) ways to arrange the twelve edges. The halve from \$12!\$ is because edges must always be in an even permutation exactly when the corners are. Eleven edges can be flipped independently, with the flip of the twelfth/final edge depending on the preceding eleven, given \$2^{11}\$ (\$2,048\$) possibilities.

Putting this together, we have the following formula:
$$8!×3^7×\frac{12!}{2}×2^{11} = 43,252,003,274,489,856,000$$
Source: Wikipedia - Rubik's Cube Permutations
Although this may already look pretty complex, it's still rather straight-forward for a 3x3x3 Cube. For even cubes the formula is slightly different; this is the formula for a 4x4x4 Cube for example:
$$\frac{8!×3^7×24!^2}{24^7} = 7,401,196,841,564,901,869,874,093,974,498,574,336,000,000,000$$
Which is approximately 7.40 quattuordecillion on the short scale.
And for larger NxNxN Cubes (i.e. the current World Record 33x33x33) the formula will be extended quite a bit. To not make this introduction too long however, I put these links here instead, where the permutations of the 4x4x4 Cube and some other sized NxNxN Cubes are explained with a resulting formula:

2x2x2
4x4x4
5x5x5
6x6x6
7x7x7
33x33x33

You may be wondering by now: is there a general formula based on \$N\$ for any \$N\$x\$N\$x\$N\$ Cube? There certainly is. Here are three completely different algorithms, all giving the exact same results based on \$N\$:
1: Chris Hardwick's Formula:
$$\frac{(24×2^{10}×12!)^{N\pmod2}×(7!×3^6)×(24!)^{\lfloor{\frac{1}{4}×(N^2-2×N)}\rfloor}}{(4!)^{6×\lfloor{\frac{1}{4}×(N-2)^2}\rfloor}}$$
Try it on WolframAlpha.
2: Christopher Mowla's trig Formula:
$$8!×3^7×\left(\frac{24!}{(4!)^6}\right)^{\frac{1}{4}×((N-1)×(N-3)+\cos^2(\frac{N×\pi}{2}))}×(24!)^{\frac{1}{2}×(N-2-\sin^2(\frac{N×\pi}{2}))}×(12!×2^{10})^{\sin^2(\frac{N×\pi}{2})}×\frac{1}{24^{\cos^2(\frac{N×\pi}{2})}}$$
Try it on WolframAlpha.
3: Christopher Mowla's primes Formula:
$$2^{\frac{1}{2}×(2×N×(N+7)-17-11×(-1)^N)}×3^{N×(N+1)+2}×5^{\frac{1}{2}×(2×N×(N-2)+1+(-1)^N)}×7^{\frac{1}{8}×(6×N×(N-2)+3+5×(-1)^N)}×11^{\frac{1}{4}×(2×N×(N-2)-1+(-1)^N)}×96577^{\frac{1}{8}×(2×N×(N-2)-3+3×(-1)^N)}$$
where \$96577\$ is \$(13×17×19×23)\$.
Try it on WolframAlpha.
Source: Cubers-reddit - Mathematical Counting Formulas of Number of Positions, God's Number, etc.
Challenge:
Choose and implement one of these three formulas (or your own derivative), which given an input-integer \$N\$ in the range \$[2,100]\$, outputs the correct result.
Challenge rules:

You are free to use another formula besides these three, but keep in mind that these three are proven to be correct. If you use another formula, please add a link of where you got it from (or if you come up with it yourself add an in-depth explanation). And I will check for all integers in the range if the output is correct. Perhaps inspiration could be found in the oeis for this sequence: A075152.
If your language automatically outputs a scientific output (i.e. \$1.401...×10^{45}\$ instead of the number after the 4x4x4 formula) this is allowed. But please add additional code to your answer to convert this scientific rounding to an exact output so the results can be verified, since rounding errors due to floating point precision during the execution of the formula in your code is not allowed - the actual result should be exact.
Your program/function should be correct for at least the inputs in the range \$[2,100]\$ (although, since \$N=100\$ already results in a huge-ass number, any larger \$N\$ will probably work as well if you're able to output this one correctly).
You are not allowed to loop over all possible permutations with a counter, since that would never output anything in a reasonable amount of time. Only the implementation of a formula (either one of the three provided, a derivative of one of those, or a completely new formula), or another method that will give the correct results in a reasonable amount of time (without hard-coding of course) is allowed. I thought about adding a restricted-time to enforce this, but I'm personally against restricted-time in combination with code-golf, so I won't. Still, please make sure your program gives the answers, and if it's too slow for TIO for some reason, add some screenshots with the output from your local machine as verification.

General rules:

This is code-golf, so shortest answer in bytes wins.
Don't let code-golf languages discourage you from posting answers with non-codegolfing languages. Try to come up with an as short as possible answer for 'any' programming language.
Standard rules apply for your answer with default I/O rules, so you are allowed to use STDIN/STDOUT, functions/method with the proper parameters and return-type, full programs. Your call.
Default Loopholes are forbidden.
If possible, please add a link with a test for your code (i.e. TIO).
Also, adding an explanation for your answer is highly recommended.

Test cases:
Here the test cases for \$N\$ in the range \$[2,10]\$ (feel free to use the WolframAlpha links above for larger test cases):
n=2
3674160

n=3
43252003274489856000

n=4
7401196841564901869874093974498574336000000000

n=5
282870942277741856536180333107150328293127731985672134721536000000000000000

n=6
157152858401024063281013959519483771508510790313968742344694684829502629887168573442107637760000000000000000000000000

n=7
19500551183731307835329126754019748794904992692043434567152132912323232706135469180065278712755853360682328551719137311299993600000000000000000000000000000000000

n=8
35173780923109452777509592367006557398539936328978098352427605879843998663990903628634874024098344287402504043608416113016679717941937308041012307368528117622006727311360000000000000000000000000000000000000000000000000

n=9
14170392390542612915246393916889970752732946384514830589276833655387444667609821068034079045039617216635075219765012566330942990302517903971787699783519265329288048603083134861573075573092224082416866010882486829056000000000000000000000000000000000000000000000000000000000000000

n=10
82983598512782362708769381780036344745129162094677382883567691311764021348095163778336143207042993152056079271030423741110902768732457008486832096777758106509177169197894747758859723340177608764906985646389382047319811227549112086753524742719830990076805422479380054016000000000000000000000000000000000000000000000000000000000000000000000000000000000

NOTE: Since this is a code-golf challenge, it basically boils down to: implement one of these three formulas (or a derivative / your own method that still produces the correct results) as short as possible.

Why is the 4x4x4 cube considered the "odd" one instead of the 3x3x3? :P

@Quintec Shhh.. you didn't see anything! ;p Fixed, thanks. :)

Doing this in x86-64 will be a fun challenge. I'll have to roll my own bigint (likely just a 256-bit or 512-bit int), and make it golfy.

Off topic, but are there formula's for "super cubes", where the center pieces are marked (usually by colors on their edges) so that they can be uniquely identified? For example 3x3x3 centers orientation would matter (but there is parity). 4x4x4 centers can be moved and oriented (and there is also parity).

@moonheart08 Note that with N=100, the result takes about 128 Kb (16 KB): https://www.wolframalpha.com/input/?i=for+n+%3D+100,+N%5BLog%5B2,+((242%5E1012!)%5EMod%5Bn,+2%5D+7!3%5E6+24!%5EFloor%5B(n%5E2+-+2+n)%2F4%5D)%2F4!%5E(6+Floor%5B(n+-+2)%5E2%2F4%5D)%5D,3%5D

Note that Mowla's "trig" formula just uses $\sin^2$ and $\cos^2$ to obfuscate floors.

@rcgldr Yes, the reddit link also has a bunch of other formulas for the permutations of a Super NxNxN Cube; N-sized Minxes (dodecahedron-shaped variations of the Cubes, like Kilominx, Megaminx, Master Kilominx, Gigaminx, etc. for 2x2x2, 3x3x3, 4x4x4, 5x5x5, etc. respectively); and N-sized Super Minxes.

@attinat, I think it's a more useful perspective to say that both the trig and the floors are obfuscating $N \bmod 2$

+1 for "huge-ass number". Puts me in mind of trolls from the Discworld.

@KevinCruijssen No... I falsely analyzed the formula. My mistake, sorry.

@attinat and Peter Taylor, I don't appreciate the insult.  If you all would have read the rest of the formulas in the list, the purpose of using trig functions was to be able to take the first derivative of the formula.  Sounds stupid, yeah?  But with that derivative, I found a formula which has given very precise predictions on the world record times for the 4x4x4, 5x5x5, and 6x6x6 cubes for the past 7 years.  I guess it was my fault that I didn't provide links to the posts in which I originally posted the formula.  But Peter Taylor was spot on with it being N mod 2 instead of floor.

@ChristopherMowla Don't take their comments personal. I'm amazed that you were able to find these formulas and made such accurate predictions in the first place, and your formulas were one of the inspirations for this challenge. This is code-golf however, so readability, performance, warnings, best practices, historical significance, and sometimes just common sense are all thrown overboard if it can save a single byte on an answer. ;) attinat and PeterTaylor simply suggested such a golf based on your formulas, since N mod 2 is quite a bit shorter to use in programming languages than trig.

@ChristopherMowla "I guess it was my fault that I didn't provide links to the posts in which I originally posted the formula." You could link them in this comment section, curious to read the reasoning behind it, and how you've used to to predict times.

I appreciate your comment, and I understand what you mean, but attinat could have used a different verb.  He should know that I didn't list the formulas with coders to implement this.  But also, thanks for the compliment.  Without my trig formula, I would not have been able to find the prime formula (in its current most-condensed form).  So even if I didn't mention the prediction formula, that would have been worth it on its own.

As far as the world record prediction formula, see this post, https://www.speedsolving.com/forum/threads/world-record-formula-for-nxnxn-cubes.36621/#post-739116.
I provide links to the origin of my first 1st derivative formula and to the trig polynomial which appears in the actual world record formula (which accurately approximates the power of 10 of the first derivative formula up to n = 31).
I later used the same "power of 10" idea for the "God's number prediction formula.  These prediction formulas just "happened to work out", but you're right.  It's interesting that they do.

I didn't have enough space to also mention that the two posts above this post https://www.speedsolving.com/forum/threads/calculating-permutations-on-nxnxn-rubiks-cube.22683/page-2#post-649727 contain the original explanations of the trig formula and that trig polynomial (the one up to n = 31).

I just realized that the LaTeX in those posts is messed up in the number of positions thread (due to the forum going through upgrades).  However, below is a link to a saved PDF which has the original format so that you can read it more easily.
https://drive.google.com/file/d/1yQuQYOIe9AMRR07lDDOuPeIo-9Z0qGwW/view?usp=sharing
(The PDF is perfectly clear if you download it and view from your device.)

I forgot to mention that I made the "Wolfram|Alpha friendly version" of my God's number for the nxnxn (lowerbound,upperbound) formula by using a ratio of the same trig polynomial expression.  (I just added this link to my Reddit post, BTW.)
(wolframalpha.com/input/?i=Table%5BTable%5BFloor%5B((Pi-3)(62n%5E2-56n-4Cos%5B(nPi)%2F2%5D%2B21Cos%5BnPi%5D-33))%2F(ProductLog%5B1%2Fx((Pi-3)(62n%5E2-56n-4Cos%5B(nPi)%2F2%5D%2B21Cos%5BnPi%5D-33))%5D)-1%5D,%7Bx,2,3%7D%5D,%7Bn,2,20%7D%5D)
(See the following topic to compare the lowerbounds.)
http://cubezzz.dyndns.org/drupal/?q=node/view/236

Wolfram Language (Mathematica), 59 bytes
f@n_:=(s=24^6)(24!/s)^(m=n-2)f@m
f@2=7!3^6
f@3=4!12!2^10f@2

Try it online!
uses Herbert Kociemba's algorithm found in OEIS page
here is the recursive formula:
a(1)=1; a(2)=7!*3^6; a(3)=8!*3^7*12!*2^10; a(n)=a(n-2)*24^6*(24!/24^6)^(n-2)
6 bytes saved by @Peter Taylor
one more byte saved by @Expired Data

65 bytes

@ExpiredData very nice! I was trying to do the same but it failed cause the order was different

The question doesn't require you to support f@1, so you can save 6 bytes. Obviously you'd also want to adjust your test framework to use Range[2,10].

@PeterTaylor great observation. thanks!

@CSM unfortynately f[3] needs to be defined.Otherwise the formula returns wrong results

x86 machine code, 119 bytes
Hexdump:
60 c6 02 02 33 db be 25 01 10 00 f6 c1 01 74 05
be 26 2a b2 36 33 ed 51 b1 06 33 ff 53 8a 04 1a
f6 e1 03 c7 b5 0a f6 f5 88 64 1a 02 66 98 8b f8
4b 79 ea 5b 43 43 f6 f5 66 89 02 84 c0 75 0c 60
8b fa 8d 72 01 8b cb f3 a4 61 4b 41 d1 ee 72 ca
75 f9 be 1d d4 0d 10 4d 79 be 59 49 49 8b e9 be
06 02 02 22 83 f9 02 73 ae c6 44 1a 01 00 80 0c
1a 30 4b 79 f9 61 c3

The function receives the number n in ecx, and a pointer to a string to fill in edx (i.e. fastcall convention).
Before I show the source code, some explanations on how it does the thing. It uses the recursive formula, which I wrote in the following way:
init = 2
m1 = 24^6 = 6*8*9*16*24*32*36
m2 = 24!/24^6 = 6*7*9*10*11*17*19*21*22*23*25*26*35
num(2) = init * 6*7*9*12*15*27
num(3) = init * 6*8*9*12*16*18*20*24*27*28*30*32*33*35*36
num(n+2) = num(n) * m1 * m2^n

So all the code should do is multiplication by small numbers. The numbers are in the range 6...36, which is small enough to be represented in a 32-bit bitmap. I actually don't store the bit which represents multiplication by 6 - this lets me arrange the code in a do-while loop, starting with unconditional multiplication by 6.
The big numbers are represented using decimal form - each byte is a value in the range 0...9, starting from the MSB.
The multiplication is performed from LSB to MSB; it assumes that the number of digits will grow by 2 for each multiplication. After doing multiplication by a small factor like 6, the number of digits may grow by only 1. So if MSB = 0, it moves the whole intermediate result left. It can actually happen that the number of digits doesn't grow at all, and then MSB will still be 0, but this problem will fix itself as the code proceeds to greater factors.
Because the multiplication code is large, I don't want to put it twice. I don't want to move it to a function either, because the machine code for calling a function is large. So I rearranged the outer loops in such a way that the multiplying code is needed only once.
C code:
void num(int n, char* x)
{
    *x = 2;
    int len = 1;
    int exp_i;
    uint32_t m32_1;
    int m1;
    int carry;
    int temp;
    int str_i;
    bool cf;

    if (n % 2 == 0)
    {
        m32_1 = 0x100125; // 6*7*9*12*15*27
    }
    else
    {
        m32_1 = 0x36b22a26; // 6*8*9*12*16*18*20*24*27*28*30*32*33*35*36
    }

    exp_i = 0;
    while (true)
    {
        for (; exp_i >= 0; --exp_i)
        {
            m1 = 6;
            cf = true;
        do_mult:
            carry = 0;
            for (str_i = len - 1; str_i >= 0; --str_i)
            {
                temp = x[str_i] * m1 + carry;
                x[str_i + 2] = temp % 10;
                carry = temp / 10;
            }
            len += 2;
            x[1] = carry % 10;
            carry /= 10;
            x[0] = carry;
            if (carry == 0)
            {
                --len;
                for (str_i = 0; str_i < len; ++str_i)
                    x[str_i] = x[str_i + 1];
            }

        shift_m1:
            ++m1;
            cf = m32_1 & 1;
            m32_1 >>= 1;
            if (cf)
                goto do_mult;
            if (m32_1)
                goto shift_m1;

            m32_1 = 0x100dd41d; // 24!/24^6 = 6*7*9*10*11*17*19*21*22*23*25*26*35
        }
        --n;
        --n;
        exp_i = n;
        if (n < 2)
            break;
        m32_1 = 0x22020206; // 24^6

    }
    x[len] = 0;
    for (str_i = len - 1; str_i >= 0; --str_i)
    {
        x[str_i] += '0';
    }
}

Disassembly:
60                     pushad;
C6 02 02               mov byte ptr [edx], 2; // edx = x
33 DB                  xor ebx, ebx; // ebx = len - 1
BE 25 01 10 00         mov esi, 0x100125; // esi = m32_1
F6 C1 01               test cl, 1;
74 05                  jz skip1;
BE 26 2A B2 36         mov esi, 0x36b22a26; // esi = m32_1
                   skip1:
33 ED                  xor ebp, ebp; // ebp = exp_i

                   loop_n:

51                     push ecx;
                   loop_exp_i:
B1 06                  mov cl, 6; // cl = m1
                   do_mult:
33 FF                  xor edi, edi; // edi = carry
53                     push ebx; // ebx = str_i
                   loop_str_i:
8A 04 1A               mov al, [edx + ebx];
F6 E1                  mul cl;
03 C7                  add eax, edi;
B5 0A                  mov ch, 10;
F6 F5                  div ch;
88 64 1A 02            mov [edx + ebx + 2], ah;
66 98                  cbw;
8B F8                  mov edi, eax;
4B                     dec ebx;
79 EA                  jns loop_str_i;

5B                     pop ebx; // ebx = len - 1
43                     inc ebx;
43                     inc ebx;
F6 F5                  div ch;
66 89 02               mov [edx], ax;

84 C0                  test al, al;
75 0C                  jnz skip2;

60                     pushad;
8B FA                  mov edi, edx;
8D 72 01               lea esi, [edx + 1];
8B CB                  mov ecx, ebx;
F3 A4                  rep movsb;
61                     popad;
4B                     dec ebx;
                   skip2:

                   shift_m1:
41                     inc ecx;
D1 EE                  shr esi, 1;
72 CA                  jc do_mult;
75 F9                  jnz shift_m1;

BE 1D D4 0D 10         mov esi, 0x100dd41d;

4D                     dec ebp;
79 BE                  jns loop_exp_i;

59                     pop ecx; // ecx = n
49                     dec ecx;
49                     dec ecx;
8B E9                  mov ebp, ecx;
BE 06 02 02 22         mov esi, 0x22020206;
83 F9 02               cmp ecx, 2;
73 AE                  jae loop_n;

C6 44 1A 01 00         mov byte ptr [edx + ebx + 1], 0;
                   loop_to_ascii:
80 0C 1A 30            or byte ptr [edx + ebx], '0';
4B                     dec ebx;
                       dec         ebx
79 F9                  jns loop_to_ascii;

61                     popad;
C3                     ret;

The running time for n = 100 is about 4 seconds, and the result is a number with 38416 digits:
23491019577617 (many digits here) ... (many zeros here) 0000000000000000

05AB1E, 38 bytes
Initial attempt.
Uses Chris Hardwick's Formula.
Will attempt to golf further and explain when I have time.
24©To12!PIÉm7!729®!InI·-4÷mP®IÍn4÷6*m÷

Try it online!

Julia 1.0, 83 76 bytes
n->^(24576*~12,n%2)*3^6*~7(~24)^((m=n-2)n÷4)/24^(m^2÷4*6)
~n=prod(big,1:n)

Try it online!
Uses Chris Hardwick's Formula. Takes input as Big integer.
Thanks to H.PWiz for -7 bytes

~=n->factorial(big(n)) -> ~n=prod(big,1:n) and (24576*~12)^(n%2) -> ^(24576*~12,n%2)

Why do you use ~=n-> instead of ~n=?

@H.PWiz, simply because I didn't even know it would work this way, and also didn't notice that in your previous comment :)

Python 2, 72 bytes
lambda n:3674160*61600**(n%2)*24**(~-n/2*6)*0xb88d4641131f0**(n*(n-2)/4)

Try it online!
Saved 4 bytes by copying n*(n-2)/4 from Neil.

Wolfram Language (Mathematica), 67 bytes
Using Chris Hardwick's Formula.
(12!24576)^Mod[#,2]7!729(24!)^⌊#(#-2)/4⌋/24^(6⌊(#-2)^2/4⌋)&

Try it online!

63 bytes. Was about to post mine based on Mowla's "trig" formula, but it turned out to be identical to this...

JavaScript (Node.js), 81 bytes
Herbert Kociemba's recursive formula. Takes a BigInt as input.
f=n=>[1n,3674160n,322252536375n<<27n][--n]||f(--n)*0xb640000n*0xb88d4641131f0n**n

Try it online!

JavaScript (Node.js),  102 98  96 bytes
Chris Hardwick's formula. Takes a BigInt as input.
n=>(n&1n?1403325n<<25n:4n)*918540n*0x83629343d3dcd1c00000n**(n*n-n-n>>2n)/24n**(6n*(n*n/4n-~-n))

Try it online!

JavaScript (Node.js), 77 75 73 bytes
n=>0xb88d4641131f0n**(n*(n-2n)/4n)*13824n**n*851558400n**(n%2n)*315n>>14n

Try it online! Based on Christopher Mowla's formula. Takes a BigInt as input. Test harness shamelessly stolen from @Arnauld. 0xb88d4641131f0n is 3246670537110000n in decimal. Explanation: I started with the last prime exponent and simplified it to n*(n-2n)/4n (this is integer division, so I don't need the adjustment for odd numbers). I then examined the other primes to see whether their exponents were related to this value (which I shall refer to as o), and found that they were after a fashion, if I allowed the use of the parity of n (which I shall refer to as p). The formulae for the exponents are as follows:
23:       o
19:       o
17:       o
13:       o
11:      2o +   p
 7:      3o +   p +  1
 5:      4o +  2p +  1
 3: 3n + 4o +  3p +  2
 2: 9n + 4o + 14p - 14

The powers can then be grouped by exponent so for example p is the exponent of 11*7*5**2*3**3*2**14.

Racket, 151 141 bytes
-7 bytes thanks to fede s.!
(λ(n[e expt])(/(*(e 11771943321600(modulo n 2))3674160(e 620448401733239439360000(floor(/(*(- n 2)n)4))))(e 24(*(floor(/(sqr(- n 2))4))6))))

Try it online!
The longest answer using Chris Hardwick's Formula :)

You could change the define for a λ (2 bytes), and use a default value for an extra parameter to save 3 more bytes from the three expt calls: (λ(n[e expt])...(e ...)...).

Python 2, 122 bytes
import math
f=math.factorial
x=lambda n:(1,f(7)*729,f(8)*3**7*f(12)*1024)[n-1]if n<4else x(n-2)*24**6*(f(24)/24**6)**(n-2)

Try it online!
Uses the Herbert Kociemba recursive method.
-2 bytes thanks to Herman L

2 bytes can be saved by replacing 3**6 with 729 and 2**10 with 1024 TIO

Hey, nice catch!

Defining your own factorial function is 3 bytes shorter

103 bytes with hardcoded factorial values

Jelly,  39  38 bytes
I feel like I've missed some golfs, but...
12!×⁽^K*Ḃɓ_2×ṭ¥⁸:4×1,6“ð¥‘!¤*:/×“9Ḟɠ’×

A monadic Link implementing Chris Hardwick's Formula.
Try it online! Or see the test-suite (n=[1..33]).

CJam (47 bytes)
qi[1_7m!Z6#*_3*Cm!*2D#*]{2-_j24_m!\6#:P/@#*P*}j

Online demo
This implements Herbert Kociemba's recursion from OEIS: $$a(n) = \begin{cases}
1 & \textrm{ if } n \in \{0,1\} \\
7! \times 3^6 & \textrm{ if } n=2 \\
a(n-1) \times 3\times 12!\times 2^{13} & \textrm{ if } n=3 \\
a(n-2) \times \left(\frac{24!}{24^6}\right)^{n-2} \times 24^6 & \textrm{ if } n>3 \end{cases}$$
using CJam's memoised recursion operator j. I've ordered the terms in the MathJax block in the same order as in the code to make the correspondence easy to verify for those who read CJam: any further dissection is not going to shed more light.

Perl 6, 85 bytes
{0xAB4DE800000**($_%2)*3674160*([*] 1..24)**($_*($_-2)div 4)/24**(($_-2)**2 div 4*6)}

Try it online!

Haskell, 86 85 74 bytes
-1 byte saved thanks to H.PWiz

-11 bytes saved thanks to Max Yekhlakov
a=24^6
r 2=3674160
r 3=r 2*a*61600
r n=r(n-2)*a*div(product[2..24])a^(n-2)

Try it online!

24576 is shorter than 2^13*3

74 bytes Try it online!

Husk, 51 48 44 bytes
-4 bytes thanks to H.PWiz
÷^*6÷4□-2⁰Π4*^÷4-D⁰□⁰Π24*729*Π7^%2⁰*24576Π12

Try it online!
This is Chris Hardwick's Formula. Also, this is my first husk program, so any tips would be well appreciated.

Here's an easy 2 bytes: ÷^*6÷4□-2⁰Π4*^÷4-D⁰□⁰Π24*729*Π7^%2⁰*24*1024Π12

Or, better yet, ÷^*6÷4□-2⁰Π4*^÷4-D⁰□⁰Π24*729*Π7^%2⁰*24576Π12

Icon, 125 110 bytes
procedure f(n)
q:=1;every q*:=1 to 24
return 11771943321600^(n%2)*5040*3^6*q^(n*(t:=n-2)/4)/24^(6*(t^2/4))
end

Try it online!

Jelly, 43 bytes
_2²:4×6*@24
²_Ḥ:4;ḂU
“€ð‘!×⁽^K,1*ÇP×“9Ḟɠ’:Ñ

Try it online!

C (gcc) -lgmp, 279 bytes
#include "gmp.h"
#define s mpz_init_set_str
#define m(X)mpz_##X
f(int N,m(t)_){m(t)x;m(init)(x);m(init_set_str)(_,N&1?"3LFbOUwC":"1",62);m(mul_si)(_,_,3674160);m(fac_ui)(x,24);m(pow_ui)(x,x,(N*N-2*N)/4);m(mul)(_,_,x);m(set_si)(x,24);N-=2;m(pow_ui)(x,x,6*N*N/4);m(tdiv_q)(_,_,x);}

Try it online!

Suggest N--*--N/4 instead of (N*N-2*N)/4 and remove N-=2 and #define s mpz_init_set_str

TI-BASIC, 63 62 bytes, (noncompeting)
{fPart(.5Ans),1,1,-6}int(4⁻¹{8,4,Ans²-2Ans,(Ans-2)²:prod({9*11!2^15,7!3^6,24!,24}^Ans

Expression which takes input as an integer on Ans. Implementation of Chris Hardwick's formula. Noncompeting because the hardware it runs on will only store up to 16 decimal places, so the answer will never be 100% accurate.
Explanation:
{fPart(.5Ans),1,1,-6}              # the list {(N (mod 2))/2,1,1,-6}
                                   # implicitly multiplied by
int(4⁻¹{8,4,Ans²-2Ans,(Ans-2)²     # the list {2,1,⌊¼(N²-2N)⌋,⌊¼(N-2)²⌋}
:                                  # store this list of the formula's exponents as Ans
     {9*11!2^15,7!3^6,24!,24}      # list of the formula's bases
                             ^Ans  # raised to their exponents
prod(                              # multiplied together
                                   # implicit print

Python 2, 92 bytes
lambda n:0xab4de800000**(n%2)*3674160*0x83629343d3dcd1c00000**(n*(n-2)/4)/24**((n-2)**2/4*6)

Try it online!

C++, 187 185 180 176 195(there was a bug) 193 175 bytes (with help from ceiling cat)
This uses the GMP C++ wrapper (GNU multi-precision library), and the formula used by @J42161217 (https://codegolf.stackexchange.com/a/183381/55953).
Use g++ -g rubix.cpp -lgmp -lgmpxx to compile and link
#include <gmpxx.h>
#define R return
using z=mpz_class;z p(z a,z b){z c=1;while(b--)c*=a;R c;}z r(z n){if(n==2)R 3674160;if(n==3)R z("pX4dIaR7jDk",62);R r(n-2)*p(24,6)*p(z("ErvSErbeq",62),n-2);}

ungolfed, with test code
#include <gmpxx.h>
#include <iostream>
mpz_class p(mpz_class a, mpz_class b) // returns a to power of b. Only works for b  = positive integer
{
    mpz_class c=1;

    while(b--)
        c*=a;

    return c;
}

mpz_class r(mpz_class n) // returns the rubix permutations for a cube of size n
{
    if(n==2)
        return 3674160; // 7!*3^6;

    if(n==3)
        return z("pX4dIaR7jDk",62); // 43252003274489856000 = 8!*3^7*12!*2^10

    return r(n-2) * p(24,6) * p(z("ErvSErbeq", 62), n-2);

    // "ErvSErbeq"base 62 = 3246670537110000 = (24!/24^6)
}

main()
{
    for(int i=2; i<34; i++)
        std::cout<<i<<'\t'<<r(i) << std::endl;
}

https://tio.run/##PZAxb4MwEIV3foWVDrETqBpARMImWZqha7t0iFQZ4xC3xrg2tJERf73UIVXfcE937zvpdEzrqGZsmu6EYrKvOKkbfbncn3dBb4WqgSsa7d6YpNZiBzR0gIYOlGhwgBUb/H0WksMyihBbFRQb3vVGAYZHB4xnFRr@Rqoo4n2SbdNN9pD7Jtk7uNCvafVEn7fvjx@LMItRbqCKYrTSME7D7OoeOpivl4Mp@eeMhFcAj//3AiJa2xlOm13QUKEgCoYAeJ1aA4XqgChiDARJUl/XazRnXrar8py1fUeIIGR57JaE@AUECLllXFUSB2Mw/bCTpLWdIjm/5ua/

Could you perhaps add a screenshot of the n=10 test case, so I can verify that it works? I guess there isn't any way to make this work on the C++ (clang) or C++ (gcc) TIO due to the used library?

argg. It's not working for odd values of n

Thanks for the screenshot, and good that you've been able to pinpoint the error and fix it. +1 from me. :)

182 bytes

Ta @ceilingcat. The #define return isn't needed any more, as there's only two return points


text:
Maximum summed subsequences with non-adjacent items
Introduction:
Inspired by these two SO questions (no doubt from the same class): print the elements in the subarray of maximum sum without adjacent elements java and Maximum sum of non adjacent elements of an array, to be printed.
Challenge:
Given a list of integers, output a subsequence consisting of non-adjacent elements that have the highest sum. Here some examples:

[1,2,3,-1,-3,2,5] would result in [1,3,5] (with a sum of 9) at the 0-based indices [0,2,6].
[4,5,4,3] would result in either [4,4] (with a sum of 8) at the 0-based indices [0,2] or [5,3] (also with a sum of 8) at the 0-based indices [1,3].
[5,5,10,100,10,5] would result in [5,100,5] (with a sum of 110) at either the 0-based indices [0,3,5] or [1,3,5].

What's most important about these examples above, the indices containing the elements are at least 2 apart from each other. If we look at the example [5,5,10,100,10,5] more in depth: we have the following potential subsequence containing non-adjacent items; with their indices below it; with their sums below that:
[[5],[10],[100],[10],[5],[5],[100,5],[10,5],[10,10],[5,5],[5,10],[5,100],[5,5],[5,10],[5,100],[5,10],[5,100,5],[5,100,5],[5,10,5],[5,10,10]]   // non-adjacent subsequences
[[5],[ 4],[  3],[ 2],[1],[0],[  3,5],[ 2,5],[ 2, 4],[1,5],[1, 4],[1,  3],[0,5],[0, 4],[0,  3],[0, 2],[1,  3,5],[0,  3,5],[0, 2,5],[0, 2, 4]]   // at these 0-based indices
[  5,  10,  100,  10,  5,  5,    105,    15,     20,   10,    15,    105,   10,    15,    105,    15,      110,      110,      20,       25]   // with these sums
                                                                                ^         ^                        // and these two maximums

Since the maximum sums are 110, we output [5,100,5] as result.
Challenge rules:

You are allowed to output key-value pairs of the index + value. So instead of [5,100,5] you can output [[0,5],[3,100],[5,5]] or [[1,5],[3,100],[5,5]] as result (or [[1,5],[4,100],[6,5]]/[[2,5],[4,100],[6,5]] when 1-based indexing is used instead of 0-based).

If you use key-value pairs, they can also be in reverse or random order, since it's clear which values are meant due to the paired index.
Outputting just the indices without values isn't allowed. It should either output the values, or the values/indices as key-value pairs (or two separated lists for 'keys' and 'values' of the same size if key-value pairs are not possible in your language of choice).

You are allowed to output all possible subsequences with the maximum sum instead of just one.
As you can see from the examples, the input-list can contain negative and duplicated values as well. You can assume the input-integers are within the range \$[-999,999]\$.
The output-list cannot be empty and must always contain at least one element (if a list would only contain negative values, a list containing the single lowest negative value would then be output as result - see last two test cases).
If there is one possible output but for multiple different indices, it's allowed to output both of them even though they might look duplicates. (i.e. the example above, may output [[5,100,5],[5,100,5]] for both possible index-combinations).

Test cases:
Input:                   Possible outputs:       At 0-based indices:     With sum:

[1,2,3,-1,-3,2,5]        [1,3,5]                 [0,2,6]                 9
[4,5,4,3]                [4,4]/[5,3]             [0,2]/[1,3]             8
[5,5,10,100,10,5]        [5,100,5]               [0,3,5]/[1,3,5]         110
[10]                     [10]                    [0]                     10
[1,1,1]                  [1,1]                   [0,2]                   2
[-3,7,4,-2,4]            [7,4]                   [1,4]                   11
[1,7,4,-2]               [7]                     [1]                     7
[1,2,-3,-4,5,6,-7]       [2,6]                   [1,5]                   8
[800,-31,0,0,421,726]    [800,726]/[800,0,726]   [0,5]/[0,3,5]/[0,2,5]   1526
[-1,7,8,-5,40,40]        [8,40]                  [2,4]/[2,5]             48
[-5,-18,-3,-1,-10]       [-1]                    [3]                     -1
[0,-3,-41,0,-99,-2,0]    [0]/[0,0]/[0,0,0]       [0]/[3]/[6]/[0,3]/
                                                  [0,6],[3,6]/[0,3,6]    0

If there is more than one identical set (but from different indices) is it ok to list all of them? e.g. [5,100,5] twice for your third example.

powersetis a set of subsets isn't it? but it looks like you are returning a set of subsequences? [4,5,4,3] would result in either [4,4] where [4,4] is clearly not a set.

@NickKennedy Yes, that's allowed, will edit the description.

@ExpiredData I initially had sublists as term in the challenge, but then I remembered that 05AB1E/Jelly/Pyth and alike call non-adjacent sublists powersets, so I changed it to that term. For example, here the powerset of the list [1,2,3,4,5] in 05AB1E, for which the resulting inner lists also aren't in order.

@Arnauld Yes, if the values are key-value pairs with their index it's clear which indexed values are meant in the input, so they can be in any order. Will also edit this into the challenge description.

Just to be sure: outputting the indices isn't an option, is it?

@Shaggy No sorry, either the values, or the key-value pairs with index/value. Sorry. I'll edit the challenge to mention it.

@KevinCruijssen: That whole output in your link is the power set of [1, 2, 3, 4, 5]. The individual sublists aren't all power sets. "Power set" is not a correct term to describe what this challenge is asking for.

@user2357112 What would be the correct term here then? Powerset-sublists or something like that? With just the term sublists I think of consecutive elements.

The classical term is "subsequence". This has the same problem of people thinking of contiguous subsequences, though. I would say "subset" if we were actually working with sets here, but these are definitely sequences - order matters and duplicates are allowed.

@user2357112 Changed it to subsequences.

FYI, the first linked SO question is no more.

@Arnauld Oh. Ah well, luckily inspirational links aren't necessary for the actual challenge, nor for it to be self-contained. Will leave it like this, even though the SO post doesn't exist anymore.

Jelly, 8 bytes
ŒPm€2SÞṪ

Try it online!
This is a function submission (it's almost usable as a full program, which is what the TIO link does, but the output format is a little weird; when used as a function, though, it just outputs the list of elements from the subsequence).
Explanation
I found a much terser algorithm (loosely inspired by my previous Brachylog answer), saving 6 bytes over the Jelly and Brachylog answers and 3 bytes over the current leading answer.
ŒPm€2SÞṪ
ŒP         On the set of all subsequences of {the input},
  m 2        remove every second element of
   €         each subsequence,
      ÞṪ     then {return} the resulting subsequence with the maximal
     S       sum

The basic idea is that we don't look for the subsequence we want directly; rather, we find a subsequence with an additional element inserted in between each of the elements we want. This is always possible if the elements of the subsequence we want are not adjacent in the original list (because there'll always be at least one in between to insert), and never possible if the elements of the subsequence we want are adjacent in the original list. So all we have to do is look at all the subsequences of the input, then take alternating elements from them, to get a list of all the non-adjacent-item subsequences of the input.
This algorithm will return many of the non-adjacent-item subsequences more than once, but that doesn't matter, because we're just looking for a maximal element, and repeating elements of the set we're taking a maximum of has no influence on the maximum.

Husk, 11 bytes
►Σ†!¹mü≈tṖŀ

Try it online!
Explanation
►Σ†!¹mü≈tṖŀ  Implicit input, say L=[4,5,3,4].
          ŀ  Indices: [1,2,3,4]
         Ṗ   Powerset: [[],[1],[2],[1,2],..,[1,2,3,4]]
        t    Tail (remove the empty list): [[1],[2],[1,2],..,[1,2,3,4]]
     m       For each,
      ü      de-duplicate by
       ≈     differing by at most 1.
             For example, [1,2,4] becomes [1,4].
  †          Deep map
   !¹        indexing into L: [[4],[5],[4],..,[5,4],[4,3]]
►            Maximum by
 Σ           sum: [5,4]

Just a comment to inform you someone outgolfed you, so I've switched the accepted answer. Sorry for the minor loss of rep as a result. :)

R, 136 125 bytes
function(l,G=unlist(Map(combn,list(y<-seq(a=l)),y,c(function(x)'if'(all(diff(x)>1),l[x],-Inf)),F),F))G[which.max(Map(sum,G))]

Try it online!
-6 bytes thanks to digEmAll, who incidentally also outgolfed me.
Returns the shortest subsequence as a list, breaking ties on lexicographically first by indices.
Brute-force generates all index subsequences, then Filters for those that are non-adjacent, i.e., where all(diff(x)>1). Then subsets [ into l using these indices, selecting [[ the first one where the sum is the max (which.max).
I'm pretty sure this is the first R answer I've ever written that uses Filter! sad, Filter is ungolfy, no wonder I've never used it...

130

Haskell, 60 bytes
snd.([]%)
r%(h:t)=max(r%t)$(r++[h])%drop 1t
r%_=(sum r<$r,r)

Try it online!
The helper function % recursively branches on choosing whether the include the first element and drop the second, or to skip the first element. It takes the maximum of all outcomes, which are tuples whose first element is the sum, and whose second element is the corresponding list which is extracted for the output.
To handle the rule that the empty list is disallowed even if it would have the smallest trick, we do a cute trick of writing sum r<$r rather than sum r.This makes a list whose elements all are sum r and whose length is that of r. That way, when we choose the maximum, we prioritize any list over an empty r, but otherwise comparisons depend on the first element which is sum r .

05AB1E, 14 bytes
Saved 1 byte thanks to Kevin Cruijssen
ā<æʒĆ¥≠W}èΣO}θ

Try it online!
or as a Test Suite
Explanation
ā<               # push [0 ... len(input)-1]
  æ              # compute powerset
   ʒ    }        # filter, keep lists where:
      ≠W         # no element is 1 in the
     ¥           # deltas
    Ć            # of the list with the head appended
         è       # index into the input with each
          ΣO}    # sort by sum
             θ   # take the last element

You may not be happy, but it's still 4 bytes shorter than my initial solution. ;) And you can golf 1 more changing ¤ª to Ć.

@KevinCruijssen: Oh yeah! For some reason I had convinced myself I needed a repeat element at the end. Thanks!

Brachylog (v2), 14 bytes
{~ba~c∋₁ᵐ}ᶠ+ᵒt

Try it online!
Function submission; input from the left, output from the right, as usual. Very slow; a five-element list is probably the maximum for testing on TIO.
{~ba~c∋₁ᵐ}ᶠ+ᵒt
 ~b              Prepend an arbitrary element to the input
   a             Take a prefix or suffix of the resulting list
    ~c           Ordered partition into contiguous sublists
      ∋₁         Take the second element
        ᵐ          of each sublist
{        }ᶠ      Find all possible ways to do this
           +ᵒ    Sort by sum
             t   Take the greatest

The results we get from prefixes aren't incorrect, but also aren't interesting; all possible results are generated via taking a suffix (which is possibly the list itself, but cannot be empty), but "suffix" is more verbose in Brachylog than "prefix or suffix", so I went with the version that's terser (and less efficient but still correct). The basic idea is that for each element we want in the output list, the partition into contiguous sublists needs to place that element and the element before into the same sublist (because the element is the second element of the sublist), so two consecutive elements can't appear in the result. On the other hand, it's fairly clear that any list without two consecutive elements can appear in the result. So once we have all possible candidate lists, we can just take the sums of all of them and see which one is largest.

J, 47 bytes
>@{:@(<@#~/:+/@#~)1(-#~0=1 1+/@E."1-)[:i.&.#.=~

Try it online!
-7 bytes thanks to FrownyFrog
original
J, 54 bytes
[:(>@{:@/:+/&>)]<@#~"1[:(#~0=1 1+/@E."1])[:#:@}.@i.2^#

Try it online!

47

Thanks you very much @FrownyFrog

Jelly, 16 14 bytes
JŒPḊf’$ÐḟịµSÞṪ

Try it online!
Thanks to @EriktheOutgolfer for saving 2 bytes!

14 bytes.

JavaScript (ES6),  138 132 130 129  126 bytes
Outputs key-value pairs.
a=>a.reduce((a,x,i)=>[...a,...a.map(y=>[[x,i],...y])],[[]]).map(m=a=>a.some(s=p=([v,i])=>p-(s=~~s+v,p=i)<2)|s<m||(r=a,m=s))&&r

Try it online!
Step 1
We first compute the powerset of the input with \$[value, index]\$ pairs.
a.reduce((a, x, i) => // for each value x at position i:
  [                   //   update a[] to a new array consisting of:
    ...a,             //     all previous entries
    ...a.map(y =>     //     for each value y in a[]:
      [[x, i], ...y]  //       append [x, i], followed by all original entries
    )                 //     end of map()
  ],                  //   end of new array
  [[]]                //   start with a = [[]]
)                     // end of reduce()

Step 2
We then look for the maximum sum \$m\$ among these sets, discarding sets with at least two adjacent elements. The best set is stored in \$r\$.
.map(m =              // initialize m to a non-numeric value
  a =>                // for each entry a[] in the powerset:
  a.some(s = p =      //   initialize s and p to non numeric values
    ([v, i]) =>       //   for each value v and each index i in a[]:
    p - (             //     compute p - i
      s = ~~s + v,    //     add v to s
      p = i           //     update p to i
    ) < 2             //     if p - i is less than 2, yield true
  ) |                 //   end of some()
  s < m ||            //   unless some() was truthy or s is less than m,
  (r = a, m = s)      //   save a[] in r[] and update m to s
) && r                // end of map(); return r[]

T-SQL, 122 119 118 bytes
Input is a table variable.
This query picks all elements from the table variable, combining these with all non-adjacent elements with higher position values and show the text generated for the highest sum of these values.
WITH C(y,j,v)as(SELECT*,x*1FROM @
UNION ALL
SELECT y+','+x,i,v+x
FROM @ JOIN C ON~-i>j)SELECT
TOP 1y FROM C ORDER BY-v

Try it online ungolfed

Haskell, 81 80 bytes
snd.maximum.map((,)=<<sum).tail.f
f(a:b:c)=f(b:c)++map(a:)(f c)
f a=[]:map(:[])a

Try it online!
f builds all valid subsequences by either skipping the next element (f(b:c)) or using it and skipping the next (map(a:)(f c)) and recursively work on the rest. For the result, build all subsequences (f), drop the empty subsequence (which occurs first in the list: tail), make pairs (<sum>,<subsequence>) (map((,)=<<sum)), find the maximum (pairs are compared in lexicographical order) -> maximum) and drop the sum (snd).
Edit: -1 byte thanks to @Lynn.

map(:[])a is a byte shorter than (pure<$>a) ^^

R, 108 107 bytes
function(v,M=-Inf){for(j in J<-seq(a=v))for(i in combn(J,j,,F))if(all(diff(i)>1)&sum(v[i])>sum(M))M=v[i]
M}

Try it online!
-1 thanks to @Giuseppe

Pyth, 19 bytes
esDm@LQdtf!q#1.+TyU

Try it online here, or verify all the test cases at once here.
esDm@LQdtf!q#1.+TyUQ   Implicit: Q=eval(input())
                       Trailing Q inferred
                  UQ   Generate range [0-len(Q))
                 y     Take the powerset of the above
         f             Filter keep elements of the above, as T, using:
              .+T        Take differences of consecutive elements of T
           q#1           Keep those differences equal to 1
          !              Logical NOT - empty lists evaluate to true, populated ones to false
                       Result of the filter is those sets without consecutive numbers
        t              Drop the first element (empty set)
   m                   Map the remaining sets, as d, using:
     L d                 For each element of d...
    @ Q                  ... get the element in Q with that index
 sD                    Order the sets by their sum
e                      Take the last element, implicit print

Wolfram Language (Mathematica), 144 bytes
If[Max[a=#]>(d=Max[m=Tr/@(g=a[[#]]&/@Select[Subsets[Range[x=Length@#],{2,x}]&@#,FreeQ[Differences@#,1]&]&@a)]),{Max@a},g[[#]]&@@@Position[m,d]]&

Try it online!

Wolfram Language (Mathematica), 70 63 bytes
MaximalBy[Select[q=Rest@Subsets@#,!FreeQ[q,#~Riffle~_]&],Tr,1]&

Try it online!
High-level search
          Select[q=Rest@Subsets@#,                     ]        (*choose nonempty subsets of the input such that*)
                                  !FreeQ[q,          ]&         (*there exists a subset of the input which matches*)
                                           #~Riffle~_           (*this list, with an item inserted between adjacent elements*)
MaximalBy[                                              ,Tr,1]& (*and return one with the greatest total*)

,1 is required so as not to inadvertently return invalid sets (otherwise, for example, an input of {1,1,1} would result in an output of {{1,1},{1,1},{1,1}})

Gaia, 24 bytes
e:w;ċz⟨ọ1>¦ẏ⟩⁇‼⁇E‡ev2%Σ⌠

Try it online!
Ugh, E‡ does some weird stuff...according to the documentation, it should do something like "given length i set of  lists X and length j set of indices Y, return X[i][Y[j]]", but instead returns [X[i][Y[j]] X[i][Y[-j]] where negative indexing represents the complement, so we have to do ev2% to extract only the ones we want.
e				| eval as a list l
 :				| dup
  w				| wrap as a list
   ;				| push l again
    ċ				| push [1..len(l)]
     z				| push all subsets of [1..len(l)] -- index powerset.
      ⟨      ⟩⁇			| filter this for:
       ọ			| deltas
        1>¦			| are greater than 1
           ẏ			| all (all deltas greater than 1)
	       ‼⁇		| filter for non-empty lists
		 E‡		| table extract elements. Given l and index set i, this pushes
				| [l[i] l[setdiff(1..l,i)]] for some reason
		   ev2%		| get the l[i] only by unlisting, reversing, and taking every other element
		       Σ⌠	| Get the one with the maximum sum

Out of curiosity, why does the output have two trailing ]] instead of one?

@KevinCruijssen Just another fun quirk of the interpreter; all lists are printed out like that, so [[1] [2]] gets printed [[1]] [2]]]] which makes it super hard to read/debug list output.

I think it's because of the expression re.sub(" ?$","]",result) in the interpreter which should instead be re.sub(" +$","]",result) but my python is super bad.

Python 2, 63 70 65 bytes
f=lambda a:a and max([a[:1],a[:1]+f(a[2:]),f(a[1:])],key=sum)or a

Try it online!
5 bytes thx to ArBo

Your test case [1, 7, 4, -2] [1, 4] 5 7 is getting the wrong answer.

@xnor: fixed now.

65 bytes

Haskell, 300 168 bytes
import Data.List
h[]=1>2
h(x:y)=fst$foldl(\a c->((fst a)&&(c-snd a>1),c))(1<2,x)y
z=snd.last.sortOn fst.map((,)=<<sum.snd).filter(h.fst).map unzip.subsequences.zip[0..]

Try it online!
-132 bytes thanks to all the feedback from @nimi :)

Original
Ungolfed (original)
import Data.List
import Data.Function

f :: [Int] -> [(Int, Int)] -- attach indices for later use
f [] = []
f xs = zip xs [0..length xs]

g :: [[(Int, Int)]] -> [([Int], [Int])] -- rearrange into list of tuples
g [] = []
g (x:xs) = (map fst x, map snd x) : g xs

h :: [Int] -> Bool -- predicate that checks if the indices are at least 2 apart from each other
h [] = False
h (x:xs) = fst $ foldl (\acc curr -> ((fst acc) && (curr - snd acc > 1), curr)) (True, x) xs
j :: [([Int], [Int])] -> [([Int], [Int])] -- remove sets that don't satisfy the condition
j xs = filter (\(elements, indices) -> h indices) xs

k :: [([Int], [Int])] -> [(Int, ([Int], [Int]))] -- calculate some of elements
k xs = map (\(elements, indices) -> (foldl1 (+) elements, (elements, indices))) xs

l :: [(Int, ([Int], [Int]))] -> ([Int], [Int]) -- grab max
l xs = snd $ last $ sortBy (compare `on` fst) xs

z -- put things together
```

Some tips: flip the element and its index within the pairs returned by f: f x=zip[0..length x]x, so f becomes f=zip[0..]. g is just g=map unzip. The function to filter with in j is h.fst (<- flipped pairs!). j=filter(h.fst). The foldl1+ from k is sum and with a pointfree pair making k=map((,)=<<sum.snd). sortBy(...) can be replaced by sortOn fst: l=snd.last.sortOn fst. Finally as you are using all functions only once, you can inline them into a single pointfree expression: z=snd.last.sortOn fst.map((,)=<<sum.snd).filter(h.fst).map unzip.subsequences.zip[0..]

.... Try it online!.

oh, and no need to import Data.Function anymore.

That's great, thanks for the feedback :)

Next h: we're looking for non-adjacent elements, i.e. the difference of adjacent indices must be >1.  zipWith(-)=<<tail builds such a list of differences, but fails for the empty list, so we need an additional tail on the subsequences to get rid of it. Inline again. Try it online!

Instead of making a pair of the sum and the list and sorting on the first element of the pair, we can directly sort on the sum: sortOn(sum.snd). Try it online!

Charcoal, 46 bytes
≔⟦υ⟧ηＦθ«≔υζ≔Ｅη⁺κ⟦ι⟧υ≔⁺ζηη»≔Φ⁺υηιη≔ＥηΣιζＩ§η⌕ζ⌈ζ

Try it online! Link is to verbose version of code. Explanation:
≔⟦υ⟧η

The variable u is predefined with an empty list. This is put in a list which is assigned to h. These variables act as accumulators. u contains the sublists that include the latest element of the input q while h contains the sublists that do not (and therefore are suitable for appending the next element of the input).
Ｆθ«

Loop over the elements of the input.
≔υζ

Save the list of sublists that contain the previous element.
≔Ｅη⁺κ⟦ι⟧υ

Take all of the sublists that do not contain the previous element, append the current element, and save the result as the list of sublists that contain the current element. (I don't use Push here as I need to clone the list.)
≔⁺ζηη»

Concatenate both previous sublists into the new list of sublists that do not contain the current element.
≔Φ⁺υηιη

Concatenate the sublists one last time and remove the original empty list (which Charcoal can't sum anyway).
≔ＥηΣιζ

Compute the sums of all of the sublists.
Ｉ§η⌕ζ⌈ζ

Find an index of the greatest sum and output the corresponding sublist.

Japt -h, 21 bytes
Ever have one of those challenges where you completely forget how to golf?!
ð¤à fÊk_än ø1Ã®mgUÃñx

Try it
ð¤à fÊk_än ø1Ã®mgUÃñx     :Implicit input of array U
ð                         :Indices of elements that return true when
 ¤                        :  Converted to a base-2 string (to account for 0s)
  à                       :Combinations
    f                     :Filter by
     Ê                    :  Length (to remove the empty combination)
      k_                  :Remove elements that return true
        än                :  Deltas
           ø1             :  Contains 1
             Ã            :End remove
              ®           :Map
               m          :  Map
                gU        :    Index into U
                  Ã       :End map
                   ñ      :Sort by
                    x     :  Sum
                          :Implicit output of last element

Japt -h, 22 bytes
Êo à k_mÄ øZÃm!gU
fÊñx

Try it

source: Stack Exchange

license: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/4.0/

sort: votes


text:
The Nth Gryphon Number
I came up with a series of numbers the other day and decided to check what the OEIS number for it was.  Much to my surprise, the sequence did not appear to be in the OEIS database, so I decided to name the sequence after myself (note that someone else who's a lot smarter than me has probably already come up with this, and if someone finds the actual name of this sequence, please comment and I'll change the question title).  As I couldn't find the sequence anywhere, I decided to name it after myself, hence "Gryphon Numbers".  EDIT: Thanks to @Surb for bringing to my attention the fact that this sequence is equal to OEIS sequence A053696 - 1.
A Gryphon number is a number of the form \$a+a^2+...+a^x\$, where both \$a\$ and \$x\$ are integers greater than or equal to two, and the Gryphon sequence is the set of all Gryphon numbers in ascending order.  If there are multiple ways of forming a Gryphon number (the first example is \$30\$, which is both \$2+2^2+2^3+2^4\$ and \$5+5^2\$) the number is only counted once in the sequence.  The first few Gryphon numbers are: \$6, 12, 14, 20, 30, 39, 42, 56, 62, 72\$.
Your Task:
Write a program or function that receives an integer \$n\$ as input and outputs the \$n\$th Gryphon number.
Input:
An integer between 0 and 10000 (inclusive).  You may treat the sequence as either 0-indexed or 1-indexed, whichever you prefer.  Please state which indexing system you use in your answer to avoid confusion.
Output:
The Gryphon number corresponding to the input.
Test Cases:
Please note that this assumes the sequence is 0-indexed.  If your program assumes a 1-indexed sequence, don't forget to increment all the input numbers.
Input:    Output:
0   --->  6
3   --->  20
4   --->  30
10  --->  84
99  --->  4692
9999 -->  87525380

Scoring:
This is code-golf, so the lowest score in bytes wins.

The Gryphon sequence is A053696 - 1. In other words, A053696 is the increasing sequence of numbers of the form $\color{red}{a^0}+a^1+\ldots+a^x$.

@Surb ah, that's why I couldn't find it.  In that case, I'll put that information in an edit, but keep the rest of the question the way it is since there doesn't seem to be a better name for the sequence.

Jelly, 9 bytes
bṖ’ḅi-µ#Ṫ

A full program which reads a (1-indexed) integer from STDIN and prints the result.
Try it online!
How?
A Gryphon number is a number which is expressible in a base less than itself such that all the digits are ones except the least significant, which is a zero. For example:
\$30=1\times 2^4+1\times 2^3+1\times2^2+1\times2^1+0\times2^0 \rightarrow 30_2 = 11110\$
\$84=1\times 4^3+1\times 4^2+1\times 4^1+0\times 4^0 \rightarrow 84_4 = 1110\$
This program takes n, then starts at v=0 and tests for this property and increments v until it finds n such numbers, then outputs the final one.
To test a base b number it subtracts one from every digit, converts from base v, and then checks if the result is \$-1\$. (Note that b is less than v)
\$30_2 \rightarrow 0\times 30^4+0\times 30^3+0\times30^2+0\times30^1+(-1)\times30^0 = -1\$
\$84_4 \rightarrow 0\times 84^3+0\times 84^2+0\times 84^1+(-1)\times 84^0 = -1\$
bṖ’ḅi-µ#Ṫ - Main Link: no arguments
       #  - set v=0 then count up collecting n=STDIN matches of:
      µ   -  the monadic link -- i.e. f(v):  e.g. v=6
 Ṗ        -    pop (implicit range of v)            [1,2,3,4,5]
b         -    to base (vectorises)                 [[1,1,1,1,1,1],[1,1,0],[2,0],[1,2],[1,1]]
  ’       -    decrement (vectorises)               [[0,0,0,0,0,0],[0,0,-1],[1,-1],[0,1],[0,0]]
   ḅ      -    from base (v) (vectorises)           [0,-1,5,1,0]
     -    -    literal -1                           -1
    i     -    first index of (zero if not found)   2
          - }  e.g. n=11 -> [6,12,14,20,30,39,42,56,62,72,84]
        Ṫ - tail         -> 84
          - implicit print

MATL, 16 13 bytes
:Qtt!^Ys+uSG)

1-based.
Try it online!
Explanation
Consider input n = 3 as an example.
:    % Implicit input: n. Range
     % STACK: [1 2 3]
Q    % Add 1, element-wise
     % STACK: [2 3 4]
tt   % Duplicate twice, transpose
     % STACK: [2 3 4], [2 3 4], [2;
                                 3;
                                 4]
^    % Power, element-wise with broadcast
     % STACK: [2 3 4], [ 4   9  16;
                         8  27  64;
                        16  81 256]
Ys   % Cumulative sum of each column
     % STACK: [2 3 4], [ 4    9  16;
                         12  36  80;
                         28 117 336]
+    % Add, element-wise with broadcast (*)
     % STACK: [ 6  12  20;
               14  39  84
               30 120 340]
u    % Unique elements. Gives a column vector
     % STACK: [  6;
                14;
                30;
                12;
               ···
               340]
S    % Sort
     % STACK: [  6;
                12
                14;
                20;
               ···
               340]
G)   % Push input again, index. This gets the n-th element. Implicit display
     % STACK: 14

The matrix obtained in step (*) contains possibly repeated Gryphon numbers. In particular, it contains n distinct Gryphon numbers in its first row. These are not necessarily the n smallest Gryphon numbers. However, the lower-left entry 2+2^+···+2^n  exceeds the upper-right entry n+n^2, and thus all the numbers in the last row exceed those in the first row. This implies that extending the matrix rightward or downward would not contribute any Gryphon number lower than the lowest n numbers in the matrix. Therefore, the matrix is guaranteed to contain the n smallest Gryphon numbers. Consequently, its n-th lowest unique element is the solution.

What the hell, this is amazing!

Haskell, 53 bytes
([n|n<-[6..],or[a^y+n==n*a+a|a<-[2..n],y<-[3..n]]]!!)

Try it online!
\$\newcommand{\ta}{{a}}\newcommand{\ti}{{i}}\newcommand{\tx}{{x}}\newcommand{\tn}{{n}}\newcommand{\ty}{{y}}\$A number \$\tn\$ is Gryphon if there exist integers \$\ta \geq 2\$ and \$\tx \geq 2\$ such that \$\tn=\sum_{i=1}^\tx \ta^i\$.
We generate an infinite list of all \$\tn \geq 6\$ such that a brute-force search shows this is the case.
The answer is a (zero-indexed) index function into this list, denoted in Haskell as (list!!).
Why is a^y+n==n*a+a correct?
From the formula for summing a geometric progression:

$$\sum_{i=1}^\nu \alpha \rho^{i-1} = \frac{\alpha(1-\rho^\nu)}{1-\rho}$$

we have, letting \$(\alpha, \rho, \nu) = (\ta, \ta, \tx)\$: $$\tn = \sum_{\ti=1}^\tx \ta^\ti = \frac{\ta(1-\ta^\tx)}{1-\ta} = \frac{\ta-\ta^{\tx+1}}{1-\ta}.$$
Rearranging the equation, we get \$n(1-a) = a - a^{x+1}\$.
Rearranging that even further, we get \$\ta^{\tx+1}+n = \tn \ta + \ta\$.
A substitution \$\ty = \tx+1\$ in the brute-force search yields the final expression a^y+n=n*a+a.
Is searching until n enough?

If \$\ta > \tn\$ (in other words, \$\ta \geq \tn+1\$), then $$\ta^\ty + \tn > \ta^2 \geq (\tn+1)\ta = \tn\ta + \ta$$
which proves \$\ta^\ty+\tn \neq \tn\ta+\ta\$. So there's no sense checking any of the values \$\ta > \tn\$.
Similarly: if \$\ty > \tn\$, then $$\ta^\ty + \tn > \ta^\tn = \ta^{\tn-1}\ta \geq 2^{\tn-1}\ta \stackrel{\color{red}*}> (\tn+1)\ta = \tn\ta + \ta,$$
again proving \$\ta^\ty + \tn \neq \tn\ta + \ta\$.
\$\color{red}{{}^{*}}\$We can assume \$2^{n-1}>n+1\$ because we know \$n \geq 6\$, the smallest Gryphon number.

Python 3.8 (Pre-release), 98 92 89 78 71 bytes
lambda n:sorted({a*~-a**x//~-a for a in(r:=range(2,n+3))for x in r})[n]

Try it online!
0-indexed. Integer division must be used here because f(10000) overflows floats.
Generates all Gryphon numbers where \$2 \leq a \leq n+2\$ and \$2 \leq x \leq n+2\$, sorts them, and selects the \$n\$th element.
-6 bytes thanks to Jonathan Allan
-3 bytes thanks to ArBo. I almost did as he suggested myself, but tried to use {*(...)} which didn't save space anyway
-11 bytes thanks to mathmandan
-7 bytes thanks to ArBo
Mathematical Proof of Validity
Using 0-indexing for the sake of this proof, even though mathematical convention is 1-indexed.

Let \$G_n\$ be the \$n\$th Gryphon number
Let \$g(a,x) = a + a^2 + ... + a^x\$ (The Gryphon number from \$a\$ and \$x\$)
Let \$A_n\$ be the set of all Gryphon numbers where \$2 \leq a \leq n+2\$ and \$2 \leq x \leq n+2\$
We know that \$A_0 = \{ g(2,2) \} = \{ 6 \} = \{ G_0 \}\$
\$A_{n+1} = \{ g(a, x), g(a+1, x), g(a, x+1), g(a+1, x+1) | g(a, x) \in A_n \} \$
\$g(a+1, x) < g(a+1, x+1)\$ for all \$a\$ and \$x\$
\$g(a, x+1) < g(a+1, x+1)\$ for all \$a\$ and \$x\$
Therefore \$G_{n+1} \neq g(a+1, x+1)\$ if \$G_n = g(a, x)\$
\$g(a+1, x) < g(a+2, x)\$ for all \$a\$ and \$x\$
\$g(a, x+1) < g(a, x+2)\$ for all \$a\$ and \$x\$
Therefore \$G_{n+1}\$ must either be \$g(a+1, x)\$ or \$g(a, x+1)\$ if \$G_n = g(a, x)\$ since no other possibilities exist.
We can use this information to conclude that \$G_{n+1} \in A_{n+1}\$ if \$G_n \in A_n\$
Since we know that \$G_0 \in A_0\$, we can use this rule to induce that \$G_n \in A_n\$ for all \$n\$
Since this can be applied from \$G_0\$ to \$G_n\$, then \$G_n\$ must be at index \$n\$ of \$A_n\$ if \$A_n\$ is ordered from smallest to largest

f= is unnecessary, and lambda n,r=range: will save 4 more (like so)

You can drop the set() and replace it by a set comprehension to get to 89

Also, you can remove the f= from the TIO link by putting it in the header, as in the TIO of my 89-byter

86 bytes with Python 3.8 and assignment expressions

At the line "Therefore Gn+1≠(a+1,x+1) if  Gn=g(a,x)" is a mistake, it should be Gn+1≠g(a+1,x+1) if ...

Why not use the geometric series formula: sum(a**b for b in r(1,x)) is equivalent to (a-a**x)/(1-a), provided a and x are at least 2.  I'd recommend lambda n,r=range:sorted({(a**x-a)/~-a for a in r(2,n+3)for x in r(3,n+4)})[n] for 77 bytes. (This will return floats, but the sorting should still work fine.  Use the // operation if returning integers is required.)

71 bytes using pre-release Python 3.8, together with mathmandan's formula and some algebra to make the ranges the same. (70 bytes if floats are ok)

J, 35 32 bytes
-3 bytes thanks to FrownyFrog
3 :'y{/:~~.,}.+/\(^~/>:)1+i.y+2'

Try it online!
Explanation is same as original.  Simply uses explicit form to save bytes be removing the multiple @.
original answer, tacit, with explanation: 35 bytes
{[:/:~@~.@,@}.[:+/\@(^~/>:)1+i.@+&2

Try it online!
Similar to Luis Mendo's approach, we create a "power table" (like a times table) with top row 2 3 ... n and left column 1 2 ... n resulting in:
 2   3    4     5     6      7
 4   9   16    25    36     49
 8  27   64   125   216    343
16  81  256   625  1296   2401
32 243 1024  3125  7776  16807
64 729 4096 15625 46656 117649

^~/ >: creates the table, and 1+i.@+&2 creates the 1... n sequences, and we add 2 (+&2) to the input to ensure we always have enough elements to create a table even for 0 or 1 inputs.
After we have the table above the solution is trivial.  We just scan sum the rows +/\, and then remove the first row, flatten, take unique, and sort /:~@~.@,@}..  Finally { uses the original input to index into that result, producing the answer.

explicit notation saves 3

thank you, nice catch.

R, 65 62 bytes
-1 byte thanks to Giuseppe.
n=scan();unique(sort(diffinv(t(outer(2:n,1:n,"^")))[3:n,]))[n]

Try it online!
1-indexed.
Generates a matrix of all values of the form \$a^i\$, takes the cumulative sum, removes the first row (0s) and the second row (entries corresponding to \$x=1\$), then takes the unique sorted values.
Note that sort(unique(...)) would not work, as unique would give unique rows of the matrix, and not unique entries. Using unique(sort(...)) works because sort converts to vector.

It takes a bit more work, but using t and diffinv is 64 bytes

@Giuseppe Thanks! I didn't know diffinv. I golfed down another 2 bytes by replacing [-1:-2,] with [3:n,].

Gaia, 18 bytes
)┅:1D¤*‡+⊣¦tv<_uȯE

Try it online!
1-based index.
This is a rather sad answer with a long nose: )┅: It probably wishes it could be golfed down further.
Copies the algorithm given by Luis Mendo's answer

Wolfram Language (Mathematica), 51 bytes
(Union@@Array[Sum[(#2+1)^k,{k,#+1}]&,{30,#}])[[#]]&

Try it online!
1-indexed
-8 bytes from @attinat

51 bytes

An explanation would be nice.

JavaScript (ES7), 76 bytes
1-indexed.
f=(n,a=[i=2])=>(n-=a.some(j=>a.some(k=>(s+=j**k)==i,s=j)))?f(n,[...a,++i]):i

Try it online!

JavaScript (ES7), 89 bytes
1-indexed.
n=>eval('for(a=[i=1e4];--i>1;)for(s=1e8+i,x=1;a[s+=i**++x]=x<26;);Object.keys(a)[n]-1e8')

Try it online!

05AB1E, 12 bytes
L¦ãε`LmO}êIè

0-indexed
Try it online or verify the first \$n\$ items.
Explanation:
L             # Create a list in the range [1, (implicit) input-integer]
              #  i.e. 4 → [1,2,3,4]
 ¦            # Remove the first item to make the range [2, input-integer]
              #  i.e. [1,2,3,4] → [2,3,4]
  ã           # Create each possible pair of this list by using the cartesian product
              #  i.e. [2,3,4] → [[2,2],[2,3],[2,4],[3,2],[3,3],[3,4],[4,2],[4,3],[4,4]]
   ε          # Map each pair to:
    `         #  Push the values of the pair separately to the stack
              #   i.e. [4,3] → 4 and 3
     L        #  Take the list [1, value] for the second value of the two
              #   i.e. 3 → [1,2,3]
      m       #  And then take the first value to the power of each integer in this list
              #   i.e. 4 and [1,2,3] → [4,16,64]
       O      #  After which we sum the list
              #   i.e. [4,16,64] → 84
        }ê    # After the map: uniquify and sort the values
              #  i.e. [6,14,30,12,39,120,20,84,340] → [6,12,14,20,30,39,84,120,340]
          Iè  # And index the input-integer into it
              #  i.e. [6,12,14,20,30,39,84,120,340] and 4 → 30
              # (after which the result is output implicitly)

Charcoal, 36 bytes
ＮθＦθＦθ⊞υ÷⁻Ｘ⁺²ι⁺³κ⁺²ι⊕ιＦ⊖θ≔Φυ›κ⌊υυＩ⌊υ

Try it online! Link is to verbose version of code. 1-indexed. Uses Luis Mendo's algorithm. Explanation:
Ｎθ

Input \$ n \$.
ＦθＦθ⊞υ

Create an \$ n \$-by-\$ n \$ grid of Gryphon numbers and push each one to the predefined list.
÷⁻Ｘ⁺²ι⁺³κ⁺²ι⊕ι

Calculate the Gryphon number using the fact that \$ \sum_1^x a^i = \frac{a^{x+1}-a}{a-1} \$.
Ｆ⊖θ≔Φυ›κ⌊υυ

Remove the lowest \$ n-1 \$ Gryphon numbers.
Ｉ⌊υ

Print the lowest remaining Gryphon number.

Japt, 23 bytes
Dear Jebus! Either I really have forgotten how to golf or the booze is finally taking its toll!
Not a direct port of Jonathan's solution but very much inspired by his observation.
@ÈÇXìZ mÉ ìZÃeÄ}fXÄ}gNÅ

Try it

all_licenses: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/4.0/


text:
New Order #6: Easter Egg
Introduction (may be ignored)
Putting all positive integers in its regular order (1, 2, 3, ...) is a bit boring, isn't it? So here is a series of challenges around permutations (reshuffelings) of all positive integers. This is the sixth challenge in this series (links to the first, second, third, fourth and fifth challenge).
This challenge has a mild Easter theme (because it's Easter). I took my inspiration from this highly decorated (and in my personal opinion rather ugly) goose egg.

It reminded me of the Ulam spiral, where all positive integers are placed in a counter-clockwise spiral. This spiral has some interesting features related to prime numbers, but that's not relevant for this challenge.

We get to this challenge's permutation of positive integers if we take the numbers in the Ulam spiral and trace all integers in a clockwise turning spiral, starting at 1. This way, we get:
1, 6, 5, 4, 3, 2, 9, 8, 7, 20, 19, 18, 17, 16, 15, 14, 13, 12, 11, 10, 25, 24, 23, etc.

If you would draw both of the spirals, you'd get some sort of an infinite mesh of (egg shell) spirals (note the New Order reference there).
This sequence is present in the OEIS under number A090861. Since this is a "pure sequence" challenge, the task is to output \$a(n)\$ for a given \$n\$ as input, where \$a(n)\$ is A090861.
Task
Given an integer input \$n\$, output \$a(n)\$ in integer format, where \$a(n)\$ is A090861.
Note: 1-based indexing is assumed here; you may use 0-based indexing, so \$a(0) = 1; a(1) = 6\$, etc. Please mention this in your answer if you choose to use this.
Test cases
Input | Output
---------------
1     |  1
5     |  3
20    |  10
50    |  72
78    |  76
123   |  155
1234  |  1324
3000  |  2996
9999  |  9903
29890 |  29796

Rules

Input and output are integers.
Your program should at least support input in the range of 1 up to 32767).
Invalid input (0, floats, strings, negative values, etc.) may lead to unpredicted output, errors or (un)defined behaviour.
Default I/O rules apply.
Default loopholes are forbidden.
This is code-golf, so the shortest answers in bytes wins

Jelly,  16 14 11 10 9  8 bytes
-1 thanks to Lynn (mod-2; logical NOT; add to self: Ḃ¬+ -> bitwise OR with 1:|1)
|1×r)ẎQi

A monadic Link accepting an integer, n, which yields an integer, a(n).
Try it online! (very inefficient since it goes out to layer \$\lceil\frac n2\rceil\$)
An 11-byte version, ½‘|1×rƲ€ẎQi, completes all but the largest test case in under 30s - Try it at TIO - this limits the layers used to \$\lceil\frac{\lfloor\sqrt n\rfloor+1}2\rceil\$.
How?
The permutation is to take the natural numbers in reversed slices of lengths [1,5,3,11,5,17,7,23,9,29,11,35,13,...] - the odd positive integers interspersed with the positive integers congruent to five modulo six, i.e [1, 2*3-1, 3, 4*3-1, 5, 6*3-1, 7, 8*3-1, 9, ...].
This is the same as concatenating and then deduplicating reversed ranges [1..x] of where x is the cumulative sums of these slice lengths (i.e. the maximum of each slice) - [1,6,9,20,25,42,49,72,81,110,121,156,169,...], which is the odd integers squared interspersed with even numbers multiplied by themselves incremented, i.e. [1*1, 2*3, 3*3, 4*5, 5*5, 6*7, 7*7,...].
Since the differences are all greater than 1 we can save a byte (the reversal) by building ranges [x..k] directly, where k is the 1-indexed index of the slice.
Due to this structure the permutation is a self-conjugate permutation, i.e. we know that \$P(n) = v \iff P(v) = n\$, so rather than finding the value at (1-indexed) index n (|1×r)ẎQị@) we can actually get the (1-indexed) index of the first occurrence of n (|1×r)ẎQi).
|1×r)ẎQi - Link: integer, n       e.g. 10
    )    - for each k in [1..n]:  vs = [ 1, 2, 3, 4, 5, 6, 7, 8, 9,10]
|1       -   bitwise-OR (k) with 1     [ 1, 3, 3, 5, 5, 7, 7, 9, 9,11]
  ×      -   multiply (by k)           [ 1, 6, 9,20,25,42,49,72,81,110]
   r     -   inclusive range (to k)    [[1],[6..2],[9..3],[20..4],...,[110..10]]
     Ẏ   - tighten                     [1,6,5,4,3,2,9,8,7,6,5,4,3,20,...,4,......,110,...,10]
      Q  - de-duplicate                [1,6,5,4,3,2,9,8,7,20,...,10,......,110,...82]
       i - first index with value (n)  20

Very nice. And you surpassed the MATL answer!

Tied now... :-)

@LuisMendo I just realised I can do something sneaky here and save one byte :)

@JonathanAllan Aww. That deserves one upvote :-)

9 bytes: |1×r)ẎQị@

@Lynn I am actually just updating to a different 9 byter. Yours will prob make 8!

@Lynn yep, now at 8, thanks!

Sweet! This is a very nice answer :)

JavaScript (ES7),  46 45  41 bytes
0-indexed.
n=>((x=n**.5+1&~1)*2-(n<x*x+x)*4+3)*x+1-n

Try it online!
How?
This is based on the 1-indexed formula used in the example programs of A090861.
\$x_n\$ is the layer index of the spiral, starting with \$0\$ for the center square:
$$x_n=\left\lfloor\frac{\sqrt{n-1}+1}{2}\right\rfloor$$
Try it online!
\$k_n\$ is set to \$6\$ for the bottom part of each layer (including the center square), and to \$-2\$ everywhere else:
$$k_n=\begin{cases}
-2&\text{if }n\le 4{x_n}^2+2x_n\\
6&\text{otherwise}
\end{cases}$$
Try it online!
Then \$a_n\$ is given by:
$$a_n=8{x_n}^2+k_nx_n+2-n$$
Try it online!
Which can be translated into:
n=>8*(x=(n-1)**.5+1>>1)*x+(n<=4*x*x+2*x?-2:6)*x+2-n

Making it 0-indexed saves 5 bytes right away:
n=>8*(x=n**.5+1>>1)*x+(n<4*x*x+2*x?-2:6)*x+1-n

The formula can be further simplified by using:
$${x'}_n=2\times\left\lfloor\frac{\sqrt{n}+1}{2}\right\rfloor$$
which can be expressed as:
x=n**.5+1&~1

leading to:
n=>2*(x=n**.5+1&~1)*x+(n<x*x+x?-1:3)*x+1-n

and finally:
n=>((x=n**.5+1&~1)*2-(n<x*x+x)*4+3)*x+1-n

Wolfram Language (Mathematica), 60 bytes
8(s=⌊(⌊Sqrt[#-1]⌋+1)/2⌋)^2-#+2+If[#<=4s^2+2s,-2,6]s&

Try it online!

MATL, 12 11 bytes
Eq1YL!tPG=)

Try it online!
Very memory-inefficient. Prepending X^k makes it more efficient.

C# (Visual C# Interactive Compiler), 67 bytes
n=>8*(x=(int)Math.Sqrt(--n)+1>>1)*x+(n<4*x*x+2*x?-2:6)*x+1-n;int x;

Try it online!

Python 3.8, 104 74 65 60 57 bytes
lambda n:(-2,6)[n>4*(x:=(n**.5+1)//2)*x+2*x]*x+2+~n+8*x*x

Edit: Thanks to Johnathan Allan for getting it from 74 to 57 bytes!
This solution uses 0-based indexing.

Save 39 avoiding the imports, removing some redundant parentheses, and using > in place of <= and x*x in place of x**2 ...like so: def f(n):x=((n-1)**.5+1)//2;return 8*x**2+(-2,6)[n>4*x*x+2*x]*x+2-n ...TIO

Awesome!  I will incorporate the edits. Made some changes before I saw your comment and got it down to 74 bytes.  Does it matter that yours returns floats?  I assume not...

Float representations of integers should be fine. Save some more using Python 3.8 assignment ...EDIT: make it zero indexed

Very cool. Feel free to make any further edits directly!

Python 3.8 (pre-release), 53 bytes
A direct port of Arnauld's JavaScript answer, go upvote that, and/or J42161217's Mathematica answer, and/or Kapocsi's Python answer :)
lambda n:((x:=int(n**.5+1)&-2)*2-(n<x*x+x)*4+3)*x+1-n

0-indexed.
Try it online!

Befunge, 67 57 bytes
This solution assumes 0-based indexing for the input values.
p&v-*8p00:+1g00:<
0:<@.-\+1*g00+*<|`
0g6*\`!8*2+00g4^>$:0

Try it online!
Explanation
We start by calculating the "radius" at which the input n is found with a loop:
radius = 0
while n > 0
  radius += 1
  n -= radius*8

At the end of the loop, the previous value of n is the offset into the spiral at that radius:
offset = n + radius*8

We can then determine if we're on the top or bottom section of the spiral as follows:
bottom = offset >= radius*6

And once we have all these details, the spiral value is calculated with:
value = ((bottom?10:2) + 4*radius)*radius + 1 - offset

The radius is the only value that we need to store as a "variable", limiting it to a maximum value of 127 in Befunge-93, so this algorithm can handle inputs up to 65024.

Japt, 15 bytes
Port of Jonathan's Jelly solution. 1-indexed.
gUòÈ²+X*v)õÃcÔâ

Try it
gUòÈ²+X*v)õÃcÔâ     :Implicit input of integer U
g                   :Index into
 Uò                 :  Range [0,U]
   È                :  Map each X
    ²               :    Square X
     +X*            :    Add X multiplied by
        v           :    1 if X is divisible by 2, 0 otherwise
         )          :    Group result
          õ         :    Range [1,result]
           Ã        :  End map
            c       :  Flatten
             Ô      :    After reversing each
              â     :  Deduplicate

I just told Jonathan that x+(1-x%2) is x|1 (saving a byte in Jelly), which this answer can also benefit from, I bet.

Husk, 11 bytes
!uṁS`…S*v1N

Try it online!
Instead of creating a range, an infinite list is created and nubbed. Hence, will be be very very slow for larger \$n\$

C++ (gcc), 88 bytes
#import<cmath>
int a(int n){int x=(sqrt(n-1)+1)/2;return x*(8*(x+(n>4*x*x+2*x))-2)+2-n;}

1-indexed; uses the formula on the OEIS page, but manipulated to save a few bytes.
Try it online!

Suggest sqrt(n-1)/2+.5 instead of (sqrt(n-1)+1)/2

Vyxal, 13 bytes
ƛ1⋎*›rṘ;fU?‹i

Try it Online!
Messy port