text:
Determinant of an Integer Matrix
Given a square integer matrix as input, output the determinant of the matrix.
Rules

You may assume that all elements in the matrix, the determinant of the matrix, and the total number of elements in the matrix are within the representable range of integers for your language.
Outputting a decimal/float value with a fractional part of 0 is allowed (e.g. 42.0 instead of 42).
Builtins are allowed, but you are encouraged to include a solution that does not use builtins.

Test Cases
[[42]] -> 42
[[2, 3], [1, 4]] -> 5
[[1, 2, 3], [4, 5, 6], [7, 8, 9]] -> 0
[[13, 17, 24], [19, 1, 3], [-5, 4, 0]] -> 1533
[[372, -152, 244], [-97, -191, 185], [-53, -397, -126]] -> 46548380
[[100, -200, 58, 4], [1, -90, -55, -165], [-67, -83, 239, 182], [238, -283, 384, 392]] -> 571026450
[[432, 45, 330, 284, 276], [-492, 497, 133, -289, -28], [-443, -400, 56, 150, -316], [-344, 316, 92, 205, 104], [277, 307, -464, 244, -422]] -> -51446016699154
[[416, 66, 340, 250, -436, -146], [-464, 68, 104, 471, -335, -442], [159, -407, 310, -489, -248, 370], [62, 277, 446, -325, 47, -193], [460, 460, -418, -28, 234, -374], [249, 375, 489, 172, -423, 125]] -> 39153009069988024
[[-246, -142, 378, -156, -373, 444], [186, 186, -23, 50, 349, -413], [216, 1, -418, 38, 47, -192], [109, 345, -356, -296, -47, -498], [-283, 91, 258, 66, -127, 79], [218, 465, -420, -326, -445, 19]] -> -925012040475554
[[-192, 141, -349, 447, -403, -21, 34], [260, -307, -333, -373, -324, 144, -190], [301, 277, 25, 8, -177, 180, 405], [-406, -9, -318, 337, -118, 44, -123], [-207, 33, -189, -229, -196, 58, -491], [-426, 48, -24, 72, -250, 160, 359], [-208, 120, -385, 251, 322, -349, -448]] -> -4248003140052269106
[[80, 159, 362, -30, -24, -493, 410, 249, -11, -109], [-110, -123, -461, -34, -266, 199, -437, 445, 498, 96], [175, -405, 432, -7, 157, 169, 336, -276, 337, -200], [-106, -379, -157, -199, 123, -172, 141, 329, 158, 309], [-316, -239, 327, -29, -482, 294, -86, -326, 490, -295], [64, -201, -155, 238, 131, 182, -487, -462, -312, 196], [-297, -75, -206, 471, -94, -46, -378, 334, 407, -97], [-140, -137, 297, -372, 228, 318, 251, -93, 117, 286], [-95, -300, -419, 41, -140, -205, 29, -481, -372, -49], [-140, -281, -88, -13, -128, -264, 165, 261, -469, -62]] -> 297434936630444226910432057

Sandbox, relevant meta post

is there a maximum size of matrix that needs to be supported or is it arbitrary?

@TaylorScott First rule listed: You may assume that all elements in the matrix, the determinant of the matrix, and the total number of elements in the matrix are within the representable range of integers for your language.

You know you have an interesting challenge when you have 4 Jelly answers consecutively out-golfing each other...

Jelly, 15 bytes
LŒ!ðŒcIṠ;ị"Pð€S

Try it online!
How it works
LŒ!ðŒcIṠ;ị"Pð€S   input
L                 length
 Œ!               all_permutations
   ð        ð€    for each permutation:
    Œc                take all unordered pairs
      I               calculate the difference between
                      the two integers of each pair
       Ṡ              signum of each difference
                      (positive -> 1, negative -> -1)
        ;             append:
         ị"             the list of elements generated by taking
                        each row according to the index specified
                        by each entry of the permutation
           P          product of everything
              S   sum

Why it works -- mathy version
The operator det takes a matrix and returns a scalar. An n-by-n matrix can be thought of as a collection of n vectors of length n, so det is really a function that takes n vectors from ℤn and returns a scalar.
Therefore, I write det(v1, v2, v3, ..., vn) for det [v1 v2 v3 ... vn].
Notice that det is linear in each argument, i.e. det(v1 + λw1, v2, v3, ..., vn) = det(v1, v2, v3, ..., vn) + λ det(w1, v2, v3, ..., vn). Therefore, it is a linear map from (ℤn)⊗n to ℤ.
It suffices to determine the image of the basis under the linear map. The basis of (ℤn)⊗n consists of n-fold tensor products of the basis elements of ℤn, i.e. e.g. e5 ⊗ e3 ⊗ e1 ⊗ e5 ⊗ e1. Tensor products that include identical tensors must be sent to zero, since the determinant of a matrix in which two columns are identical is zero. It remains to check what the tensor products of distinct basis elements are sent to. The indices of the vectors in the tensor product form a bijection, i.e. a permutation, in which even permutations are sent to 1 and odd permutations are sent to -1.
For example, to find the determinant of [[1, 2], [3, 4]]: note that the columns are [1, 3] and [2, 4]. We decompose [1, 3] to give (1 e1 + 3 e2) and (2 e1 + 4 e2). The corresponding element in the tensor product is (1 e1 ⊗ 2 e1 + 1 e1 ⊗ 4 e2 + 3 e2 ⊗ 2 e1 + 3 e2 ⊗ 4 e2), which we simplify to (2 e1 ⊗ e1 + 4 e1 ⊗ e2 + 6 e2 ⊗ e1 + 12 e2 ⊗ e2). Therefore:
det [[1, 2], [3, 4]]
= det(1 e1 + 3 e2, 2 e1 + 4 e2)
= det(2 e1 ⊗ e1 + 4 e1 ⊗ e2 + 6 e2 ⊗ e1 + 12 e2 ⊗ e2)
= det(2 e1 ⊗ e1) + det(4 e1 ⊗ e2) + det(6 e2 ⊗ e1) + det(12 e2 ⊗ e2)
= 2 det(e1 ⊗ e1) + 4 det(e1 ⊗ e2) + 6 det(e2 ⊗ e1) + 12 det( e2 ⊗ e2)
= 2 (0) + 4 (1) + 6 (-1) + 12 (0)
= 4 - 6
= -2
Now it remains to prove that the formula for finding the parity of the permutation is valid. What my code does is essentially find the number of inversions, i.e. the places where an element on the left is bigger than an element on the right (not necessarily consecutively).
For example, in the permutation 3614572, there are 9 inversions (31, 32, 61, 64, 65, 62, 42, 52, 72), so the permutation is odd.
The justification is that each transposition (swapping two elements) either adds one inversion or takes away one inversion, swapping the parity of the number of inversions, and the parity of the permutation is the parity of the number of transpositions needed to achieve the permutation.
Therefore, in conclusion, our formula is given by:
Why it works -- non-mathy version

where σ is a permutation of n the group of all permutations on n letters, and sgn is the sign of the permutation, AKA (-1) raised to the parity of the permutation, and aij is the (ij)th entry in the matrix (i down, j across).

That "non-mathy version" is pretty damn mathy.

@MDXF formulas and symbols and numbers hardly constitute maths. Mathematics is the abstraction and the generalization and the logic behind the formal manipulations of the symbols.

"15 bytes" In what encoding? UTF-8 requires 25 bytes for that set of characters. (I realize this has been discussed on meta, but it still seems weird to say bytes instead of characters when the byte count could well be variable.)

@JAB Jelly implements its own custom code page. (One of these days, TIO will include a link to the code page...)

So essentially you're using abstract algebra to construct the sum-of-diagonal-products formulas.

@Mego the "sum-of-diagonal-products" only work for 1x1, 2x2, and 3x3 matrices. For bigger matrices, you need to consider all permutations and their parity.

@LeakyNun True, but for the 2x2 and 3x3 cases, you end up deriving the simple formulas.

+1 for actually include the proof in the post instead of saying "because this formula is listed on the page abcxyz it must be true".

Jelly, 16 15 12 10 bytes
Ḣ×Zß-Ƥ$Ṛḅ-

Uses Laplace expansion. Thanks to @miles for golfing off 3 5 bytes!
Try it online!
How it works
Ḣ×Zß-Ƥ$Ṛḅ-  Main link. Argument: M (matrix / 2D array)

Ḣ           Head; pop and yield the first row of M.
      $     Combine the two links to the left into a monadic chain.
  Z         Zip/transpose the matrix (M without its first row).
   ß-Ƥ      Recursively map the main link over all outfixes of length 1, i.e., over
            the transpose without each of its rows.
            This yields an empty array if M = [[x]].
 ×          Take the elementwise product of the first row and the result on the
            right hand. Due to Jelly's vectorizer, [x] × [] yields [x].
       Ṛ    Reverse the order of the products.
        ḅ-  Convert from base -1 to integer.
                [a]          -> (-1)**0*a
                [a, b]       -> (-1)**1*a + (-1)**0*b = b - a
                [a, b, c]    -> (-1)**2*a + (-1)**1*b + (-1)**0*c = c - b + a
                etc.

R, 3 bytes
Trivial Solution
det

Try it online!
R, 94 92 89 bytes
re-implemented solution
outgolfed by Jarko Dubbeldam

d=function(m)"if"(x<-nrow(m),m[,1]%*%sapply(1:x,function(y)(-1)^y*-d(m[-y,-1,drop=F])),1)

Try it online!
Recursively uses expansion by minors down the first column of the matrix.
f <- function(m){
 x <- nrow(m)                 # number of rows of the matrix
 if(sum(x) > 1){              # when the recursion reaches a 1x1, it has 0 rows
                              # b/c [] drops attributes
  minor <- function(y){
   m[y] * (-1)^y *
   -d(m[-y,-1])                # recurse with the yth row and first column dropped
   }
  minors <- sapply(1:x,minor) # call on each row
  sum(minors)                 # return the sum
 } else {
  m                           # return the 1x1 matrix
 }
}

32 bytes

Haskell, 71 bytes
-3 bytes thanks to Lynn. Another one bytes the dust thanks to Craig Roy.
f[]=1
f(h:t)=foldr1(-)[v*f[take i l++drop(i+1)l|l<-t]|(i,v)<-zip[0..]h]

Try it online! Added -O flag for optimization purposes. It is not necessary.
Explanation (outdated)
f recursively implements cofactor expansion.
f[[x]]=x

This line covers the base case of a 1 × 1 matrix, in which case the determinant is mat[0, 0].
f(h:t)=

This uses Haskell's pattern matching to break the matrix into a head (the first row) and a tail (the rest of the matrix).
          [                                     |(i,v)<-zip[0..]h]

Enumerate the head of the matrix (by zipping the infinite list of whole numbers and the head) and iterate over it.
           (-1)*i*v

Negate the result based on whether its index is even since the calculation of the determinant involves alternating addition and subtraction.
                     [take i l++drop(i+1)l|l<-t]

This essentially removes the ith column of the tail by taking i elements and concatenating it with the row with the first (i + 1)th elements dropped for every row in the tail.
                   *f

Calculate the determinant of the result above and multiply it with the result of (-1)*i*v.
       sum

Sum the result of the list above and return it.

