Just one number must be removed from the set of integers between 1 and 100 (inclusive) so that the product of the factorials of the remaining numbers is a perfect square.
What is the least number of integers that must be removed from the same set of numbers so that the product of the factorials of the remaining numbers is a perfect cube?
The least number to remove is $12$, and the optimal solution is unique: $$ \frac{\prod_{k=1}^{100} k!}{37!39!43!59!60!69!71!72!81!91!93!99!} =(2^{1321}3^{648}5^{314}7^{204}11^{115}13^{95}17^{69}19^{61}23^{48}29^{35}31^{32}37^{25}41^{22}43^{20}47^{17}53^{13}59^{11}61^{11}67^{9}71^{8}73^{8}79^{6}83^{5}89^{3}97^{1})^3$$
I solved the problem via integer linear programming, with a binary decision variable $x_k$ for each $k\in\{1,\dots,100\}$ to indicate whether $k!$ appears in the product. For each prime $p \le 97$, let nonnegative integer decision variable $y_p$ represent the exponent of $p$ in the base of the cube. Let $n_{kp}$ be the exponent of $p$ in $k!$. That is, $k!=\prod_p p^{n_{kp}}$. The problem is to maximize $\sum_{k=1}^{100} x_k$ subject to linear constraints $$\sum_{k=1}^{100} n_{kp} x_k = 3 y_p \quad \text{for all primes $p \le 97$} \tag1\label1$$
The resulting optimal objective value is $88$, and the factorials $k!$ to remove are the ones where $x_k=0$.
To instead solve the perfect square problem, just change the $3$ to a $2$ in \eqref{1}. The optimal solution is again unique: remove $50!$.
_Not antispoilered since you'd have to accidentally read it quite closely to get spoilered.
Let us see how far we can get with pencil and paper:
Spoiler: We'll whittle it down to 1216 4864 (I can't count) cases and then leave the rest to a computer.
Warmup: perfect square
$1!2!...99!100! = 1^{100} 2^{99} ... 99^2 100^1$ is a square times $2 \times 4 \times ... \times 98 \times 100 = 2^{50} 50!$ As $50!$ is not a perfect square (for example, it is divisible by $47$ but not $47^2$) it is necessary to remove at least one factorial $\ge 47!$ and sufficient to remove $50!$.
Perfect cube
The full product is a perfect cube times $1\times 3^2\times 4\times6^2 \times ... \times 99^2 \times 100$. Removing a factorial $k!$ reduces all the exponents for the numbers $1,2...,k-1,k$ by one modulo 3.
The primes between 50 and 100 in this formula are $53^0 \times 59^0 \times 61^1 \times 67^1 \times 71^0 \times 73^1 \times 79^1 \times 83^0 \times 89^0 \times 97^1$. The multiplicity of $p$ in this list is affected only by the factorials $p!,p+1!,p+2!...$ and only by one incidence per factorial. To make all the exponents 0 modulo 3 which is a necessary but not a sufficient condition we can directly see that
Now it is time to look at a few more primes: $p_i=37,41,43,47$, each of those occurs in 2 numbers within the range. As the sum of these two always is a multiple of 3 they originally all have the same total exponent $e_{1/3}=1$ modulo 3. For them to simultaneously go to zero the number of factorials removed between any two $p_j,p_k$ of them must cancel (modulo 3) with that between $2p_j,2p_k$. In addition to compensating for the differential within the group $p_i$ we must also balance the shared effect on the group $p_i$ as a whole by removing a suitable number of factorials between 47 and 53.
In particular, depending on whether $f_4<$ or $\ge 82$ there must be two factorials removed either both between 37 and 41 or both between 41 and 43. $f_1$ and $f_5$ both affect all 4 $p_i$ in the same way, moreover they cancel each other. Similarly, the pairs $59,60$ and $71,72$ cancel $e_{1/3}$ for all 4 $p_i$.
If $f_2,f_3$ are both $<94$ their shared effect on all $4_pi$ cancels that of $f4$ while the differential effect must be compensated for by removing a factorial between 43 and 47.
If $f_2,f_3$ are both $\ge 94$ they have no differential effect but their shared effect doubles and must be compensated for by removing a factorial between 47 and 53.
If $f2 \ge 94 > f_3$ we'd have to remove two between 43 and 47 and another two between 47 and 53.
This puts the minimum number of factorials we must remove at 12 in the first two cases and 15 in the last.
As we are confident of being able to do it with 12 we disregard the last case.
To further constrain let's look at the next group of primes, $29,31$, that occurs 3 times each. To cancel the differential effect of $59,60$ on this pair (without removing any new factorials) we must have one of $f2,f3$ between $87$ and $93$ i.e. the first case and $f_2=93$.
Summary so far:
At this point we have reduced the problem to just 4864 cases. We could further reduce this number by counting individual primes, for example, 19 and its multiples are affected only by our choices for $f_{11},f_{12}$. But that is not very interesting, so I was happy to have my computer loop through the cases and it confirmed @Rob's answer.
Python script (requires numpy):
import itertools as it
import numpy as np
p = {}
for j in range(2,101):
if any(x[j]for x in p.values()):
continue
p[j] = np.zeros(101,int)
J = j
while J<100:
p[j][::J] += 1
J *= j
primes = np.array([*p.keys()])
factorizer = np.array([*p.values()])
full_prod = np.roll(np.r_[:101][::-1],1)
fpfactorized = factorizer @ full_prod
defect = fpfactorized % 3
corrector = np.tril(np.r_[0,100*[1]]) @ factorizer.T
candidates = np.zeros(12,int)
count = 0
candidates[:5] = 59,60,71,72,93
for candidates[5:10] in it.product(*(range(i,i+4) for i in [43,67,79,89,97])):
for candidates[10:] in it.combinations(
[41,42] if candidates[7]==82 else [37,38,39,40],2):
count += 1
if all(sum(corrector[candidates],0) % 3 == defect):
print(np.sort(candidates))
print(count)