you could save 1 byte if you replace the sum[(-1)^i*... with foldr(-)0[...

Wolfram Language (Mathematica), between 14 and 42 bytes
We've had a 3-byte built-in and a 53-byte solution that completely avoids built-ins, so here are some weirder solutions somewhere in between.
The Wolfram Language has a lot of very intense functions for decomposing matrices into products of other matrices with simpler structure. One of the simpler ones (meaning I've heard of it before) is Jordan decomposition. Every matrix is similar to a (possibly complex-valued) upper triangular matrix made of diagonal blocks with a specific structure, called the Jordan decomposition of that matrix. Similarity preserves determinants, and the determinant of a triangular matrix is the product of the diagonal elements, so we can compute the determinant with the following 42 bytes:
1##&@@Diagonal@Last@JordanDecomposition@#&

The determinant is also equal to the product of the eigenvalues of a matrix, with multiplicity. Luckily, Wolfram's eigenvalue function keeps track of multiplicity (even for non-diagonalisable matrices), so we get the following 20 byte solution:
1##&@@Eigenvalues@#&

The next solution is kind of cheating and I'm not really sure why it works. The Wronskian of a list of n functions is the determinant of the matrix of the first n-1 derivatives of the functions. If we give the Wronskian function a matrix of integers and say that the variable of differentiation is 1, somehow it spits out the determinant of the matrix. It's weird, but it doesn't involve the letters "Det" and it's only 14 bytes…
#~Wronskian~1&

(The Casoratian determinant works as well, for 1 more byte: #~Casoratian~1&)
In the realm of abstract algebra, the determinant of an n x n matrix (thought of as the map k → k that is multiplication by the determinant) is the nth exterior power of the matrix (after picking an isomorphism k → ⋀n kn). In Wolfram language, we can do this with the following 26 bytes:
HodgeDual[TensorWedge@@#]&

And here's a solution that works for positive determinants only. If we take an n-dimensional unit hypercube and apply a linear transformation to it, the n-dimensional "volume" of the resulting region is the absolute value of the determinant of the transformation. Applying a linear transformation to a cube gives a parallelepiped, and we can take its volume with the following 39 bytes of code:
RegionMeasure@Parallelepiped[Last@#,#]&

The solution I had along these lines was Exp@*Tr@*MatrixLog, but unfortunately this doesn't work for singular matrices.

@MishaLavrov Ooh, that's clever! I think you can fix it with Check[E^Tr@MatrixLog@#,0]&.

That's great! I hadn't been aware of Check before.

I created a challenge for Jordan Decomposition a while ago. You might be interested in it, too. What a great answer!

Haskell, 59 bytes
p%(l:r)=l!!0*f(tail<$>p++r)-(p++[l])%r
[]%_=1
p%_=0
f=([]%)

Try it online!

Proton, 98 bytes
f=m=>(l=len(m))-1?sum((-1)**i*m[0][i]*f([[m[k][j]for k:1..l]for j:0..l if j-i])for i:0..l):m[0][0]

Try it online!
-3 bytes thanks to Mr. Xcoder
-3 bytes thanks to Erik the Outgolfer
-1 byte thanks to S.S. Anne
Expansion over the first row

Just because Proton doesn't have builtin for determinant.

103 Bytes. ((~i%2)*2-1) -> ((-i%2)|1)

Also 102 bytes by replacing j!=i with j-i or i-j.

99 bytes

Funny that the first variables in the program are f, m l :P

@Challenger5 that was um totally an intentional hidden joke thank you very much for noticing um. :P that's interesting I guess

A function operating on a matrix of lists... very fml indeed. :P

@totallyhuman -1 it's not a matrix of lists xD :P

98

R, 32 bytes
function(m)Re(prod(eigen(m)$va))

Uses Not a Tree's algorithm of taking the eigenvalues of the matrix and taking the real part of their product.
Try it online!

C,  176  125 bytes
Thanks to @ceilingcat for golfing 42 bytes, and thanks to both @Lynn and @Jonathan Frech for saving a byte each!
d(M,n)int*M;{int i=n--,s=*M*!n,c,T[n*n];for(;i--;s+=M[i]*(1-i%2*2)*d(T,n))for(c=n*n;c--;T[c]=M[n-~c+c/n+(c%n>=i)]);return s;}

Calculates the determinant using the Laplace expansion along the first row.
Try it online!
Unrolled:
d(M, n)int*M;
{
    int i=n--, s=*M*!n, c, T[n*n];
    for (; i--; s+=M[i]*(1-i%2*2)*d(T,n))
        for (c=n*n; c--;)
            T[c] = M[n-~c+c/n+(c%n>=i)];
    return s;
}

(i%2*-2+1) → (1-i%2*2) saves one more byte.

n+1+c can be n-~c.

Suggest i=s instead of return s

Jelly, 43 bytes
Finally I've done writing my non-builtin solution in a golfing language!
ḣ⁹’’¤;ṫḊ€Ç×⁸ị@⁹’¤¤ḷ/¤
çÐ€⁸J‘¤µJ-*×NS
ÇḢḢ$Ṗ?

Thanks to HyperNeutrino for saving a byte!
Try it online! (spaced code for clarity)
terribly long way to remove n'th elements from a list, will improve later

This answer had been outgolfed by answers of HyperNeutrino, Dennis and Leaky Nun. Jelly is very popular as a golfing language.
Quick explanation:
ÇḢḢ$Ṗ?    Main link.
     ?    If
    Ṗ     after remove the last element, the value is not empty (truthy)
Ç         then execute the last link
 ḢḢ$      else get the element at index [1, 1].

çÐ€⁸J‘¤µJ-*×NS     Helper link 1, take input as a matrix.
çÐ€                Apply the previous link, thread right argument to
   ⁸J‘¤            the range [2, 3, ..., n+1]
       µ           With the result,
        J-*        generate the range [-1, 1, -1, 1, ...] with that length
           ×N      Multiply by negative
             S     Sum

ḣ⁹’’¤;ṫḊ€Ç×⁸ị@⁹’¤¤ḷ/¤    Helper link 2, take left input as a matrix, right input as a number in range [2..n+1]
ḣ
 ⁹’’¤                    Take head ρ-2 of the matrix
     ;                   concatenate with
      ṫ                  tail ρ (that is, remove item ρ-1)
       Ḋ€                Remove first column
         Ç               Calculate determinant of remaining matrix
          ×         ¤    multiply by
                  ḷ/     the first column,
            ị@           row #
              ⁹’¤        ρ-1 (just removed in determinant calculation routine) of
           ⁸     ¤       the matrix.

Jelly, 24 bytes
œcL’$ṚÑ€
J-*×Ḣ€×ÇSµḢḢ$Ṗ?

Try it online!
Explanation
œcL’$ṚÑ€         Helper Link; get the next level of subdeterminants (for Laplace Expansion)
œc               Combinations without replacement of length:
  L’$            Length of input - 1 (this will get all submatrices, except it's reversed)
     Ṛ           Reverse the whole thing
      Ñ€         Get the determinant of each of these
J-*×Ḣ€×ÇSµḢḢ$Ṗ?  Main Link
              ?  If the next value is truthy
             Ṗ   Remove the last element (truthy if the matrix was at least size 2)
J-*×Ḣ€×ÇSµ       Then expand
          ḢḢ$    Otherwise, get the first element of the first element (m[0][0] in Python)
J                [1, 2, ..., len(z)]
 -*              (-1 ** z) for each z in the length range
   ×             Vectorizing multiply with
    Ḣ€           The first element of each (this gets the first column); modifies each row (makes it golfier yay)
      ×Ç         Vectorizing multiply with the subdeterminants
        S        Sum

-2 bytes thanks to user202729's solution

CJam (50 45 bytes)
{:A_,{1$_,,.=1b\)/:CAff*A@zf{\f.*1fb}..-}/;C}

This is an anonymous block (function) which takes a 2D array on the stack and leaves an integer on the stack.
Online test suite
Dissection
This implements the Faddeev-LeVerrier algorithm, and I think it's the first answer to take that approach.

The objective is to calculate the coefficients \$c_k\$ of the characteristic polynomial of the \$n\times n\$ matrix \$A\$,
  $$p(\lambda )\equiv \det(\lambda I_{n}-A)=\sum _{k=0}^{n}c_{k}\lambda ^{k}$$
  where, evidently, \$c_n = 1\$ and \$c_0 = (−1)^n \det A\$.
The coefficients are determined recursively from the top down, by dint of the auxiliary matrices \$M\$,
  $$\begin{aligned}M_{0}&\equiv 0&c_{n}&=1\qquad &(k=0)\\M_{k}&\equiv AM_{k-1}+c_{n-k+1}I\qquad \qquad &c_{n-k}&=-{\frac {1}{k}}\mathrm {tr} (AM_{k})\qquad &k=1,\ldots ,n~.\end{aligned}$$

The code never works directly with \$c_{n-k}\$ and \$M_k\$, but always with  \$(-1)^k c_{n-k}\$ and \$(-1)^{k+1}AM_k\$, so the recurrence is
$$(-1)^k c_{n-k} = \frac1k \mathrm{tr} ((-1)^{k+1} AM_{k}) \\
(-1)^{k+2} AM_{k+1} = (-1)^k c_{n-k}A - A((-1)^{k+1}A M_k)$$
{               e# Define a block
  :A            e#   Store the input matrix in A
  _,            e#   Take the length of a copy
  {             e#     for i = 0 to n-1
                e#       Stack: (-1)^{i+2}AM_{i+1} i
    1$_,,.=1b   e#       Calculate tr((-1)^{i+2}AM_{i+1})
    \)/:C       e#       Divide by (i+1) and store in C
    Aff*        e#       Multiply by A
    A@          e#       Push a copy of A, bring (-1)^{i+2}AM_{i+1} to the top
    zf{\f.*1fb} e#       Matrix multiplication
    ..-         e#       Matrix subtraction
  }/
  ;             e#   Pop (-1)^{n+2}AM_{n+1} (which incidentally is 0)
  C             e#   Fetch the last stored value of C
}

Octave, 28 bytes
@(x)round(prod(diag(qr(x))))

Try it online!
This uses the QR decomposition of a matrix X into an orthgonal matrix Q and an upper triangular matrix R. The determinant of X  is the product of those of Q and R. An orthogonal matrix has unit determinant, and for a triangular matrix the determinant is the product of its diagonal entries. Octave's qr function called with a single output gives R.
The result is rounded to the closest integer. For large input matrices, floating-point inaccuracies may produce an error exceeding 0.5 and thus produce a wrong result.

That's an interesting way to evade the det builtin. ;)

@tomsmeding :-) Also, it had already been "used" in Stewie's answer

MATL, 3 bytes / 5 bytes
With built-in function
&0|

Try it online!
Without built-in
Thanks to Misha Lavrov for pointing out a mistake, now corrected
YvpYo

Try it online!
This computes the determinant as the product of the eigenvalues, rounded to the closest integer to avoid floating-point inaccuracies.
Yv       % Implicit input. Push vector containing the eigenvalues
p        % Product
Yo       % Round. Implicit display

Wouldn't the product of the singular values only tell you the absolute value of the determinant?

@MishaLavrov You're totally right! Thanks for noticing. I corrected it by using eigenvalues instead of singular values... and that saved 4 bytes \o/

Wolfram Language (Mathematica), 53 52 bytes
1##&@@@(t=Tuples)@#.Signature/@t[Range@Tr[1^#]&/@#]&

Try it online!
Unfortunately, computing the determinant of an n by n matrix this way uses O(nn) memory, which puts large test cases out of reach.
How it works
The first part, 1##&@@@(t=Tuples)@#, computes all possible products of a term from each row of the given matrix. t[Range@Tr[1^#]&/@#] gives a list of the same length whose elements are things like {3,2,1} or {2,2,3} saying which entry of each row we picked out for the corresponding product.
We apply Signature to the second list, which maps even permutations to 1, odd permutations to -1, and non-permutations to 0. This is precisely the coefficient with which the corresponding product appears in the determinant.
Finally, we take the dot product of the two lists.

If even Signature is too much of a built-in, at 73 bytes we can take
1##&@@@(t=Tuples)@#.(1##&@@Order@@@#~Subsets~{2}&/@t[Range@Tr[1^#]&/@#])&

replacing it by 1##&@@Order@@@#~Subsets~{2}&. This computes Signature of a possibly-permutation by taking the product of Order applied to all pairs of elements of the permutation. Order will give 1 if the pair is in ascending order, -1 if it's in descending order, and 0 if they're equal.

-1 byte thanks to @user202729

52 bytes (in case you didn't know this Mathematica golfing tip)

I did, but somehow forgot about it here. Thanks!

Haskell, 62 bytes
a#((b:c):r)=b*d(a++map tail r)-(a++[c])#r
_#_=0
d[]=1
d l=[]#l

Try it online! (Footer with test cases taken from @totallyhuman's solution.)
d computes the determinant using a Laplace expansion along the first column. Needs three bytes more than the permanent.

Python 2, 75 bytes
f=lambda m,p=[]:m[0][0]*f(zip(*p+m[1:])[1:])-f(m[1:],p+m[:1])if m else[]==p

Try it online!

Java (OpenJDK 8), 195 192 177 176 bytes
long d(int[][]m){long D=0;for(int i=m.length,l=i-1,t[][]=new int[l][l],j,k;i-->0;D+=m[0][i]*(1-i%2*2)*(l<1?1:d(t)))for(j=l*l;j-->0;)t[j/l][k=j%l]=m[1+j/l][k<i?k:k+1];return D;}

Try it online!
Like many other answers, this also uses the Laplace formula. A slightly less golfed version:

long d(int[][]m){
  long D=0;
  int i=m.length,l=i-1,t[][]=new int[l][l],j,k
  for(;i-->0;)
    for(j=l*l;j-->0;)
      t[j/l][k=j%l]=m[1+j/l][k<i?k:k+1];
    D+=m[0][i]*(1-i%2*2)*(l<1?1:d(t));
  return D;
}

Octave, 30 bytes
@(x)-prod(diag([~,l,~]=lu(x)))

Try it online!
or, the boring 4 byte solution (saved 6 bytes thanks to Luis Mendo (forgot the rules regaring builtin functions)):
@det

Explanation:
Coming up! :)

TI-Basic, 2 bytes
det(Ans

Ah, well.
Please don't upvote trivial answers.
As a high school student (who's forced to own one of these calculators), this function is hella useful so...

It's still hella useful in college - linear algebra does not go away

@TaylorScott In fact, it comes back with a vengeance in differential equations.

@Mego - you are right on that; though for some reason they let me take all of calc and that before linear :/

@TaylorScott Due to an oversight by my university's math department, linalg was not a prerequisite for diffeq when I took it. When my professor realized that, he quickly gave us a 3-day crash course in linalg.

Python 2, 95 bytes
-12 bytes thanks to Lynn.
Port of my Haskell answer.
f=lambda m:sum((-1)**i*v*f([j[:i]+j[i+1:]for j in m[1:]])for i,v in enumerate(m[0]))if m else 1

Try it online!

Here, too, you can use [] as a base case: f=lambda m:sum((-1)**i*v*f([j[:i]+j[i+1:]for j in m[1:]])for i,v in enumerate(m[0]))if m else 1 for 95 bytes!

m==[]or sum(...) gives 92 bytes.

Java 8, 266 261 259 258 bytes
long d(int[][]m){long r=0;int i=0,j,k,l=m.length,t[][]=new int[l-1][l-1],q=m[0][0];if(l<3)return l<2?q:q*m[1][1]-m[0][1]*m[1][0];for(;i<l;r+=m[0][i]*(1-i++%2*2)*d(t))for(j=0;++j<l;)for(k=l;k-->0;){q=m[j][k];if(k<i)t[j-1][k]=q;if(k>i)t[j-1][k-1]=q;}return r;}

Look mom, no build-ins.. because Java has none.. >.>
-7 bytes thanks to @ceilingcat.
Explanation:
Try it here. (Only the last test case is too big to fit in a long of size 263-1.)
long d(int[][]m){             // Method with integer-matrix parameter and long return-type
  long r=0;                   //  Return-long, starting at 0
  int i=0,j,k,                //  Index-integers
      l=m.length,             //  Dimensions of the square matrix
      t[][]=new int[l-1][l-1],//  Temp-matrix, one size smaller than `m`
      q=m[0][0];              //  The first value in the matrix (to reduce bytes)
  if(l<3)                     //  If the dimensions are 1 or 2:
    return l<2?               //   If the dimensions are 1:
      q                       //    Simply return the only item in it
     :                        //   Else (the dimensions are 2):
      q*m[1][1]-m[0][1]*m[1][0];
                              //    Calculate the determinant of the 2x2 matrix
                              //  If the dimensions are 3 or larger:
  for(;i<l;                   //  Loop (1) from 0 to `l` (exclusive)
      r+=                     //    After every iteration: add the following to the result:
         m[0][i]              //     The item in the first row and `i`'th column,
         *(1-i++%2*2)         //     multiplied by 1 if `i` is even; -1 if odd,
         *d(t))               //     multiplied by a recursive call with the temp-matrix
    for(j=0;                  //   Reset index `j` to 0
        ++j<l;)               //   Inner loop (2) from 0 to `l` (exclusive)
      for(k=l;k-->0;){        //    Inner loop (3) from `l-1` to 0 (inclusive)
        q=m[j][k];            //     Set the integer at location `j,k` to reduce bytes
        if(k<i)               //     If `k` is smaller than `i`:
          t[j-1][k]=q;        //      Set this integer at location `j-1,k`
        if(k>i)               //     Else-if `k` is larger than `i`:
          t[j-1][k-1]=q;      //      Set this integer at location `j-1,k-1`
                              //     Else: `k` and `i` are equals: do nothing (implicit)
      }                       //    End of inner loop (3)
                              //   End of inner loop (2) (implicit / single-line body)
                              //  End of loop (1) (implicit / single-line body)
  return r;                   //  Return the result-long
}                             // End of method

JavaScript (ES6), 91
Recursive Laplace
q=(a,s=1)=>+a||a.reduce((v,[r],i)=>v-(s=-s)*r*q(a.map(r=>r.slice(1)).filter((r,j)=>j-i)),0)

Less golfed
q = (a,s=1) => // s used as a local variable
  a[1] // check if a is a single element array
       // if array, recursive call expanding along 1st column
  ? a.reduce((v,[r],i) => v-(s=-s)*r*q(a.map(r=>r.slice(1)).filter((r,j)=>j-i)),0)
  : +a // single element, convert to number

Test

q=(a,s=1)=>+a||a.reduce((v,[r],i)=>v-(s=-s)*r*q(a.map(r=>r.slice(1)).filter((r,j)=>j-i)),0)

TestCases=`[[42]] -> 42
[[2, 3], [1, 4]] -> 5
[[1, 2, 3], [4, 5, 6], [7, 8, 9]] -> 0
[[13, 17, 24], [19, 1, 3], [-5, 4, 0]] -> 1533
[[372, -152, 244], [-97, -191, 185], [-53, -397, -126]] -> 46548380
[[100, -200, 58, 4], [1, -90, -55, -165], [-67, -83, 239, 182], [238, -283, 384, 392]] -> 571026450
[[432, 45, 330, 284, 276], [-492, 497, 133, -289, -28], [-443, -400, 56, 150, -316], [-344, 316, 92, 205, 104], [277, 307, -464, 244, -422]] -> -51446016699154`
.split('\n')

TestCases.forEach(r=>{
  [a,k] = r.split (' -> ')
  a = eval(a)
  d = q(a)
  console.log('Test '+(k==d ? 'OK':'KO')+
    '\nMatrix '+a.join('|')+
    '\nResult '+d+
    '\nCheck  '+k)
})

83 bytes with the same behavior

Or 85 bytes to support the empty matrix (whose determinant should be 1).

(I used the same optimizations in this answer, which is derived from yours.)

SageMath, various
Here are a bunch of methods for computing the determinant that I found interesting, all programmed in SageMath. They can all be tried here.
Builtin, 3 bytes
det

This one isn't too interesting. Sage provides global-level aliases to many common operations that would normally be object methods, so this is shorter than lambda m:m.det().

Real Part of Product of Eigenvalues, 36 bytes
lambda m:real(prod(m.eigenvalues()))

Unfortunately, eigenvalues is not one of those global-level aliases. That, combined with the fact that Sage doesn't have a neat way to compose functions, means we're stuck with a costly lambda. This function symbolic values which are automatically converted to numeric values when printed, so some floating point inaccuracy may be present in some outputs.

Product of Diagonal in Jordan Normal Form, 60 bytes
lambda m:prod(m.jordan_form()[x,x]for x in range(m.nrows()))

In Jordan Normal form, an NxN matrix is represented as a block matrix, with N blocks on the diagonal. Each block consists of either a single eigenvalue, or a MxM matrix with a repeated eigenvalue on the diagonal and 1s on the super-diagonal (the diagonal above and to the right of the "main" diagonal). This results in a matrix with all eigenvalues (with multiplicity) on the main diagonal, and some 1s on the super-diagonal corresponding to repeated eigenvalues. This returns the product of the diagonal of the Jordan normal form, which is the product of the eigenvalues (with multiplicty), so this is a more roundabout way of performing the same computation as the previous solution.
Because Sage wants the Jordan normal form to be over the same ring as the original matrix, this only works if all of the eigenvalues are rational. Complex eigenvalues result in an error (unless the original matrix is over the ring CDF (complex double floats) or SR). However, this means that taking the real part is not necessary, compared to the above solution.

Product of Diagonal in Smith Decomposition
lambda m:prod(m.smith_form()[0].diagonal())

Unlike Jordan normal form, Smith normal form is guaranteed to be over the same field as the original matrix. Rather than computing the eigenvalues and representing them with a block diagonal matrix, Smith decomposition computes the elementary divisors of the matrix (which is a topic a bit too complicated for this post), puts them into a diagonal matrix D, and computes two matrices with unit determinant U and V such that D = U*A*V (where A is the original matrix). Since the determinant of the product of matrices equals the product of the determinants of the matrices (det(A*B*...) = det(A)*det(B)*...), and U and V are defined to have unit determinants, det(D) = det(A). The determinant of a diagonal matrix is simply the product of the elements on the diagonal.
Laplace Expansion, 109 bytes
lambda m:m.nrows()>1and sum((-1)**j*m[0,j]*L(m[1:,:j].augment(m[1:,j+1:]))for j in range(m.ncols()))or m[0,0]

This performs Laplace expansion along the first row, using a recursive approach. det([[a]]) = a is used for the base case. It should be shorter to use det([[]]) = 1 for the base case, but my attempt at that implementation had a bug that I haven't been able to track down yet.

Leibniz's Formula, 100 bytes
L2 = lambda m:sum(sgn(p)*prod(m[k,p[k]-1]for k in range(m.ncols()))for p in Permutations(m.ncols()))

This directly implements Leibniz's formula. For a much better explanation of the formula and why it works than I could possibly write, see this excellent answer.

Real Part of e^(Tr(ln(M))), 48 bytes
lambda m:real(exp(sum(map(ln,m.eigenvalues()))))

This function returns symbolic expressions. To get a numerical approximation, call n(result) before printing.
This is an approach that I haven't seen anyone use yet. I'm going to give a longer, more-detailed explanation for this one.
Let A be a square matrix over the reals. By definition, the determinant of A is equal to the product of the eigenvalues of A. The trace of A is equal to the sum of A's eigenvalues. For real numbers r_1 and r_2, exp(r_1) * exp(r_2) = exp(r_1 + r_2). Since the matrix exponential function is defined to be analogous to the scalar exponential function (especially in the previous identity), and the matrix exponential can be computed by diagonalizing the matrix and applying the scalar exponential function to the eigenvalues on the diagonal, we can say det(exp(A)) = exp(trace(A)) (the product of exp(λ) for each eigenvalue λ of A equals the sum of the eigenvalues of exp(A)). Thus, if we can find a matrix L such that exp(L) = A, we can compute det(A) = exp(trace(L)).
We can find such a matrix L by computing log(A). The matrix logarithm can be computed in the same way as the matrix exponential: form a square diagonal matrix by applying the scalar logarithm function to each eigenvalue of A (this is why we restriced A to the reals). Since we only care about the trace of L, we can skip the construction and just directly sum the exponentials of the eigenvalues. The eigenvalues can be complex, even if the matrix isn't over the complex ring, so we take the real part of the sum.

The last part is a fascinating idea, but the header and explanation don't match the code, which doesn't take a matrix logarithm. It's just real(prod(m.eigenvalues())) ungolfed.

Wolfram Language (Mathematica), 3 bytes
Det

Try it online!
Per the meta consensus, mainly upvote nontrivial solutions that take effort to write.

Python 3, 238 bytes, 227 bytes, 224 bytes, 216 bytes
from functools import*
from itertools import*
r=range;n=len;s=sum
f=lambda l:s(reduce(lambda p,m:p*m,[l[a][b]for a,b in zip(r(n(l)),j)])*(-1)**s(s(y<j[x]for y in j[x:])for x in r(n(l)))for j in permutations(r(n(l))))

Try it online!
My solution uses the definition of a determinant for calculations. Unfortunately, the complexity of this algorithm is n! and I can not show the passage of the last test, but in theory this is possible.

Pari/GP, 6 bytes
matdet

Try it online!

Pari/GP, 52 bytes
f(a)=if(#a,sum(n=1,#a,-(-1)^n*a[1,n]*f(a[^1,^n])),1)

Try it online!

Python 2 + numpy, 29 bytes
from numpy.linalg import*
det

Try it online!

Julia, 3 bytes
det

Try it online!

MY, 3 (4?) bytes
ω∥↵

Try it online!
Precision is a little weird (the second test case is 5.000000000000001), this can be fixed at a cost of one byte: ω∥⌊↵

Right argument, parallel, return... what?

Right argument, determinant, output with new line.

TI-Basic (83 series), 137 133 bytes
Ans→[A]
For(N,max(dim(Ans)),2,-1
 Matr▶list(Ans,N,C
 1+sum(not(cumSum(not(not(∟C→R
 If R>N
 Then
  0[A]→[A]
 Else
  rowSwap([A],R,N
  *row(cos(π(R≠N))Ans(N,N),Ans,1
  *row(-Ans(N,N)⁻¹,Ans,N
  For(I,1,N
   *row+(Ans(I,N),Ans,N,I→[A]
  End
 End
End
Ans(1,1

(Indentation is purely decorative and doesn't appear on the calculator.)
Instead of the built-in det(, does row reduction by hand. The pivots are built backwards, starting at the (n,n)th entry and going back to the (1,1)th (or is that (1,1)st?).
Every time we need to divide by a number to set the pivot to 1, we multiply the first row by the same number (possibly negated, if we needed to swap rows), so that the determinant of the matrix stays the same. (If we didn't find a pivot, we zero out the entire matrix.) As a result, at the end, the (1,1)th (1,1)st top left entry holds the determinant.

Python 3, 153 149 Bytes
def d(n):
 l=len(n);r=range(0,l)
 if l==1:return n[0][0]
 else:return sum([(-1)**i*n[0][i]*d([[n[a][b]for b in r if b!=i]for a in r[1:]])for i in r])

I'm new to code golf, so I wouldn't be surprised if this could be improved.
It's calculated in a pretty basic way; by summing the product of each element of the first row with the determinant of its respective minor, with a sign change based on the element's position.
Edit: Realized I could shorten it by using a single character for the function name :/

Suggest range(l) instead of range(0,l) and replace the if l==1: with if l>1: and swap the resulting expressions.

Clojure, 141 bytes
(fn D[A](apply +(map *(cycle[1 -1])(first A)(if(next A)(for[R[(range(count A))]i R](D(for[a(rest A)](for[j(remove #{i}R)](nth a j)))))[1]))))

Laplace's formula, longer than I'd like...

Jelly, 2 bytes
ÆḊ

Try it online!
Currently only built-in solution. :( Looks like others have managed to find out interesting non-built-in solutions.

Excel VBA, 58 50 Bytes
Returns the determinant of the matrix entered into the the upper left corner of the default Sheet1 object.
?Evaluate("MDeterm("+Sheet1.UsedRange.Address+")")

-8 Bytes thanks to Greedo

Save 7 with ?Evaluate("MDeterm("&Sheet1.UsedRange.Address &")")

J, 5 bytes
-/ .*

Try it online!

APL(NARS), 64 chars, 128 bytes
{1≥k←↑⍴w←⍵:⍵⋄2=k:-/⌽w[2;]×⌽w[1;]⋄v←⍳k⋄-/w[1;]×∇¨{w[v∼1;v∼⍵]}¨⍳k}

This would be the formula used in the schools, it is easy, much easier than i tought...
if w is a 2x2 matrix -/⌽w[2;]×⌽w[1;] rotate first row multiply with second row than reverse
and make difference in pratice it does: (w[1;1]×w[2;2])-w[1;2]×w[2;1]
else if it is a matrix nxn with n>2 multiply the first row with its adjunt Matrix (if i remember well the name), and make -/;
test:
  D←{1≥k←↑⍴w←⍵:⍵⋄2=k:-/⌽w[2;]×⌽w[1;]⋄v←⍳k⋄-/w[1;]×∇¨{w[v∼1;v∼⍵]}¨⍳k}
  D ⊃(1 9 3)(4 ¯5  6)(7 8 9)
162
  D 3 3 ⍴ ⍳20
0
  D ⊃(1 0)(0 1)
1
  D ,89
89
  D  1
1
  D 'ac'
RANK ERROR

05AB1E, 23 bytes
ā<œε©2.ÆíÆ.±Iε®Nèè}«P}O

Port of @LeakyNun's Jelly answer.
Try it online or verify all but one test cases (the biggest test case times out on TIO).
Explanation:
ā             # Push a list in the range [1, (implicit) input-length]
 <            # Decrease it by 1 to make the range [0, input-length)
  œ           # Get all permutations of this list
   ε          # Map each permutation to:
    ©         #  Store the permutation in variable `®` (without popping)
     2.Æ      #  Get all 2-element combinations of this list
        í     #  Reverse each inner list
         Æ    #  Reduce each by subtraction
          .±  #  And take their signum (-1 if a<0; 0 if a==0; 1 if a>0)
    I         #  Push the input
     ε        #  Map each row to:
      ®       #   Get the current permutation from variable `®`
       Nè     #   Use the map-index to index into this permutation
         è    #   And use that to index into the row
     }«       #  After the map: merge the two lists together
       P      #  And take its product
   }O         # After the outer map: sum all value
              # (after which it is output implicitly as result)

Maple, 11 bytes
Documentaiton
Determinant

LinearAlgebra[Determinant](Matrix([[2, 3, 1], [3, 2, 3], [0, 3, 2]]));
LinearAlgebra[Determinant](LinearAlgebra[HilbertMatrix](4));


text:
The cyclic sequence of even digits, with odds in between
Consider the following sequence:
1, 0, 1, 2, 4, 1, 6, 8, 0, 1, 2, 4, 6, 8, 1, 0, 2, 4, 6, 8, 1, 0, 2, 4, 6, 8, 0, 1, ...

The even digits start from 0 and are grouped into runs of increasing length. They are arranged cyclically, meaning that they are sorted in ascending order until 8 is reached, and then cycled back from 0. 1 separates the runs of even digits, and it also starts the sequence. Let's visualize how this sequence is formed:
                 1, 0, 1, 2, 4, 1, 6, 8, 0, 1, 2, 4, 6, 8, 1, 0, 2, 4, 6, 8, 1,  ...

                    -     ----     -------     ----------     -------------
run length:         1      2          3            4                5            ...
position of 1:   X     X        X           X              X                 X   ...
even sequence:      0,    2, 4,    6, 8, 0,    2, 4, 6, 8,    0, 2, 4, 6, 8      ...

Acceptable Input and Output methods:

Receive an integer N as input and output the Nth term of this sequence.
Receive an integer N as input and output the first N terms of this sequence.
Print the sequence indefinitely.

You can choose either 0 or 1-indexing for the first two methods.
You can compete in any programming language, while using the standard input and output methods. Standard loopholes are forbidden. This is code-golf, so the shortest code in each language wins.

This challenge has been sandboxed.

Haskell, 50 46 bytes
1#cycle[0,2..8]
n#r=1:take n r++(n+1)#drop n r

Try it online!
1#cycle[0,2..8] returns the sequence as infinite list.
-4 bytes thanks to Ørjan Johansen!

You can save four bytes by making the first argument of # just a number. Try it online!

Husk, 12 11 10 bytes
ṁ:1CN¢mDŀ5

Try it online!
Prints the sequence indefinitely.
Alternatively:
J1CΘN¢mDŀ5

Try it online!
Explanation
        ŀ5   Start from [0, 1, 2, 3, 4]
      mD     Double each value to get [0, 2, 4, 6, 8]
     ¢       Repeat this list indefinitely, [0, 2, 4, 6, 8, 0, 2, ...]
   CN        Cut it into chunks of increasing lengths,
             [[0], [2, 4], [6, 8, 0], ...]
ṁ:1          Prepend 1 to each sublist and concate the resulting lists.

For the alternative solution:
     ¢mDŀ5   Again, get [0, 2, 4, 6, 8, 0, 2, ...].
  CΘN        This time we prepend a zero to the natural numbers, which
             prepends an empty list to the resulting chunks.
J1           Join all the sublists with 1.

We could also do ...ΘCN..., because Θ does "prepend default element", which prepends a zero for lists of integers and an empty list for lists of lists.

Jelly, 10 bytes
5ḶḤṁR€1pFḣ

Returns the first n items of the sequence.
Try it online!
How it works
5ḶḤṁR€1pFḣ  Main libk. Argument: n

5           Set the return value to 5.
 Ḷ          Unlength; yield [0, 1, 2, 3, 4].
  Ḥ         Unhalve; yield [0, 2, 4, 6, 8].
    R€      Range each; yield [[1], [1, 2], [1, 2, 3], ..., [1, ..., n]].
   ṁ        Mold; in the result to the left, replace [1] with [0], [1, 2] with
            [2, 4], [1, 2, 3] with [6, 8, 0], and so forth.
      1p    Take the Cartesian product of [1] and the result.
        F   Flatten the result.
         ḣ  Head; take the first n items of the result.

ಠ_ಠ Unhalve... Isn't that just Double?

It's just a mnemonic. H is halve, so Ḥ is unhalve. ÆA is arccosine, so ÆẠ is unarccosine.

Python 3, 52 bytes
def f(n):t=(8*n+1)**.5+1;return 0==t%1or(n-t//2)%5*2

Takes a 1-based index and returns True or an integral float.
Try it online!

APL (Dyalog Unicode), 52 59 56 bytes
r←c k
r←~n←0
:For j :In ⍳k
n+←j-1
r,←1,⍨j⍴n⌽0,2×⍳4
:End
r←k⍴r

Try it online!
This is a tradfn (traditional function) taking one argument k and returning the first k items of the sequence.
Thanks to @GalenIvanov for pointing out an error in the function.
Thanks to @Adám for 3 bytes.
How it works:
r←c k             ⍝ The function c takes the argument k and results in r
r n←1 0           ⍝ Initializing the variables r and n, and setting them to 1 and 0, respectively.
:For j :In ⍳k     ⍝ For loop. ⍳k yields [1, 2, 3, ..., k], and j is the control variable.
n+←j-1            ⍝ Accumulates j-1 into n, so it follows the progression (0, 1, 3, 6, 10, 15...)
r,←1,⍨j⍴n⌽0,2×⍳4  ⍝ This line is explained below.
:End              ⍝ Ends the loop
r←k⍴r             ⍝ return the first k items of r.
                  ⍝ ⍴ actually reshapes the vector r to the shape of k;
                  ⍝ since k is a scalar, ⍴ reshapes r to a vector with k items.
            2×⍳4  ⍝ APL is 1-indexed by default, so this yields the vector 2 4 6 8
          0,      ⍝ Prepend a 0 to it. We now have 0 2 4 6 8
        n⌽        ⍝ Rotate the vector n times to the left.
      j⍴          ⍝ Reshape it to have j items, which cycles the vector.
   1,⍨            ⍝ Append a 1, then
r,←               ⍝ Append everything to r.

Below are a Dfn(direct function) and a tacit function that also solve the challenge, both kindly provided by @Adám.

Dfn: {⍵⍴1,∊1,⍨¨j⍴¨(+\¯1+j←⍳⍵)⌽¨⊂0,2×⍳4} Try it online!
Tacit: ⊢⍴1,∘∊1,⍨¨⍳⍴¨(⊂0,2×⍳4)⌽⍨¨(+\¯1+⍳) Try it online!

I'm interested in the explanation of the tacit vesrion. Thanks!

@GalenIvanov I'll add one later today.

I noticed that the 3 functions yield incorrect answers - the subsequences always start from 0 after 1 - they should continue from the last even digit of the previous subsequence.

@GalenIvanov You're right. I'll see if I can fix it and add the explanations today.

Python 2, 51 bytes
k=n=0
while 1:print 1;k+=1;exec'print n%10;n+=2;'*k

Prints the whole sequence.
Try it online!

APL, 25 bytes
Returns nth term.
1,(⌽n↑⌽(+/⍳n←⎕)⍴0,2×⍳4),1

Explanation
n←⎕     Prompts for screen input of integer
+/⍳      Creates a vector of integers of 1 to n and sums
⍴0,2×⍳4  Creates a vector by replicating 0 2 4 6 8 to the length of sum
⌽n↑⌽   Rotates the vector, selects first n elements and rotates result
        (selects last n elements}
1,...,1 Concatenates 1s in front and behind result

C (gcc), 84 bytes
i;j;f(){for(i=1;;i++){printf("%d ",1);for(j=0;j<i;)printf("%d ",(2*j+++i*i-i)%10);}}

Try it online!
A function (f()) that prints the sequence infinitely, separated by spaces.
i is the length of the current even-run.
j is the index in the current even-run
(2*j+++i*i-i)%10 gives the correct even number, given i and j (and increments j), equivalent to ((j+Tr(i))%5)*2, where Tr(x) is the xth triangular number (which is the number of even numbers that have been printed before the current even run;

77 bytes

Java 8, 96 bytes
v->{for(int i=0,j=-2,k;;i++,System.out.println(1))for(k=0;k++<i;System.out.println(j%=10))j+=2;}

Prints indefinitely, each number on a new line.
Explanation:
Try it here.
v->{                               // Method with empty unused parameter and no return-type
  for(int i=0,                     //  Index integer, starting at 1
          j=-2,                    //  Even index integer, starting at -2
          k;                       //  Inner index integer
      ;                            //  Loop (1) indefinitely
                                   //    After every iteration:
       i++,                        //     Increase index `i` by 1
       System.out.println(1))      //     And print a 1
    for(k=0;                       //   Reset index `k` to 0
        k++<i;                     //   Inner loop (2) from 0 to `i`
                                   //     After every iteration:
       System.out.println(j%=10))  //      Set `j` to `j` modulo-10, and print it
      j+=2;                        //    Increase `j` by 2
                                   //   End of inner loop (2) (implicit / single-line body)
                                   //  End of loop (1) (implicit / single-line body)
}                                  // End of method

Batch, 85 bytes
@set/an=%2+1,t=n*-~n/2
@if %t% lss %1 %0 %1 %n%
@cmd/cset/a(%1-n)%%5*2*!!(t-=%1)+!t

Outputs the Nth term of the sequence. Works by calculating the next triangular number.

J, 46 42 40 bytes
-6 bytes thanks to cole
{.({.[:;[:#:&.>2^i.);@(1,&.><;.1)10|2*i.

Outputs the first N terms of this sequence.
How it works:
10|[:+:i. - generates a list of length N of 0 2 4 6 8 0 2 4... It simply takes  mod 10 of the doubled items of a list of integers starting at 0.
[:;[:#:&.>2^i. - generates a bit mask for cutting the above list.
(1 means start):  1 1 0 1 0 0 1 0 0 0 1 0 0 0 0 1 0 0 0 0 ... It finds 2 to the power of the consecutive non-negative integers, converts them into binary, flattens the list and takes only the first N items, so that the length of the both lists is the same.
;@(1,&.><;.1) - splits (cuts) the list of the even digits into sublistsaccording to the map of ones and zeroes, appends the sublist to 1 and finally flattens the resulting list
]{.  - takes only the first N items, getting rid of the additional numbers in the list due to the added 1s.
Try it online!

42 bytes: {.({.[:;[:#:&.>2^i.);@(1,&.><;.1)(10|2*i.). The changes I did were to use hooks and refactor the right tine of the fork to take advantage of how forks work. I like the 2^i. trick.
I'm trying to work on the left tine of the fork now.

@cole Thanks, I need to learn to use forks better. Apparently the fork 2*i. is better then than the capped hook  [:+:i.

You can also drop the parens on the right tine (10|2*i.) ->10|2*i.

Perl 5, 44 bytes
do{say 1;map{say(($r+=2)%10)}0..$i++}while 1

Try it online!
Infinite output

-4 bytes {...;redo} instead of do{...}while 1

Perl 5, 35 bytes
{say$i--?$r++%5*2:($i=++$j)>0;redo}

try it online

Common Lisp, 74 bytes
(do((b 1(1+ b))(k -2))(())(print 1)(dotimes(m b)(print(mod(incf k 2)10))))

Try it online!
Prints the sequence indefinitely.

JavaScript (ES6), 62 54 52 bytes
Returns the Nth term of the sequence, 0-indexed.
n=>(g=e=>n--?g(e+=k++<l?2:k=!++l):k<l?e%10:1)(k=l=0)

Demo

let f =

n=>(g=e=>n--?g(e+=k++<l?2:k=!++l):k<l?e%10:1)(k=l=0)

for(n = 0; n < 30; n++) {
  console.log('a(' + n + ') = ' + f(n))
}

Java (OpenJDK 8), 107 bytes
i->{String b="";for(int l=i,r=0,z=0;i>0;b+=l-i--==r?(l=i)<1+r++*0?"1":"1":"02468".charAt(z++%5));return b;}

Try it online!

Mathematica, 68 bytes
Returns the Nth term
Insert[Join@@Table[{0,2,4,6,8},#^2],1,Array[{(2-#+#^2)/2}&,#]][[#]]&

Try it online!

Wolfram Language (Mathematica), 61 bytes
Flatten[TakeList[2Mod[r=0~Range~#,5],UpTo/@r]~Riffle~1][[#]]&

Try it online!

bash, 42 bytes
for((;;)){ echo $[i--?r++%5*2:(i=++j)>0];}

or 34 if valid
echo $[i--?r++%5*2:(i=++j)>0];. $0

try it online

JavaScript, 45 bytes
1 indexed:
f=(x,t=0,p=0)=>p<x?f(x-1,t+1,p+t):x-p?x%5*2:1

0 indexed:
g=(x,p=0)=>p<x?g(x-1,p?++t+p:t=1):x-p?x%5*2:1

    f=(x,t=0,p=0)=>p<x?f(x-1,t+1,p+t):x-p?x%5*2:1

g=(x,p=0)=>p<x?g(x-1,p?++t+p:t=1):x-p?x%5*2:1

a.value=[...Array(100)].map((_,i)=>f(i+1))
b.value=[...Array(100)].map((_,i)=>g(i))
<output id=a></output>
<output id=b></output>

Ruby, 48 46 bytes
a=b=c=0;loop{a>b&&(p 2*b%10;b+=1)||a+=c+=p(1)}

Try it online!
Print the sequence indefinitely

Jelly, 17 bytes
9Ḷm2ẋ⁸ṁJR$$j1 1;ḣ

Try it online!

Proton, 55 bytes
i=0 j=-2while1{for k:0..i print(j=(j+2)%10)print(1)i++}

Try it online!
Prints the sequence indefinitely

ಠ_ಠ What is that, new syntax? :p i=0 j=-2while1{...

@Mr.Xcoder :P -2while is just like Python and while1 is because I made identifiers not able to be a keyword followed by a number

Pip, 15 bytes
W P1L UiP2*Uv%t

Prints indefinitely. Try it here! Or, here's a 16-byte equivalent in Pip Classic: Try it online!
Explanation
W                While this expression is truthy--
  P1             print 1 (evaluates to 1, which is truthy)
                 --loop:
      Ui          Increment i (initially 0)
    L             and loop that many times:
           Uv      Increment v (initially -1)
         2*        Multiply by 2
             %t    Mod 10
        P          Print

source: Stack Exchange

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

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

sort: votes


text:
The lowest initial numbers in a Fibonacci-like sequence
Given a positive integer input N, output the two non-negative numbers, a and b, where a < b, with the lowest possible mean value that will result in the number N being part of the recurring relation sequence:
f(0) = a
f(1) = b
f(n) = f(n-2)+f(n-1)

In case there are more than one solution where the mean of a and b are minimal, then you should output the one with the lowest b. Proof that this is not possible is provided by cardboard_box.
Your can assume N is within the representative range of integers in your language / system.
Test cases
N = 1
a = 0, b = 1

N = 15
a = 0, b = 3

N = 21
a = 0, b = 1

N = 27
a = 0, b = 9   <- Tricky test case. [3, 7] is not optimal and [4, 3] is not valid

N = 100
a = 4, b = 10

N = 101
a = 1, b = 12

N = 102
a = 0, b = 3

N = 1000
a = 2, b = 10

If a>=0 and a<b are there ever multiple solutions?

I can't guarantee that there are or aren't multiple solutions. Both 1,4 and 2,3 would give 5, and they have the same mean. I guess it's possible to find cases similar to that one, where these are the lowest mean values. If you can show/prove that there aren't multiple solutions then you don't need to check for that condition.

I checked up to 7457, and I didn't find any multiple solutions, although I haven't thought of any obvious proof.

totally inspired by https://codegolf.stackexchange.com/q/147200/67961

... But certainly not a dupe.

@ØrjanJohansen What do you means by "multiple solution"? Two pairs (a1, b1) and (a2, b2) such that a1+b1=a2+b2 and a1≠a2 and b1≠b2?

@user202729 Yes, but still with minimal sum.

The corresponding OEIS sequence for the lowest possible mean, A249783, has a wild looking graph.

@ØrjanJohansen I added a proof to my answer that there are no duplicate solutions (since my answer depended on it).

Husk, 19 18 16 14 13 15 bytes
Thanks Zgarb for saving 1 byte.
ḟö£⁰ƒẊ++ÖΣṖ2Θḣ⁰

Try it online!
Explanation:
Disclaimer: I really don't understand the ȯƒẊ++ section of the code.
Edit: It appears to translate to the Haskell fix.(mapad2(+).).(++), where mapad2 is apply function to all adjacent pairs in a list. (Although, knowing Husk, in the context of this program it could mean something else)
            Θḣ⁰    Create the list [0..input]
          Ṗ2       Generate all possible sublists of length 2
        ÖΣ         Sort them on their sums
ḟ                  Find the first element that satisfies the following predicate.
    ƒẊ++             Given [a,b], magically generate the infinite Fibonacci-like
                     sequence from [a,b] without [a,b] at the start.
 ö£⁰                 Is the input in that list (given that it is in sorted order)?

Save a byte with ö

I'm sure I tried that...

APL (Dyalog), 75 71 64 59 53 48 44 43 bytes
2 bytes saved thanks to @Adám
12 bytes saved thanks to @ngn
⊃o/⍨k∊¨+\∘⌽⍣{k≤⊃⍺}¨o←a/⍨</¨a←,⍉|-\¨⍳2⍴1+k←⎕

Try it online!
Uses ⎕IO←0.
How? This had gone real nuts.
k←⎕ - assign input to k
⍳2⍴1+k←⎕ - Cartesian product of the range 0 to k with itself
|-\¨ - substract each right pair element from the left, and get absolute values
a←,⍉ - transpose, flatten and assign to a
o←a/⍨</¨a - keep only pairs where the left element is smaller than the right, and assign to o
o now contains list of all pairs with a < b, ordered by their arithematical mean
+\∘⌽⍣{k≤⊃⍺}¨o - for each pair in o, apply fibonacci (reverse the pair and cumsum) until k or higher term is reached
k∊¨ - then decide if k is this last term (meaning it is contained in the sequence)
o/⍨ - and keep pairs in o where the previous check applies
⊃ - return the first result.

JavaScript (Node.js), 92 90 89 91 83 82 bytes
-3 bytes -1 byte thanks to ThePirateBay
-8 -9 bytes thanks to Neil.
f=(n,a=1,b=0,c=(a,b)=>b<n?c(a+b,a):b>n)=>c(a,b)?b+2<a?f(n,a-1,b+1):f(n,b-~a):[b,a]

Try it online!
Note: this solution relies on the fact that there are never multiple minimal solutions.
Proof that there are never multiple solutions:
Let FIB(a,b,k) be the Fibonacci-like sequence starting with a,b:
FIB(a,b,0) = a
FIB(a,b,1) = b
FIB(a,b,k) = FIB(a,b,k-1) + FIB(a,b,k-2)

Lemma 1
The difference between Fibonacci-like sequences is itself Fibonacci-like, i.e. FIB(a1,b1,k) - FIB(a0,b0,k) = FIB(a1-a0,b1-b0,k).
The proof is left to reader.
Lemma 2
For n >= 5, a solution a,b exists satisfying a+b < n:
if n is even, FIB(0,n/2,3) = n
if n is odd, FIB(1,(n-1)/2,3) = n
Proof
Cases where n < 5 can be checked exhaustively.
Suppose we have two minimal solutions for n >= 5, a0,b0 and a1,b1 with a0 + b0 = a1 + b1 and a0 != a1.
Then there exist k0,k1 such that FIB(a0,b0,k0) = FIB(a1,b1,k1) = n.

Case 1: k0 = k1
WLOG assume b0 < b1 (and therefore a0 > a1)
Let DIFF(k) be the difference between the Fibonnaci-like sequences starting with a1,b1 and a0,b0:
DIFF(k) = FIB(a1,b1,k) - FIB(a0,b0,k) = FIB(a1-a0,b1-b0,k) (Lemma 1)
DIFF(0) = a1 - a0 < 0
DIFF(1) = b1 - b0 > 0
DIFF(2) = (a1+b1) - (a0+b0) = 0
DIFF(3) = DIFF(1) + DIFF(2) = DIFF(1) > 0
DIFF(4) = DIFF(2) + DIFF(3) = DIFF(3) > 0
Once a Fibonnaci-like sequence has 2 positive terms, all subsequent terms are positive.
Thus, the only time DIFF(k) = 0 is when k = 2, so the only choice for k0 = k1 is 2.
Therefore n = FIB(a0,b0,2) = a0 + b0 = a1 + b1
The minimality of these solutions contradicts Lemma 2.
Case 2: k0 != k1:
WLOG assume k0 < k1.
We have FIB(a1,b1,k1) = n
Let a2 = FIB(a1,b1,k1-k0)
Let b2 = FIB(a1,b1,k1-k0+1)
Then FIB(a2,b2,k0) = FIB(a1,b1,k1) = FIB(a0,b0,k0) (exercise for the reader)
Since FIB(a1,b1,k) is non-negative for k >= 0, it is also non-decreasing.
This gives us a2 >= b1 > a0 and b2 >= a1+b1 = a0+b0.
Then let DIFF(k) = FIB(a2,b2,k) - FIB(a0,b0,k) = FIB(a2-a0,b2-b0,k) (Lemma 1)
DIFF(0) = a2 - a0 > 0
DIFF(1) = b2 - b0 >= (a0 + b0) - b0 = a0 >= 0
DIFF(2) = DIFF(0) + DIFF(1) >= DIFF(0) > 0
DIFF(3) = DIFF(1) + DIFF(2) >= DIFF(2) > 0
Once again, DIFF has 2 positive terms and therefore all subsequent terms are positive.
Thus, the only time when it's possible that DIFF(k) = 0 is k = 1, so the only choice for k0 is 1.
FIB(a0,b0,1) = n
b0 = n
This contradicts Lemma 2.

89 bytes

77 bytes

@Neil That minimizes b instead of minimizing a+b, and thus your solution gives f(27) = [3,7] but the optimal solution is f(27)=[0,9].  After reverting the breaking changes we're down to 83 bytes.

I think you can save another byte by using b-~a instead of a+b+1.

There's a small error in your second case: a2 >= a1 + b1 isn't correct when k1-k0=1. Instead you can use a2 >= b1 > a0 and b2 >= a1+b1 = a0+b0, and then the rest follows.

@ØrjanJohansen Thanks. Should be fixed now.

Haskell, 76 72 74 bytes
EDIT:

-4 bytes: @H.PWiz suggested using / instead of div, although this requires using a fractional number type.
+2 bytes: Fixed a bug with Enum ranges by adding -1.

f takes a value of Double or Rational type and returns a tuple of same.
Double should suffice for all values that are not large enough to cause rounding errors, while Rational is theoretically unlimited.
f n|let a?b=b==n||b<n&&b?(a+b)=[(a,s-a)|s<-[1..],a<-[0..s/2-1],a?(s-a)]!!0

Try it online! (with H.PWiz's header adjustments to input/output Rationals in integer format)
How it works

? is a locally nested operator in the scope of f. a?b recursively steps through the Fibonacci-like sequence starting with a,b until b>=n, returning True iff it hits n exactly.
In the list comprehension:

s iterates through all numbers from 1 upwards, representing the sum of a and b.
a iterates through the numbers from 0 to s/2-1. (If s is odd the end of range rounds up.)
a?(s-a) tests if the sequence starting with a,s-a hits n. If so, the list comprehension includes the tuple (a,s-a). (That is, b=s-a, although it was too short to be worth naming.)
!!0 selects the first element (hit) in the comprehension.

Python 2, 127 109 107 bytes
-2 bytes thanks to ovs (changing and to *)
g=lambda x,a,b:a<=b<x and g(x,b,a+b)or b==x
f=lambda n,s=1,a=0:g(n,a,s-a)*(a,s-a)or f(n,s+(a==s),a%s+(a<s))

Try it online!
Any bonus points for n,a,s-a?
Explanation:

The first line declares a recursive lambda, g, which verifies whether a, b expanded as a Fibonacci sequence will reach x. It also checks that a <= b, one of the criteria of the question. (This would allow cases where a == b, but in such a case 0, a would have already been discovered and returned).

The chained inequality a<=b<x performs two handy tasks at once: verifying a <= b, and that b < x.
If b < x yields True, the function calls itself again with the next two numbers in the Fibonacci sequence: b, a+b. This means the function will keep working out new terms until...
If b < x yields False, then we have reached the point where we need to check if b==x. If so, this will return True, signifying that the initial pair a, b will eventually reach x. Otherwise, if b > x, the pair is invalid.

The second line declares another recursive lambda, f, which finds the solution for a given value n. It recursively tries new initial pairs, a, b, until g(n, a, b) yields True. This solution is then returned.

The function recursively counts up initial Fibonacci pairs using two variables, s (initially 1) and a (initially 0). On each iteration, a is incremented, and a, s-a is used as the first pair. However, when a hits s, it is reset back to 0, and s is incremented. This means the pairs are counted up in the following pattern:
s=1 (0, 1) (1, 0)
s=2 (0, 2) (1, 1)  (2, 0)
s=3 (0, 3) (1, 2), (2, 1), (3, 0)
 Obviously, this contains some invalid pairs, however these are eliminated instantly when passed to g (see first bullet point).
When values a and s are found such that g(n, a, s-a) == True, then this value is returned. As the possible solutions are counted up in order of 'size' (sorted by mean, then min value), the first solution found will always be the smallest, as the challenge requests.

R, 183 bytes 160 bytes
n=scan();e=expand.grid(n:0,n:0);e=e[e[,2]>e[,1],];r=e[mapply(h<-function(n,a,b,r=a+b)switch(sign(n-r)+2,F,T,h(n,b,r)),n,e[,1],e[,2]),];r[which.min(rowSums(r)),]

Try it online!
Thanks to Giuseppe for golfing off 23 bytes
Code explanation
n=scan()                        #STDIO input
e=expand.grid(n:0,n:0)          #full outer join of integer vector n to 0
e=e[e[,2]>e[,1],]               #filter so b > a
r=e[mapply(
  h<-function(n,a,b,r=a+b)switch(sign(n-r)+2,F,T,h(n,b,r)),
                                #create a named recursive function mid-call
                                #(requires using <- vs = to denote local variable creation
                                #rather than argument assignment
  n,e[,1],e[,2]),]              #map n, a and b to h() which returns a logical
                                #which is used to filter the possibilities
r[which.min(rowSums(r)),]       #calculate sum for each possibility,
                                #get index of the minimum and return
                                #because each possibility has 2 values, the mean and
                                #sum will sort identically.

160 bytes -- in general, you should save bytes wherever you can, so saving 4 bytes by removing nice naming is not only acceptable or encouraged but in some sense required by [tag:code-golf]. Even so, nice answer, +1.

GolfScript - 88 77 bytes
~:N[,{1+:a,{[.;a]}/}/][{[.~{.N<}{.@+}while\;N=]}/]{)1=\;},{(\;~+}$(\;);~~' '\

I did not check for multiple solutions, thanks to cardboard_box!
Explanation
~:N                           # Reads input
[,{1+:a,{[.;a]}/}/]           # Creates an array of pairs [a b]
[{[.~{.N<}{.@+}while\;N=]}/]  # Compute solutions
{)1=\;},         # Pairs that are not solutions are discarded
{(\;~+}$         # Sorts by mean
(\;);~~' '\      # Formats output

Try it online!

Jelly, 19 bytes
ṫ-Sṭµ¡³e
0rŒcÇÐfSÐṂ

Try it online!
-1 byte thanks to the proof by cardboard_box. In case it's disproved, you can append UṂṚ to the end of the second line for 22 bytes in total.

...a leading increment should resolve @StewieGriffin's observation.

I have a feeling you can drop the Ḣ

We only need to find the seed that makes the input, x, appear latest. If x were found at the third index for multiple then it works for 0,x and would therefore also work at either 1,(x-1)/2 (x odd) or 2,x/2-1 (x even) whereupon x would appear later in the result, so that wont happen. For a later collision the mean can only be the same if the third terms are the same too, but then one must have a lower difference between the initial terms (else they'd be the same) and hence will have x found at a later index. As such we can do ṫ-Sṭµ¡i³¶ḶŒcÇÐṀ saving four bytes.

...oops, plus the increment: ṫ-Sṭµ¡i³¶‘ḶŒcÇÐṀ

@StewieGriffin That test-case didn't exist when I answered :p

@JonathanAllan In case there are more than one solution then you should output the one with the lowest maximum value. That means I must add bytes! D: time to find a proof EDIT: looks like you've got something that doesn't need the proof...

The test case didn't exist, but the possibility of such an input did ;-)

@StewieGriffin And my tiredness also did ;-P

@JonathanAllan Looks like your version was invalid :/

Why so? Maybe deleted comments?

@JonathanAllan Your version outputs [[0, 9], [3, 7]] on input 27.

@cardboard_box thanks I see! Could add SÐṂ for 19 though I suppose (or remove the Ḣ from Erik's, due to your proof).

... in fact the Ḣ could stay and you could still have 19 with ...ÐfSÞḢ

@JonathanAllan I don't think that works (sort by sum and take the head? looks like that's getting the lowest minimum). In fact, I think I need a 22-byter where I replace Ḣ with UṂṚ.

True dat. (+6.)

Mathematica, 117 bytes
If[#==1,{0,1},#&@@SortBy[(S=Select)[S[Range[0,s=#]~Tuples~2,Less@@#&],!FreeQ[LinearRecurrence[{1,1},#,s],s]&],Mean]]&

Try it online!

Batch, 160 158 bytes
@set/aa=b=0
:g
@if %a% geq %b% set/ab-=~a,a=0
@set/ac=a,d=b
:l
@if %c% lss %1 set/ad+=c,c=d-c&goto l
@if %c% gtr %1 set/aa+=1,b-=1&goto g
@echo %a% %b%

This (also) gives 3 7 on input 27. The correct solution is 0 9.

@cardboard_box Still not seeing where the question requires that...

In the first sentence: "with the lowest possible mean value".

@cardboard_box Ah, sorry, that was too easy to overlook.

@cardboard_box OK should be fixed now.

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

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


text:
Codeium Golfide

In this challenge, your task is to take an anion and a cation, and output the chemical formula of the compound.
Input rules

Take in 2 strings (in any order) representing the anion and cation, e.g. F, NH_4, or Al.
To take in the charge of each ion, you can either have it as part of the string separated by a caret (e.g. F^-1) or take in additional numerical arguments.
Note: As long as your numeric input type is signed, then the anion's charge will be passed in as a negative number.
The symbols will always be real, and charges accurate.

Output rules

Use _ for subscripts: Fe2O3 would be Fe_2O_3.
Cation first: NaCl, not ClNa.
Neutral molecule: Li2O, not LiO or LiO-.
Lowest possible coefficients: Fe2O3, not Fe4O6.
No subscript ones: NaCl, not Na1Cl1.
No shifting: NH4OH, not NH5O.
Conditional parentheses:
Do not use parentheses on a single-atom ion: MgCl2, not Mg(Cl)2.
Do not use parentheses if there is only one of the ion per molecule: KClO3, not K(ClO3).
DO use parentheses if there are two or more of a polyatomic ion: Be3(PO4)2, not Be3PO42 or Be3P2O8.

Testcases
Input               Output
Fe^+3, O^-2         Fe_2O_3
Fe^+2, O^-2         FeO
H^+1, SO_4^-2       H_2SO_4
Al^+3, SO_4^-2      Al_2(SO_4)_3
NH_4^+1, SO_4^-2    (NH_4)_2SO_4
Hg_2^+2, PO_4^-3    (Hg_2)_3(PO_4)_2
NH_4^+1, OH^-1      NH_4OH
Hg_2^+2, O_2^-2     Hg_2O_2

This is code-golf, the shortest answer in bytes wins!

Recommended test case: Fe^+2, OH^-1: Fe(OH)_2 for a polyatomic ion with 1 of each element (OH^-1).

Recommended test case: C(NH_2)_3^+1, NO_3: C(NH_2)_3NO_3

@Adám second ion is charged: NO_3^-1. Also another test case should be the first one paired up with a ^-2, so it would make (C(NH_2)_3)_2.... Or a case where the ion that is needed more than once begins with a bracket.

@Heimdall The missing ^-1 is a typo. Yes, I wanted such a compound, but  didn't find any real-life examples.

@Adám Fe_4(Fe(CN)_6)_3 for Prussian blue.

Recommended test case: Fe^+3, Fe(CN)_6^-4: Fe_4(Fe(CN)_6)_3

This might be the funniest challenge title I've seen on this website so far, and that's saying a lot for a ppcg question. Thanks for the good laugh

@osuka_ Have you seen "what's the fastest way to kill my family" over on Arqade? Or my other PPCG challenge, "confirm the Illuminati"?

@StephenLeppik I meant on ppcg itself, maybe I should’ve used “community” instead of “website” - but ppcg is probably the funniest community on SE, so the meaning is mostly the same. I have seen the Illuminati question though, I got a hearty laugh out of that too

Can I assume the absolute value of a charge is no more than ... maybe 4 or 6?

@Heimdall ions do exist with charges of up to +9… but no more than that. Anions go up to -4 (VERY rarely) and cations go up to +9 (one element only).

Thanks. When you say one element only, I assume you mean multi-atom ions have smaller upper limit. So, can I safely assume down to -4 and up to +9?

@Heimdall yes. Also, what I really meant was "there's only one cation in total with +9". (And the only polyatomic cation I know is ammonium (NH4) which is +1.)

APL (Dyalog), 60 59 61 bytes
+2 since charges must be given signed.
Anonymous infix function. Takes list of ions (anion, cation) as left argument and list of corresponding charges as right argument.
{∊(⍺{⍵∧1<≢⍺∩⎕D,⎕A:1⌽')(',⍺⋄⍺}¨m),¨(m←s≠1)/¨'_',∘⍕¨s←⍵÷⍨∧/⍵}∘|

Try it online!
{…}∘| function where ⍺ is left argument and ⍵ is right argument's magnitude:
 ∧/⍵ LCM of the charges
 ⍵÷⍨ divide the charges by that
 s← store in s (for subscripts)
 '_',∘⍕¨ format (stringify) and prepend underbar to each
 (…)/ replicate each letter of each with the corresponding value from:
  s≠1 Is s different from 1? (gives 1 or 0)
  m← store in m (for multiple)
 (…),¨ prepend the following respectively to those:
  ⍺{…}¨m for each, call this function with ions and m as arguments:
   ⎕D,⎕A Digits followed by uppercase Alphabet
   ⍺∩ intersection of ion and that
   ≢ tally the number of characters in that
   1< Is one less than that? (i.e. do we have a multi-element ion?)
   ⍵∧ and do we need multiple of that ion?
   : if so, then:
    ')(',⍺ prepend the string to the ion
    1⌽ cyclically rotate one step to the left (puts ) on the right)
   ⋄ else
    ⍺ return the ion unmodified
 ∊ ϵnlist (flatten)

C, 208 205 175 169 bytes
argv[1] : cation
argv[2] : anion
Charges of ions are taken on stdin.
#define z(b)for(i=*++v;*++i>95;);printf(b>1?*i?"(%s)_%d":"%s_%d":"%s",*v,b);
main(c,v,d,e,g,h,i)char**v,*i;{scanf("%d%d",&c,&d);for(e=c,h=-d;g=h;e=g)h=e%g;z(-d/e)z(c/e)}

OP specified the ions can be given in any order.

Retina, 86 80 bytes
Thanks to Neil for saving 6 bytes.
^
(
\^.+
)_$&$*
(1+)(\1|¶.+)+_(\1)+$
$#3$2_$#2
_1$

m)T`()``.+\)$|\(.[a-z]?\)
¶

Try it online!
Input is linefeed-separated (test suite uses comma-separation for convenience).
Explanation
^
(

We start by prepending a ( to each molecule. The ^ matches on line beginnings because the m) towards the end of the program sets multiline mode for all preceding stages.
\^.+
)_$&$*

We replace the ^[-+]n part with )_, followed by n copies of 1 (i.e. we convert the charges to unary, dropping the signs).
(1+)(\1|¶.+)+_(\1)+$
$#3$2_$#2

This stage does three things: it divides both charges by their GCD, converts them back to decimal and swaps them. The GCD can be found quite easily in regex, by matching the longest 1+ that lets us match both charges using only the backreference \1. To divide by this, we make use of Retina's "capture count" feature, which tells us how often a group has been used. So $#2 is the first charge divided by the GCD and $#3 is the second charge divided by the GCD (both in decimal).
_1$

We remove _1s from the ends of both parts.
m)T`()``.+\)$|\(.[a-z]?\)

And we drop the parentheses from lines which end in a ) (i.e. those that had a _1 there), as well as lines that only contain a single atom.
¶

Finally, we concatenate the two molecules by dropping the linefeed.

Haskell, 101 97 bytes
(s#t)n m|let x#i|j<-div(lcm n m)i=l$x:[l(x:['(':x++")"|l x<':'])++'_':show j|j>1]=s#n++t#m
l=last

Try it online! Example usage: Fe^+3, O^-2 is taken as ("Fe"#"O")3 2.

RPL (HP48 S/SX), 294.5 bytes
Yes, ridiculously large submission, not sure how competitive it will be...
DIR
  M
    « S ROT S SWAP ABS 4 PICK ABS
      DUP2 WHILE DUP2 REPEAT MOD SWAP END DROP2
      SWAP OVER / 4 ROLL J 3 ROLLD / ROT J
      ROT 0 > IF THEN SWAP END + »
  S
    « DUP "^" POS DUP2 1 + OVER SIZE SUB OBJ
      3 ROLLD 1 - 1 SWAP SUB »
  J
    IF OVER 1 ==
    THEN SWAP DROP
    ELSE DUP SIZE DUP2 DUP SUB "Z" > - 1 >
      IF
      THEN "(" SWAP + ")" +
      END
      "_" + SWAP +
    END
END

3 routines packaged neatly in a directory. M is the main one. It expects 2 strings on the stack formatted as ions and pushes a molecule string onto the stack.
S splits the ion into charge as a number and the element formula as a string. For example, "PO_4^-3" would be taken from the stack and -3 and "PO_4" pushed onto the stack.
J joins the number of ions with the formula and decides whether to wrap the formula into brackets. The bit before ELSE deals with 1 ion, leaving the string as it is. For example, if 1 and "PO_4" are on the stack, they are replaced by "PO_4". 1 and "H" gives "H".
The rest deals with multiple ions; if it's a single atom it's not in brackets, otherwise it is. To decide whether it is, I check the length of the string and check if the last character is >"Z". Boolean expressions return 1 for true and 0 for false. By subtracting the result of this comparison from the length of the string, I get 1 or less when it's single atom, otherwise more: length 1 is a single atom; length 2 will have a letter as the last character; for a single atom it's lower case, so >"Z", making result 1, otherwise 2; length 3 or more means more than 1 atom and with 0 or 1 subtracted from the length the result will be at least 2. For example, 3 and "PO_4" gives "(PO_4)_3". 3 and "Al" gives "Al_3".
M first splits each ion using S. After the first line the level 5 of the stack (so the deepest buries object) contains charge of the second ion, level 4 second ion's formula, level 3 first ion's formula, level 2 absolute value of the first ion's charge and level 1 absolute value of the second ion's charge again. For example, if given ions on the stack are "Al^+3", "SO_4^-2", we get -2, "SO_4", "Al", 3, 2.
Second line calculates gcd of the 2 charges (leaving the charges intact).
Third line divides each charge by the gcd (to calculate the multiples) and joins it with the ion's formula using J. So we have two strings each with one given ion with charge removed (or a multiple of it) and a charge of the second one buried behind them. For example, -2, "Al_2", "(SO_4)_3" (-2 is the charge of SO_4).
Fourth line unburies the charge and if it's positive it swaps the two strings (so that the cation comes first) before joining them. So in example above, because it's negative, they are joined in the order as they are: "Al_2(SO_4)_3".

Python 3, 131 bytes
lambda c,a,C,A:g(c,A/gcd(C,A))+g(a,C/gcd(C,A))
from math import*
g=lambda s,n:[s,[s,'(%s)'%s][sum(c<'`'for c in s)>1]+'_%d'%n][n>1]

Try it online!
Python 2, 196 174 170 155 149 140 136 bytes
lambda c,a,C,A:g(c,A/gcd(C,A))+g(a,C/gcd(C,A))
from fractions import*
g=lambda s,n:[s,[s,'(%s)'%s][sum(c<'`'for c in s)>1]+'_'+`n`][n>1]

Try it online!

Python 3 for 131 bytes

Python 3, 129 bytes
lambda E,e,I,i,m=lambda e,i:(len(e)>2<=i)*"("+e+(len(e)>2<=i)*")"+"_%d"%i*(i>1):m(E,i/gcd(i,I))+m(e,I/gcd(i,I))
from math import*

Try it online!

If we need to handle negative charges of the anions, then 153 bytes:
lambda E,e,I,i,a=abs,m=lambda e,i:(len(e)>2<=i)*"("+e+(len(e)>2<=i)*")"+"_%d"%i*(i>1):m(E,a(i)/gcd(a(i),a(I)))+m(e,a(I)/gcd(a(i),a(I)))
from math import*

Try it online!

JavaScript, 145 bytes
(c,a,p,n,g=(a,b)=>b?g(b,a%b):a)=>`(${c})_${n/-g(p,-n)}(${a})_${p/g(p,-n)}`.replace(/\(((\w+)\)(?=_1)|([A-Z][a-z]*)\))/g,"$2$3").replace(/_1/g,"")

Takes arguments c is cation, a is anion, p is positive charge, n is negative charge.


text:
Largest Deep Sum
Given an n-dimensional array filled with numbers, output which p-dimensional array has the largest sum.
Input

An n-dimensional array of numbers, or your language's equivalent. You can assume this is a rectangular array, or not.
n, the number of dimensions in the array. (doesn't need to be called n)
p, the number of dimensions in the arrays you are adding up. (doesn't need to be called p)

Output
The respective indices of the p-dimensional array with the largest sum. You may have undefined behavior in the case of a tie. The dimensions of the array that you find should be the last p dimensions of the original array.
Test Cases
Input:
[[[4,2],
  [6,8]],

 [[0,0],
  [500,327]]], 3, 2

Output:
1

Input:
[[1,2],
 [4,3]], 2, 0

Since p=0, you're just finding the largest number in the array, which in this case is 4. The coordinates are:
1 0

Input:
[[1,2],
 [4,-1]], 2, 1

You might see the 1 and 4 in the same column, and that that would have the largest sum of 5. However, because those share the second index, not the first, they are ignored. Both of the rows in this array add to 3, resulting in a tie, so you can return whatever you want here.

Input:
[[1,2],
 [3,4]], 2, 2

Since p=n, there are no indices of the array with the largest sum, so the output would be empty.
Additional Rules

Standard loopholes are forbidden.
Output can be in any reasonable format.
This is code-golf, so shortest answer in bytes wins.

I think your first example has an extraneous set of [...] around it. It might be worth clarifying exactly which dimensions we're adding. From the spec I would have assumed that we'd need to look at any subset of p dimensions, but judging by the last test case, we're only looking at adding up trailing dimensions (because I can add p=1 leading dimension and get 1+4 = 5).

@MartinEnder thanks for noticing that. As for your other note, that might make for a good sequel :)

Can we take two inputs as the array and n - p? It's way better for some languages.

Will the input array always have the same size in all dimensions?

@LuisMendo the array can be rectangular if you want.

What is the definition (being used here) of a p-dimensional [sub-]array? Do we only consider those in the same orientation as the input perspective? (is this why [[1,2], [4,-1]], 2, 1 returns either [1,2] or [4,-1] rather than [1,4]?)

APL (Dyalog), 21 bytes
(⍴⊤⊃∘⍒∘,)+/,⍤⎕⊢⎕

First input is the array, second is p.
Try it online!
How?
+/,⍤⎕ - apply sum over sub arrays of rank p
⊃∘⍒∘, - flatten, grade down and take first
⍴⊤ - encode in the mixed base derived by the shape of the previous result

Wolfram Language (Mathematica), 36 bytes
#~Position~Max@#&@Total[#,{-#2,-1}]&

Try it online!
Returns an array containing all the indices where the maximum value can be found. So for normal inputs that's a singleton array, containing the list of indices to it. If there's a tie, you'll get all the lists of indices for where that value shows up.
The function doesn't actually read the n parameter, because it isn't needed, but you could pass it in as a third argument if you like (it would just be ignored anyway).
One catch here is that for some reason the dimension specification {0,-1} for Total never sums anything. I don't really understand why, but luckily that's just what we need for this challenge.

Python 3 + numpy, 92 bytes
lambda a,n,p:unravel_index(argmax(sum(a,(*-arange(n-p,n),))),a.shape[p:])
from numpy import*

Try it online!
Quite a bit longer than I would have hoped. Errors out on cases with no solution.

Python 3, 138 bytes
There is probably a better way of doing this
lambda x,n,p:f(x,n,p)[:0:-1]
f=lambda x,n,p:n>p and max([*f(l,n-1,p),i]for i,l in enumerate(x))or[sum(f(l,n-1,p)[0]for l in x)if n else x]

Try it online!


text:
Take a byte out of it!
Your task is to, given an unsigned integer n, find the largest number which can be created by removing a single byte (8 consecutive bits) of data.

Example
Given the number 7831, we first convert it to binary (removing any leading zeroes):
1111010010111

We then find the consecutive group of 8 bits which, when removed, will yield the largest new result. In this case, there are 3 solutions, shown below
1111010010111
  ^      ^
   ^      ^
    ^      ^

Removing this any of these yields 11111, which we then convert back to its decimal value 31 for the answer.

Test Cases
256        ->   1
999        ->   3
7831       ->   31
131585     ->   515
7854621    ->   31261
4294967295 ->   16777215 (if your language can handle 32 bit integers)

Rules

It is guaranteed that the bit length of n will be larger than 8.
Your solution should theoretically work for any bit length of n larger than 8, but in practice, needs only work for integers 255 < n < 216
Input/Output should be in decimal.
You may submit a full program or a function.
This is code-golf, so the shortest program (in bytes) wins!

I don't understand why people put exclamation points in challenge titles! I think it might be a character limit thing! Might be just so people notice the challenge though!

@Mendeleev It's an imperative sentence. Those are usually terminated with exclamation points. It's only correct punctuation, why does it upset you so?

@Mendeleev People often use an exclamation mark to indicate a joke. The OP is highlighting the fact that he's making a pun. F. Scott Fitzgerald didn't like it, but in this context, it seems fine to me. If it wasn't there, you'd probably get people complaining about his spelling.

@Mendeleev because it's a bad pun...

@bornfromanegg I feel like people would notice the pun

@FlipTack Yeah, worse than mine

Jelly, 6 bytes
BḄ-8ƤṀ

A monadic link taking a number and returning a number.
Try it online!
How?
Uses a nice quick, Ƥ, developed by miles...
BḄ-8ƤṀ - Link: number
B      - convert to a binary list
    Ƥ  - for loop over some slices to be determined...
  -8   - this is a negative nilad, therefore: use overlapping outfixes of length 8
       -   (exactly what the specification asks us to inspect)
 Ḅ     -   convert from a binary list to an integer (vectorises)
     Ṁ - maximum

_> ... Wow you beat my by 10 bytes

J, 12 bytes
[:>./8#.\.#:

Try it online!
          #:     to binary
     8  \.       remove consecutive groups of eight
      #.         convert each result to decimal
  >./            maximum
[:               do nothing, this lets me avoid parentheses

What nifty algorithm do you have in there? Could you add an explanation?

@Mr. Xcoder FrownyFrog converts the number to a list of binary digits ( #: ), then converts all the 8-outfixes, or the list with consecutive 8-bit infixes removed back to decimal number system ( 8#..) and finally takes the largest one. [: simply caps the the two previous verbs, by making >./ to be executed monadically (only with right argument)

You taught me about outfix today, thanks for that! It’s a shame it doesn’t appear to be shorter to use under-&.; this is the perfect kind of problem for it.

Python 2, 41 bytes
f=lambda n:n>>8and max(n>>8,2*f(n/2)+n%2)

Try it online!

JavaScript (ES6), 54 bytes

f=(n,v=n>>8,b=1,m=0)=>b>v?m:f(n,(v^n)&b^v,b+b,v>m?v:m)
<input type=number min=256 max=2147483647 oninput=o.textContent=f(this.value)><pre id=o>

Works up to 2**31-1. Because someone asked for a bit-twiddling answer...

Pyth, 21 bytes
L&K.>b8eS,K+*2y/b2%b2

This is a recursive function (must be called with y, or see the link).
Try it here!

Mathematica, 69  bytes
Max@Array[Drop[#&@@s,#;;#+7]~FromDigits~2&,Last[s=#~RealDigits~2]-7]&

Try it online!
This solution works for large numbers Try it online!
-3 bytes from KellyLowder

Save 3 more bytes: Max[c=#~RealDigits~2;Array[Drop[c[[1]],#;;#+7]~FromDigits~2&,Last@c-7]]&

@KellyLowder nice! I golfed your solution a bit more

Python 3, 68 66 60 bytes
 -2 bytes thanks to Mr. Xcoder!
 -4 bytes thanks to ovs!
 -2 bytes thanks to Erik the Outgolfer!
lambda n:max(n%2**i|n>>i+8<<i for i in range(len(bin(n))-9))

Try it online!

66 bytes. n.bit_length()-8 is equivalent to len(bin(n))-10.

62 bytes

61 bytes

Wolfram Language (Mathematica), 46 bytes
Floor@If[#<256,0,Max[#/256,2#0[#/2]+#~Mod~2]]&

Try it online!
This version can only handle inputs up to 2518-1, otherwise we run into Mathematica's stack size limit. (The bound may vary between Mathematica installations.) The second solution in this answer avoids that.
How it works
A recursive approach based on the following logic:

The maximal value should be 0 for any input less than 256, since taking a byte out of the number eats the whole number. This is our base case, which is why it's included even though the specs promise us we won't have to handle such inputs.
Otherwise, we take the Max of two options: eat the lowest byte (giving us the input divided by 256) or chop off the lowest bit, recurse on the remaining integer, and append the lowest bit back when we're done.

Wolfram Language (Mathematica), 55 bytes
Max@Table[Mod[#,m=2^k]+Floor[#/m/2^8]m,{k,0,Log2@#-8}]&

Try it online!
An alternate version that builds a table instead of recursion, so it works for numbers of any size that Mathematica can handle.

The recursion depth is exceeded for numbers bigger than 10^160 although mathematica can handle bigger numbers. But I guess OP is fine with it

Java (OpenJDK 8), 138 134 bytes
i->{int l=0,m=0,x;String y=i.toString(i,2);for(;l<y.length()-7;m=x>m?x:m)x=i.valueOf(y.substring(0,l)+y.substring(l+++8),2);return m;}

Try it online!

i.toBinaryString(i) can be i.toString(i,2).

Jelly, 12 bytes
BJṡ8ḟ@€ƊịBḄṀ

Try it online!
Saved bytes thanks to Jonathan Allan!

A "filter out sublists then index into" version comes in at 13 bytes

C (gcc), 91 bytes
j;m;t;f(x){for(j=m=0;t=x>>j+8;m<t?m=t:j++)t=t<<j|x%(1<<j);return m;}

-23 bytes from Colera Su
Supports up to 2**31-1
Try it online!
Starts with the low 8 bits (j=0), then goes up, changing output if the number with bits [j,j+8) cut out is bigger than our current, and continuing until x has no bits above j+8

Store x>>j+8 and x>>j+8<<j|x%(1<<j) in a variable (parentheses removed) will reduce it to 68 bytes.

Retina, 71 67 64 bytes
.+
$*
+`(1+)\1
$+0
01
1
.
$`_$'¶
_.{7}

A`_
O^`
1G`
+1`\B
:$`:
1

Try it online! Link only includes the faster test cases, so as not to unduly overload Dennis's server. Edit: Saved 3 bytes thanks to @MartinEnder. Explanation:
.+
$*
+`(1+)\1
$+0
01
1

Convert from decimal to binary.
.
$`_$'¶
_.{7}

A`_

Construct a list of strings obtained by deleting 8 consecutive digits in all possible ways.
O^`
1G`

Sort them in reverse order and take the first (largest).
+1`\B
:$`:
1

Convert back to decimal. (See @MartinEnder's explanation.)

I came up with this shorter binary to decimal conversion a while ago. I've explained how it works in this answer.

ReRegex, 294 275 bytes
Saved 19 bytes by using better 'function' definitions
I'd say this is pretty good for a Regex only language.
The base lib does allow for conversion between Unary and Decimal (Which is needed as the challenge spec explicitly states decimal), but does not support Binary; So I had to write that as part of the script adding 120 bytes to it.
#import base
b(\d*):(_*)\2_b/b1$1:$2b/b(\d*):(_+)\2b/b0$1:$2b/b(\d+):b/$1/b:b/0/B(_*):1/B$1$1_:/B(_*):0/B$1$1:/B(_*):B/$1/j(\d*),\1(\d)(\d{7})(\d*):/j$1$2,$1$2$3$4:,B:$1$4B/j(\d*),\1\d{0,7}:,?(.*)/,$2,/,((_+)_+),(\2),/,$1,/,(_+),(\1_*),/,$2,/^,(_*),$/d<$1>/j,b:u<(?#input)>b:

Try it online!
By Individual Regexes.
#import base
b(\d*):(_*)\2_b/b1$1:$2b/
b(\d*):(_+)\2b/b0$1:$2b/
b(\d+):b/$1/
b:b/0/
B(_*):1/B$1$1_:/
B(_*):0/B$1$1:/
B(_*):B/$1/
j(\d*),\1(\d)(\d{7})(\d*):/j$1$2,$1$2$3$4:,B:$1$4B/
j(\d*),\1\d{0,7}:,?(.*)/,$2,/
,((_+)_+),(\2),/,$1,/
,(_+),(\1_*),/,$2,/
^,(_*),$/d<$1>/
j,b:u<(?#input)>b:

Steps
Firstly, we import the 'base' library, which gives two regexes. One which converts u<numbers> into unary. And one which converts d<unary_underlines> back into decimal. This is because the challenge requires IO in base10.
Then we define a handful of regexes which convert unary into binary.
b(\d*):(_*)\2_b/b1$1:$2b/
b(\d*):(_+)\2b/b0$1:$2b/
b(\d+):b/$1/
b:b/0/

The first of these, b(\d*):(_*)\2_b/b1$1:$2b/ searches for b, optionally followed by some binary digits, then a :, Then any amount of underlines, followed by the exact same amount of underlines plus one, and finally another b.
We then replace that with b1 followed by the binary digits from before, :, and just the first half of the underscores, and finally the last b.
So this checks if the unary is not divisible by two, and if so, prepends 1 to it's binary digits, then divides it minus one by two.
The second one, b(\d*):(_+)\2b/b0$1:$2b/ is almost idendical, however does not check for an extra _, meaning it only matches if it is divisible by two, and in this case prepends a 0 instead.
The third one checks if we're out of unary digits, and if so, strips away the padding to just leave the binary digits.
The last one checks if there never was any binary digits supplied, and in that case just leaves 0.
The next group of Regexes we define are to convert binary back into unary, and are slightly more simple.
B(_*):1/B$1$1_:/
B(_*):0/B$1$1:/
B(_*):B/$1/

The first of this group, B(_*):1/B$1$1_:/, much like its antithesis, detects a B, followed by any amount of Unary digits, then :1. It doesn't check for the matching B in this case, as it is only searching for one digit at a time. If this is matched, it doubles the previously matched amount of unary digits and adds one, then removes the one.
The second, B(_*):0/B$1$1:/, is almost idendical to the first, except matches a 0 rather than a 1, and does not add an additional unary digit.
The last of these, B(_*):B/$1/, checks if there are no more binary digits, and if so unwraps the unary. Unlike its antithesis, this does not need a special 0 case.
Next we define the j regexes, which act as a splitting function.
j(\d*),\1(\d)(\d{7})(\d*):/j$1$2,$1$2$3$4:,B:$1$4B/
j(\d*),\1\d{0,7}:,?(.*)/,$2,/

The first, j(\d*),\1(\d)(\d{7})(\d*):/j$1$2,$1$2$3$4:,B:$1$4B/ does most of the heavy lifting. It searches for j, optionally followed by binary digits which are the "incrementer", then a comma followed by the incrementer then exactly 8 binary digits followed by the rest of the binary number, then a :. The first of the 8 digits is appended to the incrementer, thus incrementing it, then everything but those 8 digits from the binary input is appended after the : following a ,. So (If we were using 2 digits instead of 8) j,1001: would become j1:1001:,01 then j10:1001,01,11. Additionally, the appended array elements are wrapped in Bs, to convert them back to unary.
The other, j(\d*),\1\d{0,7}:,?(.*)/,$2,/ checks if there are less than 8 binary digits left to check after the incrementer, and if so, removes everything other than the array wrapped in ,s. Eg. ,_,___,
During and after the creation of the array we define the comparison regexes.
,((_+)_+),(\2),/,$1,/
,(_+),(\1_*),/,$2,/

The first of these, ,((_+)_+),(\2),/,$1,/ checks a comma followed by some amount of underscores, then some more, followed by a comma, then the first amount of underscores, than a comma. It then replaces it with the total amount of underscores in the first element surrounded by ,s.
The latter, ,(_+),(\1_*),/,$2,/, checks for a comma followed by some amount of underscores followed by another comma, then the same amount or more underscores, and a last comma. This will instead leave the right element.
Finally, when there is on element left thus matching ^,(_*),$, we remove the surrounding commas and convert back to decimal via d<>. Then no-more regexes can fire and the output is presented.
The input is initially placed into the template j,b:u<(?#input)>b:, which first converts the decimal input to unary, eg 5 -> j,b:_____b:, then the resulting unary to binary, j,101: Then splits the binary (which doesn't work for the example), gets the largest element, converts back to decimal, and done.

JavaScript (ES6), 94 91 bytes
-3 bytes thanks to Justin Mariner
f=(n,d='',c=n.toString(2).match(`(${d}).{8}(.*)`))=>c?Math.max('0b'+c[1]+c[2],f(n,d+'.')):0

Just throwing out a JavaScript string-based solution, but I'm hoping someone will post a separate bitwise-based solution so I might learn something.
My solution recursively grabs an 8-bit chunk from the string, taking the maximum value that is found.

f=(n,d='',c=n.toString(2).match(`(${d}).{8}(.*)`))=>c?Math.max('0b'+c[1]+c[2],f(n,d+'.')):0

console.log(f(256));
console.log(f(999));
console.log(f(7831));
console.log(f(131585));
console.log(f(7854621));
console.log(f(4294967295));

I think you can drop the +(...) that converts '0b'+c[1]+c[2] to a number, because Math.max already does that. Try it online! (spec for future reference)

@JustinMariner, sweet, thanks!

C# (.NET Core), 122+13=135 120+13=133 bytes
n=>{int m=0,i=0,t;for(var b=Convert.ToString(n,2);i<b.Length-7;m=t>m?t:m)t=Convert.ToInt32(b.Remove(i++,8),2);return m;}

Try it online!
+13 for using System;
I imagine there is a way of doing this with out using Convert. Either way, I'm sure this could be reduced.
Acknowledgements
-2 bytes thanks to Kevin Cruijssen
UnGolfed
n=>{
    int m=0,
        i=0,
        t;

    // convert n to a binary string,
    // go through removing each possible byte,
    // check if this is the biggest int so far
    for (var b=Convert.ToString(n,2); i<b.Length-7; m=t>m?t:m)
        t=Convert.ToInt32(b.Remove(i++,8),2); // remove 8 bits from position i, then convert from binary string to int

    return m;
}

You can save a byte by changing the while to for and putting the var b in it: for(var b=Convert.ToString(n,2);i<b.Length-7;)

And you can save 1 more byte by adding a new int variable ,t and not using Math.Max, but a manual check instead: n=>{int m=0,i=0,t;for(var b=Convert.ToString(n,2);i<b.Length-7;m=t>m?t:m)t=Convert.ToInt32(b.Remove(i++,8),2);return m;} (120 + 13 bytes)

MATL, 23 21 bytes
Bn8-:"GB[]@8:q+(XBvX>

Try it online!
B                       % Implicitly grab input, convert to binary
 n8-:                   % Create list of 1,2,... n-8, with n the size of the binary string
     "                  % Loop over this list
      GB                % Grab the input again, convert to binary once more.
        @8:q+           % Create indices of a slice of length 8
             [](        % Index into binary string, delete the slice
                XB    % Convert the remainder from binary to integer
                  vX> % Get the maximum number so far.

Sadly, Bn8-:8:!+q&) only produces the slices to be removed, not the remainder we'd like to keep.

This saves a couple of bytes: Bn8-:"GB[]@8:q+(XBvX> (assign [] with ( instead of using &), and replace &: by : and addition)

@LuisMendo Thanks. I misread the docs, which said somewhere that you can only use a single index for null assignments, but for a different situation.

PHP, 67+1 bytes
do$r=max($r,$argn&($x=2**$i++-1)|$z=$argn>>8&~$x);while($z);echo$r;

Run as pipe with -nR or try it online.

Pyth, 19 bytes
eSmi.>>.<.BQd8d2a6l

Alternative answer:
eSmi.<<8.>.BQdd2a6l

Explanation:
eSmi.>>.<.BQd8d2a6lQ | Implicit Q at the end, where Q = input
  m             a6lQ | Map the over [0, 1, 2, ... , floor(log base 2 of Q) - 7]
         .BQ         |  Convert Q to binary string
       .<   d        |  Cyclically rotate left by d
      >      8       |  Get string from position 8 to end.
    .>        d      |  Cyclically rotate right by d
   i           2     |  Convert from binary string to integer
eS                   | Find the last element of sorted list (maximum value)

The other answer uses a similar approach, except that it rotates right first, and gets all bits except the last 8.

Ruby, 55 bytes
->n{(8../$/=~"%b"%n).map{|r|n>>8&(a=-1<<r-8)|n&~a}.max}

Try it online!

Java (OpenJDK 8), 73 bytes
n->{int m=0,i=1;for(;i<n>>7;i*=2)m=Math.max(m,n&~-i|(n>>8&-i));return m;}

Try it online!

Suggest n>>8&-i instead of (n>>8&-i)

Octave, 81 80 bytes
@(x)max(bin2dec(dec2bin(x*(c=2.^(0:(b=nextpow2(x+1)-8))))(:,[1:b b+9:end]))'./c)

Try it online!
This is a different solution to my original attempt, saving a further 14 bytes.
The code is broken down as follows:
@(x)max(                                                                       )
        bin2dec(                                                          )'./c
                                                         (:,[1:b b+9:end])
                dec2bin(x*(                            ))
                           c=2.^(0:                   )
                                   (b=nextpow2(x+1)-8)

On the sixth line the number of groups is calculated by finding the exponent of the next power of two larger than the input (number of bits in input number), and subtracting 7 as we are removing 8 bits from each group - the resulting number is stored in b for later.
We then calculate an array of powers of two on the fifth line which is large enough for all the possible groups that can be removed.  We save this in variable c for later.
On the next line up (forth line), we multiply the input by the array of powers of two (essentially bit shifting each number up), and convert the result to binary. If we take the example of 7831, this results in a 2D array containing:
000001111010010111
000011110100101110
000111101001011100
001111010010111000
011110100101110000
111101001011100000

If we then chop out the centre 8 bits, that is equivalent to removing each of the groups of 8 bits. This is done by the third line.
The resulting array is converted back to decimal on the second line. We also have to divide by c to undo the scaling that was done to each group initially.
Finally on the first line an anonymous function is declared, and the maximum value from all groups is calculated.

Save 1 byte by using nextpow2(x+1) rather than nnz(bin2dec(x))

Original attempt - 120 98 95 bytes
@(x)max(bin2dec(reshape(repmat(a=(d=@dec2bin)(x)',1,b=nnz(a)-7)(d(255*2.^(0:b-1))'<49),[],b)'))

Try it online!
The code is split up as follows:
@(x)max(                                                                                )
        bin2dec(                                                                             )
                reshape(                                                              ,[],b)'
                        repmat(a=(d=@dec2bin)(x)',1,b=nnz(a)-7)(                     )
                                                                d(255*2.^(0:b-1))'<49

Basically it calculates a matrix containing the possible groups of values that can be removed, and then works out which ends up giving the largest number.
Working line by line, the fifth line calculates the groups that can be removed. For example, take 7831. This is a 13bit number, giving the groups:
1111100000000
1111000000001
1110000000011
1100000000111
1000000001111
0000000011111

The result of the fifth line is a 2D logical array which can be used for indexing.
The fourth line of the code converts the input to an array of the bits (represented as characters '0' and '1'), and then replicates it n-7 times (where n in number of bits) giving one line for each possible grouping. The group mask above is used to remove each of the possible groups.
On the third line, the result is then reshaped to undo the unwanted flattening that resulted from applying the group mask. The second line converts back to an array of resulting decimal numbers. And the first line defines the anonymous function to be the maximum value of array of possible groups.

Saved 22 bytes by generating the matrix of groups using maths.
Saved 3 bytes in conversion from binary strings to logical group mask.

Perl 5, 78 + 1 (-p) = 79 bytes
map{$\=$"if($"=oct$s=~s/.{$_}\K.{8}//r)>$\}2..length($s=sprintf'0b%b',$_)-10}{

Try it online!

Perl, 53 bytes
(the use 5.10.1 to bring perl to lanugage level 5.10.1 is free)
Give the input number on STDIN.
Will run out of memory for big numbers, but the 32-bit number in the input is not yet a  problem
#!/usr/bin/perl
use 5.10.1;
$_=sprintf"%b",<>;/.{8}(?{\$F[oct"0b$`$'"]})^/;say$#F